function fh = kernelize(grid,sample,varargin);
% kernel density estimator.
%
% kernelize(grid,sample,h,method) returns a vector of the density estimate at
% each point of the "grid" vector, calculated from the vector of
% observations "sample", using the kernel specified by "method".
%
%
% "grid" is a column vector of points at which the density is to be
% evaluated.
%
% "sample" is a column vector of i.i.d. observations from the distribution.
%
% "h" is the bandwidth. The default is the IMSE-minimizing bandwidth for
% the case of a normally distributed sample. Cf. Silverman (1986), section
% 3.4.2.
%
% method is one of the following:
% 	---'normal' or 'n': normal kernel (default)
% 	---'epanechnikov' or 'e': Epanechnikov  kernel
%	---'nk': normal kernel, for distributions with compact support.
%          Applies a reflection around highest and lowest elements of the GRID.

% copyright 2001 Leonardo Rezende.
% This program can be freely copied, used and modified, as long as this
% copyright notice is kept in place. Please inform the author of any
% useful modifications, at lrezende@stanford.edu. 

if isempty(varargin{2})
	met='n';
else 
	met=varargin{2};
end;

if isempty(varargin{1})
	switch lower(met)
	case {'normal','n','nk'}
        	d=4/3;
	case {'epanechnikov','e'}
		d=40*sqrt(pi);
	end;
	h = (d/size(sample,1))^.2*std(sample);
else
	h = varargin{1};
end;

if met == 'nk'
	lb = find(sample>max(grid)-2*h);
	b = 2*max(grid)-sample(lb);
	la = find(sample<min(grid)+2*h);
	a = 2*min(grid)-sample(la);
	adj = (size(sample,1)+size(a,1)+size(b,1))/size(sample,1);
	sample = [a;sample;b];
else
	adj = 1;
end;

X=(ones(size(sample))*grid'-sample*ones(size(grid')))/h;

switch lower(met)
case {'normal','n','nk'}
	k = exp(-(X).^2/2)/sqrt(2*pi); 
case {'epanechnikov','e'}
	S = X.*X;
	k = 3/4*(1-S).*(S<=1);
end

fh=adj*mean(k/h)';




