# Central limit theorem demonstration.
				# 
				# Pick any discrete distribution and 
				# visualize the distribution of a sum of
				# t iid variables.
				# For t greater than 20 or so, the
				# distbn looks very Gaussian indeed.
				# However a close look (for example  on
				# log scale) shows persisting and
				# interesting differences. These graphs
				# would also be good for a discussion of
				# large deviation theory.
#
				# David MacKay    22 November 2003

				# To use, run octave, then type convolve2

gset autos xy;
verbose=0;

### one lop sided probability distribution
p0 = [  0.30000
  0.10000
  0.15000
  0.00000
  0.00000
  0.20000
  0.00000
  0.00000
  0.00000
  0.00000
  0.25000
];

### another lop-sided probability distribution, lop-sided the other way
p0 = [
  0.25000
  0.00000
  0.00000
  0.00000
  0.00000
  0.10000
  0.15000
  0.00000
  0.00000
  0.20000
  0.30000
]; 

#p0=[0.5 ; 0.5];
#p0=[0.9 ; 0.1];

gset lmargin 12
gset bmargin 2
gset rmargin 2
gset tmargin 1
gset nokey
gset tics out

T = 400 ; # maximum number of variables from p0 to add together
y = p0' ;
  gset origin 0,0
  gset size 0.8,0.8
for t=1:T
  s=size(y);
##  r = 1:(s(2));
  
  r = 0:(s(2)-1); ## the list of possible values of r, where 
                  ## r = sum_{n=1}^N x_n
  
  mean = r * y' ;
  rsquared = r.*r ;
  meansquare = rsquared * y' ;
  myvar = meansquare - mean *mean ;
  g = exp( - (r-mean).**2 / (2.0 * myvar) ) / sqrt(2.0*3.14159 * myvar ) ;

  figure(1);
  gset nologscale y ;
  command = sprintf ( "gset title 'sum of %d' " , t ) ;
  eval(command) ;
  command = sprintf ( "gset xrange [ 0-0.5:%d+0.5 ] " , (s(2)-1) ) ;
  eval(command) ;
  gmin = min(g)/10;
  if (gmin < 1e-200)
    gmin = 1e-200;
  endif
  command = sprintf ( "gset yrange [ %g:%g] " , gmin ,  max(y)*1.1 ) ;
  eval(command) ;
  gplot g' w l 4 , y' w impulse 3 ;

  figure(2);
  gset logscale y ;
  gplot g' w l 4 , y' w impulse 3 ;

  y = conv( y , p0 ) ;
  if((t<2))
    disp ( "This is the distribution of one variable, and its \
 approximation by a Gaussian. (Upper figure: same, on log scale)") ;
    disp ( "The remaining pictures show the distbn of the sum of t such variables") ;
  endif
  if(verbose || (t<210))
    ans = input ( "ready" ) ;
  endif
endfor