function [kde,xgrid,mker] = kdeSM(data,paramstruct)
% KDESM, Kernel Density Estimate (1-d, Gaussian Kernel)
% Steve Marron's matlab function
% Does 1-d kernel density estimation, using binned (default) or
% direct (either matrix, or loops for bigger data sets),
% implementations, with the bandwidth either user specified
% (can be vector), or data driven (SJPI, Normal Reference,
% Silverman's ROT2, Oversmoothed).
%
% Inputs:
% data - either n x 1 column vector of 1-d data
% or vector of bincts, when imptyp = -1
%
% paramstruct - a Matlab structure of input parameters
% Use: "help struct" and "help datatypes" to
% learn about these.
% Create one, using commands of the form:
%
% paramstruct = struct('field1',values1,...
% 'field2',values2,...
% 'field3',values3) ;
%
% where any of the following can be used,
% these are optional, misspecified values
% revert to defaults
%
% fields values
%
% vh vector of bandwidths, or specifies data driven:
% 0 (or not specified) - Sheather Jones Plug In
% -1 - Simple Normal Reference
% -2 - Silverman's Rule Of Thumb 2
% (20% Smaller than min of sd and IQR)
% -3 - Oversmoothed
% Note: <0 only works for imptype = 0
% >0 - Use input number (numbers if vector)
% Note: this MUST be >0 for imptyp >= 1
% (the direct implementations)
%
% vxgrid vector of parameters for, or values of, grid to evaluate at:
% 0 (or not specified) - use endpts of data and 401 bins
% [le; lr] - le is left end, re is right, 401 bins
% (get error message and no return if le > lr)
% [le; lr; nb] - le left, re right, and nb bins
% xgrid - Use input values
% Note: need to have more than 3 entries,
% and only works when imptyp = 1 or 2
%
% imptyp flag indicating implementation type:
% -1 - binned version, and "data" is assumed to be
% bincounts of prebinned data
% 0 (or not specified) - linear binned version
% and bin data here
% 1 - Direct matrix implementation
% 2 - Slow looped implementation (only useful when
% 1 creates matrices that are too large)
%
% eptflag endpoint truncation flag (only has effect when imptyp = 0):
% 0 (or not specified) - move data outside range to
% nearest endpoint
% 1 - truncate data outside range
%
% ibdryadj index of boundary adjustment
% 0 (or not specified) - no adjustment
% 1 - mirror image adjustment
% 2 - circular design
%
% ndataoverlay overlay raw data as a jitterplot (requires imptyp >= 0)
% (only has effect when plot is made here)
% 0 - (default) no data plot
% 1 - overlay up to 1000 points (random choice, when more)
% 2 - overlay full data set
% n > 2 - overlay n random points
%
% dolcolor data overlay color
% (only has effect when plot is made here)
% default is 'g'
%
% ibigdot 0 (default) use Matlab default for dot sizes
% 1 force large dot size in prints (useful since some
% postscript graphics leave dots too small)
% (Caution: shows up as small in Matlab view)
%
% dolhtseed seed for random heights used in data overlay jitter plot
% default is [] (for using current Matlab seed)
% (should be an integer with <= 8 digits)
%
% linewidth width of lines (only has effect when plot is made here)
% default is 2, for length(vh) <= 3,
% default is 0.5, for length(vh) > 3,
%
% linecolor string with color for lines
% (only has effect when plot is made here)
% default is 'b'
% use the empty string, '', for standard Matlab colors
%
% titlestr string with title (only has effect when plot is made here)
% default is empty string, '', for no title
%
% titlefontsize font size for title
% (only has effect when plot is made here,
% and when the titlestr is nonempty)
% default is empty [], for Matlab default
%
% xlabelstr string with x axis label
% (only has effect when plot is made here)
% default is empty string, '', for no xlabel
%
% ylabelstr string with y axis label
% (only has effect when plot is made here)
% default is empty string, '', for no ylabel
%
% labelfontsize font size for axis labels
% (only has effect when plot is made here,
% and when a label str is nonempty)
% default is empty [], for Matlab default
%
% plotbottom bottom of plot window,
% use [] to get 0 - .05 * range (default)
%
% plottop top of plot window,
% use [] to get max + .05 * range (default)
%
% iplot 1 - plot even when there is numerical output
% 0 - (default) only plot when no numerical output
%
% Outputs:
% (none) - Draws a graph of the result (in the current axes)
% kde - col vector of heights of kernel kernel density estimate,
% unless vh is a vector (then have matrix, with
% corresponding cols as density estimates)
% xgrid - col vector grid of points at which estimate(s) are
% evaluted (useful for plotting), unless grid is input,
% can also get this from linspace(le,re,nb)'
% mker - matrix (vector) of kernel functions, evaluated at xgrid,
% which can be plotted to show "effective window
% sizes" (currently scaled to have mass 0.05, may
% need some vertical rescaling)
%
% Assumes path can find personal functions:
% vec2matSM.m
% lbinrSM.m
% bwsjpiSM.m
% bwosSM.m
% bwrotSM.m
% bwsnrSM.m

% Copyright (c) J. S. Marron 1996-2001

% First set all parameters to defaults
vh = 0 ; % use default SJPI
vxgrid = 0 ;
imptyp = 0 ;
eptflag = 0 ;
ibdryadj = 0 ;
ndataoverlay = 0 ;
dolcolor = 'g' ;
ibigdot = 0 ;
dolhtseed = [] ;
linewidth = 2 ;
linecolor = 'b' ;
titlestr = '' ;
titlefontsize = [] ;
xlabelstr = '' ;
ylabelstr = '' ;
labelfontsize = [] ;
plotbottom = [] ;
plottop = [] ;
iplot = 0 ;

% Now update parameters as specified,
% by parameter structure (if it is used)
%
if nargin > 1 ; % then paramstruct has been added

if isfield(paramstruct,'vh') ; % then change to input value
vh = getfield(paramstruct,'vh') ;
end ;

if isfield(paramstruct,'vxgrid') ; % then change to input value
vxgrid = getfield(paramstruct,'vxgrid') ;
end ;

if isfield(paramstruct,'imptyp') ; % then change to input value
imptyp = getfield(paramstruct,'imptyp') ;
end ;

if isfield(paramstruct,'eptflag') ; % then change to input value
eptflag = getfield(paramstruct,'eptflag') ;
end ;

if isfield(paramstruct,'ibdryadj') ; % then change to input value
ibdryadj = getfield(paramstruct,'ibdryadj') ;
end ;

if isfield(paramstruct,'ndataoverlay') ; % then change to input value
ndataoverlay = getfield(paramstruct,'ndataoverlay') ;
end ;

if isfield(paramstruct,'dolcolor') ; % then change to input value
dolcolor = getfield(paramstruct,'dolcolor') ;
end ;

if isfield(paramstruct,'ibigdot') ; % then change to input value
ibigdot = getfield(paramstruct,'ibigdot') ;
end ;

if isfield(paramstruct,'dolhtseed') ; % then change to input value
dolhtseed = getfield(paramstruct,'dolhtseed') ;
end ;

if isfield(paramstruct,'linewidth') ; % then change to input value
linewidth = getfield(paramstruct,'linewidth') ;
end ;

if isfield(paramstruct,'linecolor') ; % then change to input value
linecolor = getfield(paramstruct,'linecolor') ;
end ;

if isfield(paramstruct,'titlestr') ; % then change to input value
titlestr = getfield(paramstruct,'titlestr') ;
end ;

if isfield(paramstruct,'titlefontsize') ; % then change to input value
titlefontsize = getfield(paramstruct,'titlefontsize') ;
end ;

if isfield(paramstruct,'xlabelstr') ; % then change to input value
xlabelstr = getfield(paramstruct,'xlabelstr') ;
end ;

if isfield(paramstruct,'ylabelstr') ; % then change to input value
ylabelstr = getfield(paramstruct,'ylabelstr') ;
end ;

if isfield(paramstruct,'labelfontsize') ; % then change to input value
labelfontsize = getfield(paramstruct,'labelfontsize') ;
end ;

if isfield(paramstruct,'plotbottom') ; % then change to input value
plotbottom = getfield(paramstruct,'plotbottom') ;
end ;

if isfield(paramstruct,'plottop') ; % then change to input value
plottop = getfield(paramstruct,'plottop') ;
end ;

if isfield(paramstruct,'iplot') ; % then change to input value
iplot = getfield(paramstruct,'iplot') ;
end ;

end ; % of resetting of input parameters

% Calculate kde
%
if imptyp > 0 ; % Then do direct implementation

if min(vh) > 0 ; % Then have valid bandwidths, so proceed

n = length(data) ;

if length(vxgrid) > 3 ; % Then use input grid
xgrid = vxgrid ;
nbin = length(xgrid) ;
lend = min(xgrid) ;
rend = max(xgrid) ;
else ; % Need to generate a grid
nbin = 401 ; % Default
lend = min(data) ; % Default
rend = max(data) ; % Default
if length(vxgrid) >= 2 ; % use input endpoints
lend = vxgrid(1) ;
rend = vxgrid(2) ;
end ;
if length(vxgrid) == 3 ; % use number of grid points
nbin = vxgrid(3) ;
end ;

if lend > rend ; % Then bad range has been input
disp('!!! Error in kdeSM: invalid range chosen !!!') ;
xgrid = [] ;
else ;
xgrid = linspace(lend,rend,nbin)' ;
end ;
end ;

% do boundary adjustment if needed
%
if ibdryadj == 1 ; % then do mirror image adjustment
badata = [(lend - (data - lend)); data; (rend + (rend - data))] ;
elseif ibdryadj == 2 ; % then do circular design adjustment
badata = [(data - (rend - lend)); data; (rend - lend + data)] ;
else ;
badata = data ;
end ;

% Loop through bandwidths
kde = [] ;
for ih = 1:length(vh) ;
h = vh(ih) ;

if imptyp ~= 2 ; % Then do direct matrix implementation
kdeh = vec2matSM((badata ./ h),nbin) - vec2matSM((xgrid' ./ h),n) ;
% efficient way to divide all dif's by h
% variable name "kde" is used to avoid creating too many biggies
kdeh = exp(-(kdeh .^2) / 2) ;
% exponential part of Gaussian density
kdeh = sum(kdeh)' ;
% sum part of kde, and make result a column vector
kdeh = kdeh / (n * h * sqrt(2 * pi)) ;
% normalize, and mult by Gaussain density constant
kde = [kde kdeh] ;
else ; % Do slower looped implementation
kdeh = [] ;
for ixg = 1:nbin ; % Loop through grid points
kdehx = (badata - xgrid(ixg)) / h ;
kdehx = sum(exp(-(kdehx .^2) / 2)) ;
kdeh = [kdeh; kdehx] ;
end ;
kdeh = kdeh / (n * h * sqrt(2 * pi)) ;
kde = [kde kdeh] ;
end ;
end ;

else ; % Have invalid bandwidths

disp('!!! Error in kdeSM: A bandwidth is invalid !!!') ;
disp(' (Note: cannot use data driven, with direct impl''s)') ;

end ;

else ; % Then do binned implementation

if imptyp == -1 ; % Then data have already been binned

if (length(vxgrid) == 1) | (length(vxgrid) > 3) ;
% Then can't proceed because don't have bin ends
disp('!!! Error: kdeSM needs to know the endpoints !!!') ;
disp('!!! to use this implementation !!!') ;
bincts = [] ;
else ;
bincts = data ;

nbin = 401 ;
lend = vxgrid(1) ;
rend = vxgrid(2) ;
if length(vxgrid) == 3 ; % then use number of grid points
nbin = vxgrid(3) ;
end ;

if nbin ~= length(bincts) ; % Then something is wrong
disp('!!! Warning: kdeSM was told the wrong number of bins !!!') ;
disp('!!! will just use the number of counts. !!!') ;
nbin = size(bincts,1) ;
end ;
end ;

else ; % Then need to bin data

if length(vxgrid) > 3 ; % Then need to warn of change to default
disp('!!! Warning: kdeSM was given an xgrid, and also !!!') ;
disp('!!! asked to bin; will bin and ignore xgrid !!!') ;
end ;

% Specify grid parameters
nbin = 401 ; % Default
lend = min(data) ; % Default
rend = max(data) ; % Default
if (length(vxgrid) == 2) | (length(vxgrid) == 3) ;
% then use input endpoints
lend = vxgrid(1) ;
rend = vxgrid(2) ;
end ;
if length(vxgrid) == 3 ; % then use number of grid points
nbin = vxgrid(3) ;
end ;

if lend > rend ; % Then bad range has been input
disp('!!! Error in kdeSM: invalid range chosen !!!') ;
bincts = [] ;
else ;
bincts = lbinrSM(data,[lend,rend,nbin],eptflag) ;
end ;

% Can do data-based bandwidth selection here, if specified
if vh == -1 ; % Then use Simple Normal Reference
vh = bwsnrSM(data) ;
elseif vh == -2 ; % Then use Silverman's Rule Of Thumb 2
% (~10% Smaller than min of sd and IQR)
vh = bwrotSM(data) ;
elseif vh == -3 ; % Then use Terrell's Oversmoother
vh = bwosSM(data) ;
elseif min(vh) <= 0 ; % Then be sure to use default SJPI
% (in case an unsupported value was input)
vh = 0 ;
end ;

end ;
n = round(sum(bincts)) ;
% put this here in case of truncations during binning

% do boundary adjustment if needed
%
if ibdryadj == 1 ; % then do mirror image adjustment
babincts = [flipud(bincts); bincts; flipud(bincts)] ;
elseif ibdryadj == 2 ; % then do circular design adjustment
babincts = [bincts; bincts; bincts] ;
else ;
babincts = bincts ;
end ;

% Get bandwidth (if still not yet specified)
if vh == 0 ; % Then use SJPI bandwidth
vh = bwsjpiSM(bincts,[lend; rend; nbin],0,-1) ;
end ;

% Loop through bandwidths
kde = [] ;
for ih = 1:length(vh) ;
h = vh(ih) ;

% Create vector of kernel values, at equally spaced grid
delta = (rend - lend) / (nbin - 1) ; % binwidth
k = nbin - 1 ; % index of last nonzero entry of kernel vector
arg = linspace(0,k * delta / h,k + 1)' ;
kvec = exp(-(arg.^2) / 2) / sqrt(2 * pi) ;
kvec = [flipud(kvec(2:k+1)); kvec] ;

% Do actual kernel density estimation
kdeh = conv(babincts,kvec) ;

if ibdryadj == 1 | ibdryadj == 2 ; % then did boundary adjustment
kdeh = kdeh(nbin+k+1:k+2*nbin) / (n * h) ;
else ;
kdeh = kdeh(k+1:k+nbin) / (n * h) ;
end ;

if h < 3 * delta ; % Then need to normalize
% to make numerical integral roughly 1
kdeh = kdeh / (sum(kdeh) * delta) ;
end ;

kde = [kde kdeh] ;
end ;

xgrid = linspace(lend,rend,nbin)' ;

end ;

% Create matrix of kernels, if this is needed
%
if nargout == 3 ;
cent = mean([lend; rend]) ;
% centerpoint of evaluation grid
if length(vh) > 1 ;
mih = vec2matSM(1 ./ vh',nbin) ;
mker = vec2matSM((xgrid - cent),length(vh)) .* mih;
% argument of gaussian kernel
else ;
mih = 1 / vh ;
mker = (xgrid - cent) .* mih;
% argument of gaussian kernel
end ;
mker = exp(-mker.^2 / 2) .* mih / sqrt(2 * pi) ;
% Gaussian kernels with mass 1
mker = 0.05 * mker ;
% Make masses = 0.05
end ;

% Make plots if no numerical output requested, or if plot requested
%
if nargout == 0 | ...
iplot == 1 ; % Then make a plot

if length(vh) > 3 & ~isfield(paramstruct,'linewidth') ;
% then need to change default value of linewidth
linewidth = 0.5 ;
end ;

if isempty(linecolor) ;
plot(xgrid,kde,'LineWidth',linewidth) ;
else ;
plot(xgrid,kde,'LineWidth',linewidth,'Color',linecolor) ;
end ;

if isempty(plottop) & isempty(plotbottom) ; % then adjust top and bottom
plotbottom = 0 ;
plottop = max(max(kde)) ;
plotrange = plottop - plotbottom ;
plotbottom = plotbottom - 0.05 * plotrange ;
plottop = plottop + 0.05 * plotrange ;
elseif isempty(plottop) ; % then only adjust top
plottop = max(max(kde)) ;
plotrange = plottop - plotbottom ;
plottop = plottop + 0.05 * plotrange ;
elseif isempty(plotbottom) ; % then only adjust bottom
plotbottom = 0 ;
plotrange = plottop - plotbottom ;
plotbottom = plotbottom - 0.05 * plotrange ;
end ;

vax = [lend,rend,plotbottom,plottop] ;
axis(vax) ;

if ndataoverlay > 0 & imptyp >= 0 ;
% then overlay jitterplot of data

if ~isempty(dolhtseed) ;
rand('seed',dolhtseed) ;
end ;

if ndataoverlay == 1 ;
ndo = min(n,1000) ;
elseif ndataoverlay == 2 ;
ndo = n ;
else ;
ndo = min(n,ndataoverlay) ;
end ;

flagleft = (data < lend) ;
% ones where data below left end
flagright = (data > rend) ;
% ones where data above right end
nleft = sum(flagleft) ;
nright = sum(flagright) ;
if nleft + nright > 0 ; % then need to deal with points outside range

if eptflag == 1 ; % the truncate data outside range
datatrunc = data(~(flagleft | flagright)) ;
% keep data that is not (outside left or outside right)

else ; % then move outside points to nearest end
datatrunc = data ;
if nleft > 0 ; % then replace those points with lend
datatrunc(flagleft) = lend * ones(nleft,1) ;
end ;
if nright > 0 ; % then replace those points with rend
datatrunc(flagright) = rend * ones(nright,1) ;
end ;

end ;

else ;

datatrunc = data ;

end ;

if ndo < n ; % then need to subsample
[temp,randperm] = sort(rand(n,1)) ;
% randperm is a random permutation of 1,2,...,n
vkeep = randperm(1:ndo) ;
% indices of elements to keep for display
dataol = datatrunc(vkeep,:) ;
else ; % overlay full data set
dataol = datatrunc ;
end ;

% overlay selected data
%
vax = axis ;
plotbottom = vax(3) ;
ftop = vax(4) ;
hold on ;
yrand = plotbottom + (0.7 + 0.2 * rand(ndo,1)) ...
* (plottop - plotbottom) ;
% y coords for jittering

if ibigdot == 1 ; % plot deliberately large dots
plot(dataol,yrand,[dolcolor 'o'],'MarkerSize',1,'LineWidth',2) ;
else ; % use matlab default dots
plot(dataol,yrand,[dolcolor '.']) ;
end ;

hold off ;

% reorder, so that overlaid data appears under the smooth
%
vachil = get(gca,'Children') ;
nchil = length(vachil) ;
vachil = [vachil(2:nchil); vachil(1)] ;
set(gca,'Children',vachil) ;

end ;

if ~isempty(titlestr) ;
if isempty(titlefontsize) ;
title(titlestr) ;
else ;
title(titlestr,'FontSize',titlefontsize) ;
end ;
end ;

if ~isempty(xlabelstr) ;
if isempty(labelfontsize) ;
xlabel(xlabelstr) ;
else ;
xlabel(xlabelstr,'FontSize',labelfontsize) ;
end ;
end ;

if ~isempty(ylabelstr) ;
if isempty(labelfontsize) ;
ylabel(ylabelstr) ;
else ;
ylabel(ylabelstr,'FontSize',labelfontsize) ;
end ;
end ;

end ;