resample = 500

load( 'thymusglandweights.dat' );
x = thymusglandweights;
nsamples = length( x );
sem = std(x)/sqrt(nsamples);

mu = zeros( resample, 1 );
for i = 1:resample
  % resample:
  xr = x(randi(nsamples, nsamples, 1));
  % compute statistics on sample:
  mu(i) = mean(xr);
end
bootsem = std( mu );

hold on
hist( x, 20 );
hist( mu, 20 );
hold off

disp(['bootstrap standard error: ', num2str(bootsem)]);
disp(['standard error: ', num2str(sem)]);