First lecture on descriptive statistics

statistics/code/checkmymedian.m

@@ -0,0 +1,12 @@
+% check whether the median returned by mymedian 
+% really separates a vector into two halfs
+for i = 1:140                                    % loop over different length
+  for k = 1:10                                   % try several times
+    a = randn( i, 1 );                           % generate some data
+    m = mymedian( a )                            % compute median
+    if length( a(a>m) ) ~= length( a(a<m) )      % check
+      disp( 'error!' )
+    end
+  end
+end
+

statistics/code/diehistograms.m

@@ -0,0 +1,24 @@
+% dependence of histogram on number of rolls:
+nrolls = [ 20, 100, 1000 ];
+for i = [1:length(nrolls)]
+  d = rollthedie( nrolls(i) );
+  % plain hist:
+  % hist( d )
+
+  % check bin counts of plain hist:
+  % h = hist( d )
+
+  % force 6 bins:
+  % hist( d, 6 )
+
+  % set the right bin centers:
+  bins = 1:6;
+  %hist( d, bins )
+
+  % normalize histogram and compare to expectation:
+  hold on
+  plot( [0 7], [1/6 1/6], '-r', 'linewidth', 10 )
+  hist( d, bins, 1.0, 'facecolor', 'b' )
+  hold off
+  pause
+end

statistics/code/gaussianbins.m

@@ -0,0 +1,17 @@
+x = randn( 100, 1 );
+bins1 = -4:2:4;
+bins2 = -4:0.5:4;
+subplot( 1, 2, 1 );
+hold on;
+hist( x, bins1 );
+hist( x, bins2 );
+xlabel('x')
+ylabel('Frequeny')
+hold off;
+subplot( 1, 2, 2 );
+hold on;
+hist( x, bins1, 1.0/(bins1(2)-bins1(1)) );
+hist( x, bins2, 1.0/(bins2(2)-bins2(1)) );
+xlabel('x')
+ylabel('Probability density')
+hold off;

statistics/code/gaussianpdf.m

@@ -0,0 +1,22 @@
+% plot Gaussian pdf:
+dx=0.1
+x = [-4.0:dx:4.0];
+p = exp(-0.5*x.^2)/sqrt(2.0*pi);
+hold on
+plot(x,p, 'linewidth', 10 )
+
+% compute integral between x1 and x2:
+x1=1.0
+x2=2.0
+P = sum(p((x>=x1)&(x<x2)))*dx
+
+% draw random numbers:
+r = randn( 10000, 1 );
+hist(r,x,1.0/dx)
+
+% check P:
+Pr = sum((r>=x1)&(r<x2))/length(r)
+
+hold off
+
+

statistics/code/histogramquartiles.m

@@ -0,0 +1,24 @@
+% generate data:
+x = randn( 1, 100000 );
+
+% histogram:
+[h,b] = hist( x, 100 );
+% normalize:
+bs = b(2)-b(1);
+h = h/sum(h)/bs;
+
+% plot:
+bar( b, h );
+xlabel( 'x' );
+
+% median, quartile:
+q = quartiles( x );
+%q = quantile( x, [0.25, 0.5, 0.75 ] );
+
+% plot:
+hold on;
+bar( b(b<q(1)), h(b<q(1)), 'FaceColor', [0.5 0 0.5] );
+bar( b((b>=q(1)) & (b<q(2))), h((b>=q(1)) & (b<q(2))), 'FaceColor', [0.9 0 0] );
+bar( b((b>=q(2)) & (b<q(3))), h((b>=q(2)) & (b<q(3))), 'FaceColor', [0 0 0.9] );
+bar( b(b>=q(3)), h(b>=q(3)), 'FaceColor', [0.5 0 0.5] );
+hold off;

statistics/code/mymedian.m

@@ -0,0 +1,13 @@
+function m = mymedian( x )
+% returns the median of the vector x
+  xs = sort( x );
+  if ( length( xs ) == 0 )
+    m = NaN;
+  elseif ( rem( length( xs ), 2 ) == 0 )
+    index = length( xs )/2;
+    m = (xs( index ) + xs( index+1 ))/2;
+  else
+    index = (length( xs ) + 1)/2;
+    m = xs( index );
+  end
+end

statistics/code/quartiles.m

@@ -1,25 +1,15 @@
-% generate data:
-x = randn( 1, 100000 );
-
-% histogram:
-[h,b] = hist( x, 100 );
-% normalize:
-bs = b(2)-b(1);
-h = h/sum(h)/bs;
-
-% plot:
-bar( b, h );
-xlabel( 'x' );
-
-% median, quartile:
-xs = sort( x )
-q = [ xs(length(xs)/4), xs(length(xs)/2), xs(3*length(xs)/4) ];
-%q = quantile( x, [0.25, 0.5, 0.75 ] );
-
-% plot:
-bar( b(b<q(1)), h(b<q(1)), 'FaceColor', [0.5 0 0.5] );
-hold on;
-bar( b((b>=q(1)) & (b<q(2))), h((b>=q(1)) & (b<q(2))), 'FaceColor', [0.9 0 0] );
-bar( b((b>=q(2)) & (b<q(3))), h((b>=q(2)) & (b<q(3))), 'FaceColor', [0 0 0.9] );
-bar( b(b>=q(3)), h(b>=q(3)), 'FaceColor', [0.5 0 0.5] );
-hold off;
+function q = quartiles( x )
+  % returns a vector with the first, second, and third quartile of the vector x
+  xs = sort( x );
+  if ( length( xs ) == 0 )
+    q = [];
+  elseif ( rem( length( xs ), 2 ) == 0 )
+    index = length( xs )/2;
+    m = (xs( index ) + xs( index+1 ))/2;
+    q = [ round( xs(length(xs)/4) ), m, xs(round(3*length(xs)/4)) ];
+  else
+    index = (length( xs ) + 1)/2;
+    m = xs( index );
+    q = [ round( xs(length(xs)/4) ), m, xs(round(3*length(xs)/4)) ];
+  end
+end

statistics/code/randomwalk.m

@@ -1,4 +1,6 @@
 function x = randomwalk(n,p)
+% returns a random wolk with n steps and 
+% probability p for positive steps.
     r = rand(n,1);
     r(r<p) = -1.0;
     r(r>=p) = +1.0;

statistics/code/rollthedie.m

@@ -0,0 +1,4 @@
+function x = rollthedie( n )
+% return a vector with the result of rolling a die n times
+  x = randi( [1, 6], n, 1 );
+end