%two Gaussians, their histogram, and kernel density estimate
%example illustrating concepts of Ch. 15 in 
%Dekking et al. A Modern Introduction to Probability and Statistics. Second Edition. Springer 2007 
%Longin Jan Latecki, latecki@temple.edu, March 2012
clear all
close all

h=0.5; %bandwidth for the kernel density estimate f
binwidth=0.5; %width of histogram bins of histogram hx
sigma1=1;
sigma2=0.5; 
%g is a mixture of two Gaussians 0.5N(3,1)+0.5N(7,0.5)
g=@(v) 0.5*normpdf(v,3,sigma1)+0.5*normpdf(v,7,sigma2); 
v=[1:100]/10; %vector of points used to plot functions
gv=g(v); %compute function values at points v

%vector of 200 sample points from g
x = [3 + sigma1 * randn(1,100), 7 + sigma2 * randn(1,100)];

%normalized histogram so that the area of its bars is 1, see top of p. 211
bincenters=[0.5:binwidth:9.5];
hx=hist(x,bincenters)/(200*binwidth); 

%f is a kernel density estimate at point t, see formula (15.1) on p. 215
f = @(t) 1/(200*h)*sum(normpdf((t-x)/h,0,1)); 
%compute f values at points v
for i=1:length(v)
    fv(i)=f(v(i));
end

figure; 
bar(bincenters,hx,'y')
hold on
plot(v,gv,'LineWidth',3)
plot(v,fv,'r','LineWidth',2)
title('two guassians g=0.5N(3,1)+0.5N(7,0.5)');
leg1 = legend('histogram','g','f');
set(leg1,'Location','NorthWest')
plot(x,zeros(1,200)+0.01*rand(1,200),'.')




