%
% 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); 

%--
% Using Example 12.13 to generate samples from an inhibition process: 
%--

delta = 0.1;
% Get the vertices for the region.
rx = [0 2 2 0 0];
ry = [0 0 2 2 0];
n = 100;

X = zeros(n,2);
% Generate the first event.
X(1,:) = csbinproc(rx,ry,1);
i = 1;
% Generate the other events.
while i<n
  [sx,sy] = csbinproc(rx, ry, 1);
  xt = [sx sy ; X(1:i,:)];
  % Find the distance between the events
  dist = pdist(xt);
  % Find the distance between the candidate event
  % and the others that have been generated already.
  ind = find(dist(1:i) <= delta);
  if isempty(ind)
    % Then we keep the event.
    i = i+1;
    X(i,:) = [sx, sy];
  end
end

%--
% Using Example 12.7 to compare \hat{G} with \hat{G}_CSR (event-event NN CDF)
%--

all_dists = pdist(X); 
nw = 100; 
dw = (max(all_dists)-min(all_dists))/nw;
w  = min(all_dists):dw:max(all_dists); nw = length(w); 
nx = size(X,1); 

ghatobs = csghat(X,w);

% Get the simulations.
B = 99;
% Each row is a Ghat from a simulated CSR process.
simul = zeros(B,nw);
for b = 1:B
  [xt,yt]    = csbinproc(rx, ry, nx);
  simul(b,:) = csghat([xt,yt],w);
end

% Get the average.
ghatmu = mean(simul);
% Get the envelopes.
ghatup = max(simul);
ghatlo = min(simul);
figure; 
plot(ghatmu,ghatobs,'k',ghatmu,ghatup,'g--',ghatmu,ghatlo,'g--'); hold on; 
plot(ghatmu,ghatmu,'r' ); 

xlabel('G CSR'); ylabel('G Obs (hat)'); grid on; 
saveas( gcf, '../../WriteUp/Graphics/Chapter12/prob_12_11_gs', 'epsc' ); 

clear functions;