close all
clear all
load thymusglandweights.dat

literature_mean = 34.3;

x = thymusglandweights;
n = length(x);
y = x - mean(x) + literature_mean;

m = 2000;
me_null = zeros(m,1);
me_h1 = zeros(m,1);
for i = 1:m
    me_null(i) = mean(y(randi(n,n,1)));
    me_h1(i) = mean(x(randi(n,n,1)));
end

bins = linspace(34,35,100);

null = hist(me_null, bins);
h1 = hist(me_h1, bins);
bar(bins, null, 'FaceColor',[.3,.3,.3]);
hold on
bar(bins, h1, 'FaceColor',[.7,.7,.7]);
mu = mean(x);
plot([mu,mu],[0,200],'--r','LineWidth',3);
xlabel('thymus gland weights [g]');
ylabel('frequency');
title('bootstrapped null distribution');
hold off

% 5% significance boundaries
low = quantile(me_null,0.025);
high =  quantile(me_null,0.975);
disp(['the 5% boundaries are: ', num2str(low), ' ', num2str(high)]);

hold on
plot([low,low],[0,200],'--g','LineWidth',3);
plot([high,high],[0,200],'--g','LineWidth',3);
hold off

idx = abs(me_h1-literature_mean) > abs(literature_mean - low);
pow = mean(idx);
h1positive = hist(me_h1(idx), bins);
hold on
bar(bins, h1positive, 'FaceColor','g');
hold off