
disp('Bootstrap of variance values');
disp('Sample size:');
n=7;

% Initialize random-number generator
rand('state',0);

% Generate a row vector of n independent N(0,100)s
disp('Basic sample:');
% rand(m,n) generates a mxn matrix of U(0,1)s
% randn(m,n) generates a mxn matrix of N(0,1)s
xx=10*randn(1,n)

n   % Display n

% Compute the mean and sample variance of the row vector
xm=sum(xx)/n;  xc=xx-xm;
xvar=sum(xc.*xc)/(n-1);
disp('Mean, sample standard deviation, and sample variance of sample:');
disp(' ');
disp([xm  sqrt(xvar)  xvar]);

disp('Basic sample sorted:');
xxsorted=sort(xx)

% A bootstrap sample of xx is a sample of size n
%   randomly sampled WITH REPLACEMENT from xx
% ceil(A) is the next highest integer above A, so that
%   ceil(n.*rand(1,n)) is n independent random integers chosen from 1:n
%   and xx(ceil(n.*rand(1,n))) is a bootstrap resample
disp('Five bootstrap resamples:')
for i=1:5;  bb(i,:)=xx(ceil(n.*rand(1,n)));  end;
disp(bb)

% Compute the mean and sample variance of each row of bb
for i=1:5;  bmean=sum(bb(i,:))/n;
   dd=bb(i,:)-bmean;  bvar(i)=sum(dd.*dd)/(n-1);  end;
disp('Sample variances of five bootstrap samples:');
disp('(The second column is the first column sorted.)');
disp([bvar'  sort(bvar)']);

% Find the bootstrap percentile 95% confidence interval
%   for the true unknown variance
% First, generate nboot (e.g. 1000 or 10000) bootstrap resamples
%   and store the sample variances of the resamples in an array
disp(' ');
disp('Finding bootstrap CI for true variance:');
disp('Number of bootstrap replications (resamples):');

nboot=10000

fprintf('Finding sample variances of %d bootstrap replications:\n',nboot);
% Each bootstrap resample is generated, used, and overwritten
%   so that all the memory we need is for a one-dimensional array
%   of `nboot'' sample variances
for b=1:nboot;
   zb=xx(ceil(n.*rand(1,n)));
   zmean=sum(zb)/n;  zc=zb-zmean;
   zboot(b) = sum(zc.*zc)/(n-1);  end;

disp('Mean and Median of bootstrapped sample variances:');
disp(' ');
zbmean = sum(zboot)/nboot;
zbmed = median(zboot);
disp([zbmean zbmed]);

disp('Bootstrap median bias (proportion <= Observed value):')
mb=0;
for b=1:nboot;  if zboot(b)<=xvar mb=mb+1;  end;  end;
disp(' ');
fprintf('  %d/%d = %6.4f  (Observed Variance: %g)\n', mb,nboot, mb/nboot, xvar);

% Convert to sorted array
zb2=sort(zboot);

Smallest_7_of_Nboot = zb2(1:7)
Largest_7_of_Nboot = zb2(nboot-7+1:nboot)

disp('Bootstrap 95% (percentile) CI bracketing Median:');
disp(' ');
disp([zb2(floor(nboot*0.025))  zbmed  zb2(ceil(nboot*0.975))]);

disp('This suggests that there might be a downwards bias in the CI.');
disp('HOWEVER zbmed is close to the original true variance.');



