histw(x,varargin)

function histw(x,varargin)
% HISTW      Winsorized histogram
% INPUTS  :  x   - vector, possibly including NaN and Inf values
%            par - (optional) parameter structure, with fields (each optional)
%            pmin - lower cut-off percentile (default: 1)
%            pmax - upper cut-off percentile (default: 99)
%            ymax - upper limit of y-axis (see Note 2)
%            xmin - lower limit of x-axis (all values of x shown by default)
%            xmax - upper limit of x-axis (all values of x shown by default)
%            dx   - bin size ((q(pmax)-q(pmin))/10 by default)
%            xmrk - location of vertical line (none drawn by default)
% NOTES   : 1. HISTW draws a histogram of finite  values of X, excluding (on
%              the graph, not in the frequency calculation) values below X's
%              pmin'th  percentile,  or above pmax'th  percentile, with only
%              finite values counted.
%           2. ymax is set to min(0.01 + 1.1*fmax,1) by default, where  fmax 
%              is the maximum 'bin relative frequency' (~ height of a histo-
%              gram bar).(HISTW displays  relative frequencies, whereas HIST
%              plots absolute frequencies, i.e. counts).
%           3. As with HIST, x may not be a vector, but  must have  only one 
%              non-singleton dimension. (E.g., if A is a matrix, 'histw(A(:,
%              1))' is valid).
% EXAMPLE : A histogram of a standard-normal sample,  with top and bottom 5%
%%          of values excluded. X-axis limits are set to [-2,2]; bin size is
%%          0.1; a vertical line is drawn at x = 0.    
% par.pmin =  5;
% par.pmax = 95;
% par.xmin = -2;
% par.xmax =  2;
% par.dx   = .1;
% par.xmrk =  0;
% x = randn(1000,1);
% histw(x,par)
% AUTHOR  : Dimitri Shvorob, [e-mail], 12/1/06
x = squeeze(x);
if ~isvector(x)
   error('Input array must be a vector')
end 
if nargin > 1
   par = varargin{1};
   if ~isstruct(par)
      error('Input of structure type expected')
   end
   if ~isset('xmin'), par.xmin = -Inf; end
   if ~isset('xmax'), par.xmax = Inf;  end
   if ~isset('pmin'), par.pmin = 1;    end
   if ~isset('pmax'), par.pmax = 99;   end
   if par.pmin < 0 || par.pmin > 100
      error('Lower cut-off percentile out of [0,100] range')
   end 
   if par.pmax < 0 || par.pmax > 100
      error('Upper cut-off percentile out of [0,100] range')
   end 
   if par.pmin > par.pmax
      error('Lower cut-off percentile exceeds upper cut-off percentile')
   end
   if par.xmin > par.xmax
      error('Lower x-axis limit exceeds upper x-axis limit')
   end
else
   par.xmin = -Inf;
   par.xmax = Inf;
   par.pmin = 0;
   par.pmax = 100;
end

try
   q = prctile(x,[par.pmin par.pmax]);
catch
   error('Invalid percentiles supplied')
end
x = x(isfinite(x));
n = length(x);
x = x(x >= q(1) & x <= q(2));

if (nargin > 0 && ~isset('dx'))||(nargin == 0)
    par.dx = (q(2) - q(1))/10;
    if par.dx == 0
       par.dx = 1;
    end   
end

xmin = par.xmin; if ~isfinite(xmin), xmin = q(1); end
xmax = par.xmax; if ~isfinite(xmax), xmax = q(2); end
e = (xmin - par.dx):par.dx:(xmax + par.dx);
f = histc(x,e)/n;
bar(e,f,'histc')
hold on

if (nargin > 0 && ~isset('ymax'))||(nargin == 0)
   par.ymax = min(.01 + 1.1*max(f),1);
end

if nargin > 0 && isset('xmrk')
   line([par.xmrk par.xmrk],[0 par.ymax])
end
hold off

axis([par.xmin par.xmax 0 par.ymax])

    function[y] = isset(field)
    if ~isfield(par,field)
        y = false;
    else
        f = par.(field);
        if ~isscalar(f)
           error([field ' must be a scalar'])
        end 
        if isempty(f)||~isfinite(f)
           y = false;
        else
           y = true;
        end 
    end  
    end

end