rng('default');
rng(1);

x = [2877+213.3*random('Normal',0,1,2000000,1); ...
     3871+491.9*random('Normal',0,1,1940764,1)];

x = round(x);

figure(1);
histogram(x,'EdgeColor','auto');
xlim([0 8000]);
text(5000,30000, strcat('$\mu$   = ',sprintf('%8.2f \n', mean(x))),'FontSize',18,'Interpreter','latex');
text(5000,27000, strcat('$\sigma$  =',sprintf('%8.2f \n', std(x,1))),'FontSize',18,'Interpreter','latex');
title('Figure 7.1: Population (4 million)');
xlabel('Birthweight (g)');
set(gca,'FontSize',18);
drawnow;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

n = length(x);
indices = random('Discrete Uniform',n,100,1);
x100 = x(indices);
mp = max(180,0.7*max(x100));

figure(2);
histogram(x100,'BinWidth',200);
xlim([0 8000]);
text(4500,16, strcat('$\bar{X}$   = ',sprintf('%8.2f \n', mean(x100))),'FontSize',18,'Interpreter','latex');
text(4500,14, strcat('$S$  =',sprintf('%8.2f \n', std(x100))),'FontSize',18,'Interpreter','latex');
title('Sample of Size 100');
xlabel('Birthweight (g)');
set(gca,'FontSize',18);
drawnow;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

n = length(x);
for j=1:10000
    indices = random('Discrete Uniform',n,100,1);
    x100new = x(indices);
    x100mean(j) = mean(x100new);
end

figure(3);
histogram(x100mean,'EdgeColor','auto');
xlim([0 8000]);
text(4000,400, sprintf('mean   = %8.2f \n', mean(x100mean)),'FontSize',18);
text(4000,350, sprintf('std     = %8.2f \n', std(x100mean)),'FontSize',18);
title('Figure 7.2: Sample Mean');
xlabel('Birthweight (g)');
set(gca,'FontSize',18);
drawnow;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The is bootstrap

for j=1:10000
    indices = random('Discrete Uniform',100,100,1);
    x100resampled = x100(indices);
    x100mean(j) = mean(x100resampled);
end

figure(4);
histogram(x100mean,'EdgeColor','auto');
xlim([0 8000]);
text(4000,500, sprintf('mean   = %8.2f \n', mean(x100mean)),'FontSize',18);
text(4000,400, sprintf('std     = %8.2f \n', std(x100mean)),'FontSize',18);
title('Figure 7.3: Means of Bootstrapped Samples');
xlabel('Birthweight (g)');
set(gca,'FontSize',18);
drawnow;

figure(5);
histogram(x100mean);
xlim([3000 3800]);
text(3500,500, sprintf('mean   = %8.2f \n', mean(x100mean)),'FontSize',18);
text(3500,400, sprintf('std     = %8.2f \n', std(x100mean)),'FontSize',18);
text(3500,300, sprintf(' 2.5 pct = %0.2f \n', prctile(x100mean,2.5)),'FontSize',18);
text(3500,200, sprintf('97.5 pct = %0.2f \n', prctile(x100mean,97.5)),'FontSize',18);
title('Figure 7.3: Means of Bootstrapped Samples');
xlabel('Birthweight (g)');
set(gca,'FontSize',18);
drawnow;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The is SERIOUS bootstrap

%for j=1:300000
for j=1:30000
    indices = random('Discrete Uniform',100,100,1);
    x100resampled = x100(indices);
    x100mean2(j) = mean(x100resampled);
end

figure(6);
histogram(x100mean2);
xlim([2900 3700]);
text(3400,5000, sprintf('mean   = %8.2f \n', mean(x100mean2)),'FontSize',18);
text(3400,4500, sprintf('std     = %8.2f \n', std(x100mean2)),'FontSize',18);
text(3400,4000, sprintf(' 2.5 pct = %0.2f \n', prctile(x100mean2,2.5)),'FontSize',18);
text(3400,3500, sprintf('97.5 pct = %0.2f \n', prctile(x100mean2,97.5)),'FontSize',18);
title('Figure 7.3: Means of Bootstrapped Samples');
xlabel('Birthweight (g)');
set(gca,'FontSize',18);
drawnow;

for j=1:10000
    indices = random('Discrete Uniform',100,100,1);
    x100resampled = x100(indices);
    x100mean(j) = mean(x100resampled);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% How about bootstrapping an estimate of sigma^2

for j=1:10000
    indices = random('Discrete Uniform',100,100,1);
    x100resampled = x100(indices);
    x100std(j) = std(x100resampled);
end

figure(7);
histogram(x100std,'EdgeColor','auto');
xlim([400 900]);
text(700,500, sprintf('mean   = %8.2f \n', mean(x100std)),'FontSize',18);
text(700,400, sprintf('std     = %8.2f \n', std(x100std)),'FontSize',18);
title('Figure 7.3: Std. Devs. of Bootstrapped Samples');
xlabel('Std. Dev. of Birthweight (g)');
set(gca,'FontSize',18);
drawnow;