% 
% 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 = 5000; 
X = unifrnd( 0, 2, N, 1 );

k_values = [ 32, 64, 256 ];

n_ks = length(k_values); 
colors = 'rbk';

figure;
pt = plot( x_truth, p_truth, '-g' ); hold on; 
ph = zeros(n_ks,1); % handles for each Parzen density estimation

for ki = 1:length(k_values)
  k = k_values(ki);     
  
  % For each possible h plot the K-NN density approximation on top of the truth:
  %    
  x_grid = linspace(0,2,100); % where we will sample our density \hat{p} 
  
  p_hat = [];
  for xi = 1:length(x_grid)
    x = x_grid(xi);
    
    % find the K-NN of x and the volume around them:
    %
    [dum,inds] = sort(abs(X - x));
    nn     = inds(1:k);
    xnn    = X(nn);
    V_of_x = max(xnn) - min(xnn);
    
    p = k/(N*V_of_x);
    
    p_hat = [ p_hat, p ];  
  end
  
  ph(ki) = plot( x_grid, p_hat, ['-',colors(ki)] ); 
  
end

legend( [pt,ph(1),ph(2),ph(3)], {'truth','k=32','k=64','k=256'} );
fn = ['chap_2_prob_32.eps'];
saveas(gcf,['../../WriteUp/Graphics/Chapter2/',fn],'epsc');