done

statistics/matlab/ci_mean.m

@@ -1,6 +1,6 @@
 load thymusglandweights.dat
 
-n = 80;
+n = 16;
 x = thymusglandweights(1:n);
 
 m = 5000;

statistics/matlab/estimate_regression.m

@@ -0,0 +1,5 @@
+function param = estimate_regression(x,y, param0)
+    options = optimoptions(@fminunc,'GradObj','on');
+    myfunc = @(p)(lserr(p,x,y));
+    
+    param = fminunc(myfunc,param0, options);

statistics/matlab/lserr.m

@@ -0,0 +1,6 @@
+function [err, grad] = lserr(param, x, y)
+    err = mean( (param(1)*x + param(2) - y).^2 );
+    
+    if nargout == 2
+        grad = lserr_gradient(param, x,y);
+    end

statistics/matlab/lserr_gradient.m

@@ -0,0 +1,10 @@
+function grad = lserr_gradient(param, x, y)
+    h = 1e-6;
+    
+    grad = 0*param;
+    for i = 1:length(param)
+        paramh = param;
+        paramh(i) = param(i) + h;
+        grad(i) = (lserr(paramh,x,y) - lserr(param,x,y))/h;
+    end
+    

statistics/matlab/my_gradient_descent

@@ -0,0 +1,30 @@
+param0 = [1,1];
+step = 0.01;
+
+m = 50;
+arange = linspace(0,1,m);
+brange = linspace(.5,1.5, m);
+
+[A,B] = meshgrid(arange, brange);
+
+E = 0*A;
+
+x = linspace(-5,5,20);
+y = .5*x + 1 + randn(1, length(x));
+
+U = 0*A;
+V = 0*A;
+
+
+for i = 1:m
+    for j = 1:m
+        E(i,j) = lserr([A(i,j), B(i,j)], x, y);
+        grad = lserr_gradient([A(i,j), B(i,j)], x, y);
+        U(i,j) = grad(1);
+        V(i,j) = grad(2);   
+        
+    end
+end
+colormap('jet');
+
+surf(A,B,E, 'FaceAlpha',.5);

statistics/matlab/my_gradient_descent.m

@@ -0,0 +1,62 @@
+close all
+clear
+
+m = 50;
+arange = linspace(0,1,m);
+brange = linspace(.5,1.5, m);
+
+[A,B] = meshgrid(arange, brange);
+
+E = 0*A;
+
+x = linspace(-5,5,20);
+y = .5*x + 1 + randn(1, length(x));
+
+U = 0*A;
+V = 0*A;
+
+
+for i = 1:m
+    for j = 1:m
+        E(i,j) = lserr([A(i,j), B(i,j)], x, y);
+        grad = lserr_gradient([A(i,j), B(i,j)], x, y);
+        U(i,j) = grad(1);
+        V(i,j) = grad(2);   
+        
+    end
+end
+colormap('gray');
+
+subplot(1,2,1);
+hold on
+surf(A,B,E, 'FaceAlpha',.5);
+shading interp
+pause
+subplot(1,2,2);
+plot(x,y,'ok');
+
+%%
+
+t = linspace(-5,5,100);
+param0 = [0,0];
+step = 0.01;
+param = param0;
+
+
+for i = 1:100
+    err = lserr(param, x, y);
+    derr = lserr_gradient(param, x, y);
+    subplot(1,2,1);
+    plot3(param(1), param(2), err,'or');
+
+    subplot(1,2,2);
+    hold off
+    plot(x,y,'ok');
+    hold on
+    plot(t, param(1)*t + param(2), '--k', 'LineWidth',2);
+    pause(0.2);
+    param = param - step*derr;
+    
+end
+
+hold off

statistics/matlab/mypower.m

@@ -0,0 +1,49 @@
+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
+

statistics/matlab/plot_error_surface.m

@@ -0,0 +1,38 @@
+close all;
+clear;
+
+m = 50;
+arange = linspace(0,1,m);
+brange = linspace(.5,1.5, m);
+
+[A,B] = meshgrid(arange, brange);
+
+E = 0*A;
+
+x = linspace(-5,5,20);
+y = .5*x + 1 + randn(1, length(x));
+
+U = 0*A;
+V = 0*A;
+
+
+for i = 1:m
+    for j = 1:m
+        E(i,j) = lserr([A(i,j), B(i,j)], x, y);
+        grad = lserr_gradient([A(i,j), B(i,j)], x, y);
+        U(i,j) = grad(1);
+        V(i,j) = grad(2);   
+        
+    end
+end
+colormap('jet');
+
+surf(A,B,E, 'FaceAlpha',.5);
+%shading interp;
+hold on
+contour(A,B,E, 50, 'LineColor', 'k')
+quiver(A,B,U,V);
+xlabel('a'); 
+ylabel('b');
+zlabel('mean square error')
+axis([0,1,.5,1.5])

statistics/matlab/tests.m

@@ -11,7 +11,7 @@ y = x - mean(x) + literature_mean;
 m = 2000;
 me = zeros(m,1);
 for i = 1:m
-    me(i) = median(y(randi(n,n,1)));
+    me(i) = mean(y(randi(n,n,1)));
 end
 
 hist(me, 50);

statistics/matlab/winners_curse.m

@@ -0,0 +1,21 @@
+load thymusglandweights.dat
+
+n = 10;
+
+x = thymusglandweights;
+
+m = 5000;
+for i = 1:m
+    idx = randi(length(x), n,1);
+    y = x(idx);
+    [h,p] = ttest(y, 34.3);
+    
+    if h == 1
+        disp(['p-value: ', num2str(p)]);
+        disp(['mu: ', num2str(mean(y))]);
+        disp(['mu total: ', num2str(mean(x))]);
+        break
+    end
+end
+
+