automatic_replot = 0 ;
				# Demonstrate the central limit theorem
				# for any discrete distribution over
				# integers.

				# doit2 includes a Gaussian
				# approximation and 
				# doit3 shows everything on a log scale also.

				# To use this, run octave, then type doit3

				# Define x_n is a random variable from:
				# P(x) = p(0),p(1),p(2)....p(k)

p=[ 0.9 , 0.1 ]; 
p=[ 0.5, 0.5 ];                             use_y_for_min = 1 ;
p=[ 0.5 , 0, 0,0, 0.4 , 0.05 ,0,0,  0.05 ]; use_y_for_min = 0 ;
				# z= sum(x_n);
				# compute P(z) = p*p*p*p .... *p .
				# where "*" denotes convolution.

				# Housekeeping
verbose = 1;
gset nokey
clear y ;

Nmax = 2500;

				# Find prob distbn of r, the number of
				# heads,
				# given N tosses
y = p ;
				# Create the information required to
				# compute
				# the Gaussian approximation
s = size(y) ;
r = [0: s(2)-1 ] ;
mu = r * p' ;
rsquared = r.**2 ;
meanrsquared = rsquared * p' ;
variance = meanrsquared - mu**2 ;

for N=1:Nmax  # loop over N=1...Nmax
  s = size(y) ;

				# find the gaussian approximation
  r = [0: s(2)-1 ] ;
  muz = N*mu ; varz = N*variance ;
  g = exp( - (r - muz).**2/ (2.0* varz ) ) / sqrt( 2.0*pi *varz ) ;

  gset multiplot
  gset origin 0.1,0.1
  gset size 0.8,0.4
				# to tell gnuplot numbers that octave
				# knows,
				# we have to build a string, then
				# evaluate it.

  if( use_y_for_min )
    ourmin = min(y) ;
  else
    ourmin = min(g) ;
  endif
  if (ourmin < 1e-300 ) 
    ourmin = 1e-300 ;
  endif
  command = sprintf(" gset yrange [%g: %g]", ourmin , max(y) )  ;
  eval(command);

				# let's set the xrange in a smart way.
  command = sprintf(" gset xrange [0: %d]", s(2)-1  )  ;
  eval(command);

				# add a nice title
  command = sprintf("   gset title \"%d\" " , N )  ;
  eval(command);

  gset logscale y
				# do the plot!
  gplot y' w impulse 3, g' w l 1

  command = sprintf(" gset yrange [0: %g]", max(y) )  ;
  eval(command);

  gset title ''
  gset origin 0.1,0.54
  gset nologscale y
				# do the second plot!
  gplot y' w impulse 3, g' w l 1 



  gset nomultiplot
  if( verbose )
    ans= input("ready? (answer 1 to go faster)") ; 
    if (ans == 1) 
      verbose = 0 ;
    endif
  endif
  y= conv(y,p);  # do the convolution of y with p
endfor