%clf; 
 clear all;
 close all; 
 
 
 disp('NP Bootstrap for Hubble Correlation')
  lw = 2; 
  set(0, 'DefaultAxesFontSize', 17);
  fs = 17;
  msize = 7;
  
mpc = [.032  .034   .214  .263   .275 .275 ...
 .45   .5     .5    .63  .8   .9   ...
 .9    .9    .9     1.0    1.1    1.1 ...
 1.4   1.7   2.0    2.0   2.0  2.0 ];

v = [170  290  -130  -70 -185  -220 ...
  200       290     270    200    300     -30 ...
  650      150      500    920    450     500 ...
  500    960    500   850     800     1090 ];

pairs = [mpc' v'];

figure(1)
plot(mpc, v, '*','markersize', msize) 
xlabel('megaparsec');
ylabel('velocity');
print -depsc 'C:\NPBook\Bootstrap\hubblef1.eps'
cc = corrcoef(pairs);
ccc=cc(1,2)
% bootstrap sample and bootstrap estimators
bsam=[];
rand('state', 0 ) % setting randon number generator
B=50000;
for b = 1:B
   bs =  bootsample(pairs);
   ccbs = corrcoef(bs);
   bsam = [bsam ccbs(1,2)];
end
figure(2)
hist(bsam, 80)
%title('Histogram of correlation coefs of 10000 bootstrap samples')
bsd = std(bsam)
print -depsc 'C:\NPBook\Bootstrap\hubblef2.eps'

figure(3)
% Fisher transformation: hat zeta= 0.5*log((1+hat theta)/(1-hat theta))
% is approximatelly normally distributed with mean 0.5*log((1+theta)/(1-theta))
% and standard deviation 1/sqrt(n-3)
%
zeta = 0.5 * log( (1+bsam)./(1-bsam) );
hist(zeta,80)
print -depsc 'C:\NPBook\Bootstrap\hubblef3.eps'

% st.deviation of zeta
std(zeta)
%should be close to 1/sqrt(n-3)
1/sqrt(length(mpc)-3)