%
% Example on epage XXX
% 
% Problem on epage 517
% 
% 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');

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

load bodmin
X = [x,y];
n = length(x); 

% Get the convex hull for the boundary of our "region".
% -- K is the indices to points ON the convex hull
K    = convhull(x, y);
cvh  = [x(K), y(K)]; 

% this is a very poor way (computationally) of getting estimates of the 
% min/max of the nearest neighbor statistics 
dist = pdist(X);
% Convert to a matrix and put large numbers on the diagonal.
D = diag(realmax*ones(1,n)) + squareform(dist);
% Find the smallest distances in each row or col.
mind = min(D);

% we will use these estimates as the estimate of the min/max of the 
% point-event statistics ... again something smarter could be done ... 
% 
dw = (max(mind)-min(mind))/30; 
w  = min(mind):dw:max(mind); nw = length(w); 

% get "2*n" random points and compute their nearest event distance
% this is a surrogate for the "observed" f_hat: 
% we also get x_rnd ... the random points specified ... 
[fhatobs,x_rnd] = csfhat(X,w,n,cvh);

% Now perform a similar procedure on B bootstrap samples of random data: 
%
B = 100;
% Each row is a Fhat from a simulated CSR process.
simul = zeros(B,nw);
for b = 1:B
  % use the auxillary output from this routine to generate "n" random "event" points: 
  [dum,x_events] = csfhat(X,w,n,cvh); clear dum; 
  % the above random points go into the event location in the next call to "csfhat" 
  simul(b,:)     = csfhat(x_events,w,n,cvh);
end

% Get the average of the bootstrap samples:
fhatmu = mean(simul);
% Get the envelopes.
fhatup = max(simul);
fhatlo = min(simul);

% plot the summary statistics: 
% 
plot(fhatmu,fhatmu,'-k',fhatmu,fhatobs,'-rx',fhatmu,fhatup,'-go',fhatmu,fhatlo,'-go');
xlabel('Fhat under CSR'); ylabel('Fhat Observed'); 
axis tight; grid on; 
saveas( gcf, '../../WriteUp/Graphics/Chapter12/prob_12_7_fhat', 'epsc' ); 

clear functions;