day 2 done

statistics/matlab/bootstrap_mean.m

@@ -0,0 +1,13 @@
+load thymusglandweights.dat
+x = thymusglandweights(1:50);
+
+m = 500;
+n = length(x);
+
+mu = zeros(m,1);
+for i = 1:m
+    mu(i) = mean(x(randi(n,n,1)));
+end
+fprintf("bootstrap standard error: %.4f\n", std(mu));
+fprintf("standard error: %.4f\n", std(x)/sqrt(n));
+

statistics/matlab/ci_mean.m

@@ -0,0 +1,19 @@
+load thymusglandweights.dat
+
+n = 80;
+x = thymusglandweights(1:n);
+
+m = 5000;
+me = zeros(m,1);
+for i = 1:m
+    me(i) = median(x(randi(n,n,1)));
+end
+
+t025 = tinv(0.025, n-1);
+t975 = tinv(0.975, n-1);
+
+se = std(x)/sqrt(n);
+
+fprintf('bootstrap quantiles: %.4f, %.4f \n', quantile(me,0.025), quantile(me,0.975));
+fprintf('analytical quantile: %.4f, %.4f \n', mean(x)+t025*se, mean(x)+t975*se);
+

statistics/matlab/ci_media.m

@@ -0,0 +1,17 @@
+load thymusglandweights.dat
+x = thymusglandweights(1:50);
+
+m = 500;
+n = length(x);
+x = sort(x);
+me = zeros(m,1);
+for i = 1:m
+    me(i) = median(x(randi(n,n,1)));
+end
+
+a1 = binoinv(0.025,n,.5)-1;
+a2 = binoinv(1-0.025,n,.5);
+
+fprintf('bootstrap quantiles: %.4f, %.4f \n', quantile(me,0.025), quantile(me,1-0.025));
+fprintf('analytical quantile: %.4f, %.4f \n', x(a1),x(a2));
+

statistics/matlab/tests.m

@@ -0,0 +1,38 @@
+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 = zeros(m,1);
+for i = 1:m
+    me(i) = median(y(randi(n,n,1)));
+end
+
+hist(me, 50);
+hold on
+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,0.025);
+high =  quantile(me,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
+
+pval = mean(abs(me-literature_mean) > abs(mu - literature_mean))
+
+legend('Null distribution','measured mean','5% significance boundaries')