# 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 convolve
## November 2004: this no longer 
## works right: multiplot has been broken. So use convolve2 instead.

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
]; 

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

T = 400 ; # maximum number of variables from p0 to add together
y = p0' ;
gset size 1,1
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 ) ;

  command = sprintf ( "gset title 'sum of %d' " , t ) ;
  eval(command) ;
  command = sprintf ( "gset xrange [ 0:%d ] " , (s(2)-1) ) ;
  eval(command) ;
  command = sprintf ( "gset yrange [ %g:%g] " , min(g)/10.0 ,  max(y) ) ;
  eval(command) ;
  gset multiplot
  gset origin 0,0.5
  gset size 1,0.5
  gset nologscale y ;
  gplot g' , y' w impulse ;

  gset logscale y ;

  gset origin 0,0
  gset size 1,0.5
  gplot g' , y' w impulse ;
  gset nomultiplot

  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<21))
    ans = input ( "ready" ) ;
  endif
endfor