function [f,x] = kernel(x);

% kernel(x) returns the length(x)-point Gaussian
% kernel density estimate where h is a smoothing parameter

n=length(x);
h=1.06*std(x)*n^(-1/5);
x=sort(x);
for i=1:n
   f(i)=1/(n*h*sqrt(2*pi))*sum(exp(-.5.*((x-x(i))./h).^2));
end;