%
% Example on epage XXX
% 
% Problem on epage 518
% 
% 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; drawnow;
clc; 

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

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

%--
% First generate data from a Poisson CLUSTER process:
%--

npar = 15;
% Get the vertices for the regions.
rx = [0 1 1 0 0];
ry = [0 0 1 1 0];
rxp = [-.05 1.05 1.05 -.05 -.05];
ryp = [-.05 -.05 1.05 1.05 -.05];
% Get all of the parents.
[xp,yp] = csbinproc(rxp, ryp, npar);

lam = 15;
% Get the number of children per parent.
nchild = poissrnd(lam,1,npar);
X = [];
r = 0.05;
sig = r*eye(2);
% Locate the children.
for i = 1:npar
  xc = randn(nchild(i),2)*sig + repmat([xp(i) yp(i)],nchild(i),1);
  X = [X; xc];
end
% Find the ones that are in the region of interest.
ind = find(inpolygon(X(:,1), X(:,2), rx, ry));
% Those are the children for the sample.
x = X(ind,1);
y = X(ind,2);



%--
% Second use the statistic found in Problem~9 to test for CSR 
% we expect this test to find {\em large} values of the statistic "T" 
%--

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 "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);
% Each element in T_b is the T statistic suggested in the book 
T_b   = zeros(B,1); 
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);
  T_b(b)         = max( abs(simul(b,:) - fhatobs) );
end

% plot the density of the bootstrap samples:
% 
n_sturges = round(1 + log2(B));
h_nrr     = 3.5*std(T_b)*B^(-1/3); n_nrr = round( (max(T_b)-min(T_b))/h_nrr ); 
fh = figure; 
cshistden( T_b, n_sturges ); hold on; 
[fhat,bin] = cshistden( T_b, n_sturges ); 
plot(bin,fhat,'-k', 'linewidth', 6 ); grid on; 
xlabel('T values'); ylabel('Estimated PDF'); 
saveas( gcf, '../../WriteUp/Graphics/Chapter12/prob_12_10_T_den', 'epsc' ); 

% get the average of all of the bootstrap samples: 
T_mean = mean(T_b); 
T_std  = std(T_b); 

fprintf('the mean value of T is %10.6f\n',T_mean);
fprintf('    while the 3 sigma lower bound is %10.6f\n',T_mean-3*T_std); 

clear functions;