%
% problem 268
% example epage 250 
% 
% 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');

N_MC = 1000; 

bts_stderr = zeros(1,N_MC); bts_bias = zeros(1,N_MC); 
jkn_stderr = zeros(1,N_MC); jkn_bias = zeros(1,N_MC); 
for mci=1:N_MC,
  n = 100; 
  x = rand(1,n); 

  T = mom(x); 
  
  % compute the bootstrapped standard error (of the momiance): 
  B = 400; 
  [m] = bootstrp( B, @mom, x ); 
  bts_stderr(mci) = std(m); 
  bts_bias(mci) = mean(m)-T;

  % compute the jackknifed estimate standard eor (of the momiance): 
  jack = zeros(1,n); % <- store the jackknifed samples 
  for ii=1:n
    jack(ii) = mom( x( setdiff(1:n,ii) ) ); 
  end
  % get the mean of the jackknifed estimates
  TJ = mean(jack);  
  jackBias = (n-1)*(TJ-T);
  jackEstStdError = sqrt( ((n-1)/n)*sum( ( jack - TJ ).^2 ) ); 
  jkn_stderr(mci) = jackEstStdError;
  jkn_bias(mci) = jackBias; 
end

figure; 
[N,h] = hist( bts_stderr ); 
bar(h,N,.25,'g');
hold on; 
[N,h] = hist( jkn_stderr ); 
bar(h,N,.25,'r'); 
title( 'histogram of bootstrap (green) and jackknife (red) standard error estimates' )
saveas( gcf, 'mom_bt_jk_std_err_hist', 'eps' ); 

figure; 
[N,h] = hist( bts_bias ); 
bar(h,N,.25,'g');
hold on; 
[N,h] = hist( jkn_bias ); 
bar(h,N,.25,'r'); 
title( 'histogram of bootstrap (green) and jackknife (red) bias estimates' )
saveas( gcf, 'mom_bt_jk_bias_hist', 'eps' );