% 
% Written by:
% -- 
% John L. Weatherwax                2007-07-01
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
% Problem: Epage 71
%
%-----

close all; clc; clear; 

rand('seed',0);
randn('seed',0);

x_truth = linspace( -1, 3, 1000 );
indsGT  = x_truth>0.0;
indsLT  = x_truth<2.0;
inds    = find( indsGT & indsLT );
p_truth = zeros( 1, 1000 ); 
p_truth(inds) = 1/2;

% Generate some data from the true distribution:
%


N_values = [ 32, 256, 5000 ];
h_values = [ 0.05, 0.2 ]; 

n_hs = length(h_values); 
colors = 'rb';
for ni = 1:length(N_values)
  N = N_values(ni); 
  X = unifrnd( 0, 2, N, 1 );

  figure;
  pt = plot( x_truth, p_truth, '-g' ); hold on; 
  ph = zeros(n_hs,1); % handles for each Parzen density estimation
  
  for hi = 1:length(h_values)
    h = h_values(hi);     

    % For each possible h plot the Parzen approximation on top of the truth:
    %    
    x_grid = linspace(-1,3,100); % where we will sample our density \hat{p} 

    p_hat = [];
    for xi = 1:length(x_grid)
      x = x_grid(xi);
      p = mean( normpdf( ( repmat( x, N, 1 ) - X ) / h ) )/h; 

      p_hat = [ p_hat, p ];  
    end

    ph(hi) = plot( x_grid, p_hat, ['-',colors(hi)] ); 
    
  end
  
  legend( [pt,ph(1),ph(2)], {'truth',['N=',num2str(N),'; h=0.05'],['N=',num2str(N), '; h=0.2']} );
  fn = ['chap_2_prob_31_N_',num2str(N),'.eps'];
  saveas(gcf,['../../WriteUp/Graphics/Chapter2/',fn],'epsc');
  
end