function [bc,fhat,bins]=norm_ref_rule_w_skew_hist(x)
% NORM_REF_RULE_W_SKEW_HIST - returns the histogram based PDF for the input data 
% using the normal reference rule for selcting the bin widths and the skew factor 
% which should be introduced for skewwed distributions.
%   
% Written by:
% -- 
% John L. Weatherwax                2008-02-20
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

% specify the bin width using the normal reference rule for histograms: 
%
s = std(x); 
n = length(x); 
h = 3.5*s*(length(x)).^(-1/3);

% compute the skew factor (for histograms): 
sf = ((2^(1/3))*s) / ( exp(5*s^2/4) * ((s^2 + 2)^(1/3)) * (exp(s^2) - 1)^(1/2) );
h  = h*sf; 

% get the limits, bins, bin centers etc:
x_lim_left = min(x)-1;
x_lim_rght = max(x)+1;
t0   = x_lim_left;
tm   = x_lim_rght; 
rng  = tm-t0; 
nbin = ceil(rng/h); 
bins = t0:h:(nbin*h+t0);        % <- the bin edges ... 
bc   = bins(1:end-1)+0.5*h;     % <- the bin centers ... 

x(find(x<x_lim_left))=x_lim_left;
x(find(x>x_lim_rght))=x_lim_rght;
vk=histc(x,bins); vk(end)=[];
% normalize: 
fhat = vk/(n*h);