% kdentest.m: program to test kernel density estimator by comparing
%             maximum likelihood estimate of density of data from a normal
%             distribution 
%
%             John Rust, University of Maryland, November 2011

% first generate random sample from a N(mu,sigma) distribution

n=1000;  % number of observations to generate
mu=4;    % true mu parameter
sigma=3; % true sigma parameter

data=mu+sigma*randn(n,1);

% now compute the maximum likelihood estimates from the data

muhat=mean(data);
sighat=std(data);

% create a grid of points to plot estimated density

np=1000;  % number of points to plot
xmin=-5;  % lower bound of points to plot
xmax=15;  % upper bound of points to plot

x=(xmin:(xmax-xmin)/np:xmax)';

np1=size(x,1);
y=zeros(np1,1);
yt=zeros(np1,1);
fac=1/(sqrt(2*pi)*sighat);
fact=1/(sqrt(2*pi)*sigma);
for i=1:np1;
   y(i)=fac*exp(-(((x(i)-muhat)/sighat)^2)/2);
   yt(i)=fact*exp(-(((x(i)-mu)/sighat)^2)/2);
end;

plot(x,y,'-b','Linewidth',2);
hold on;
kde=kdensity(x,data);
plot(x,kde,'--r','Linewidth',2);
plot(x,yt,':k','Linewidth',1);
hold off;
legend('Maximum Likelihood','Nonparametric','True Density');
title('Maximum Likelihood vs Nonparametric Estimates of a Normal Density');

       locfac=.05;
       text(xmin+locfac*(xmax-xmin),max(kde)*.95,['Mean    ' num2str(mean(data))]);
       text(xmin+locfac*(xmax-xmin),max(kde)*.9,['Median  ' num2str(median(data))]);
       text(xmin+locfac*(xmax-xmin),max(kde)*.85,['Minimum ' num2str(min(data))]);
       text(xmin+locfac*(xmax-xmin),max(kde)*.8,['Maximum ' num2str(max(data))]);
       text(xmin+locfac*(xmax-xmin),max(kde)*.75,['Std dev ' num2str(std(data))]);
       text(xmin+locfac*(xmax-xmin),max(kde)*.7,['N       ' num2str(size(data,1))]);