%
% epage 325
% 
% Written by:
% -- 
% John L. Weatherwax                2008-02-20
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

clear all; 
close all; 
clc; 

addpath('../../Code/CSTool');
load snowfall; 

x = 0:140;                             % <- the points where we want to evaluate our kernel density estimate
n = length(snowfall); 
h = 1.06*sqrt(var(snowfall))*n^(-1/5); % <- optimal kernel width for kernel density estimation ...
fhat = zeros(size(x)); 
for i=1:length(x)
  xloc = x(i)*ones(1,n);               % <- repeat this x to evaluate it at each density kernel.
  arg  = ((xloc-snowfall)/h).^2; 
  fhat(i) = sum(exp(-.5*arg))/(n*h*sqrt(2*pi));
end

figure; plot( x, fhat, '-rx' ); grid on; 
xlabel( 'snowfall value' ); ylabel( 'probability density' ); 
title( 'kernel density estimate of the snowfall data' ); 
saveas( gcf, 'prob_8_11_pdf', 'epsc' ); 

fprintf('the integral of our kernel density estimate of the snowfall data is %10.5f\n', trapz(x,fhat));