November 18, 2004---Class 24---Exponential Distribution, Box-Muller Method,

Activities:

	Exponential Distribution

	The inverse transform method works well for the exponential
	distribution.  We used Mathematica to generate such numbers.

	v = Table[ -Log[ Random[] ], {1000}] ;

	Export["exponential.dat", v, "List"]

	This will put the vector v in a file called exponential.dat.
	Then use hist to plot the distribution.  (If you want to use a
	logarithmic y axis, make sure to use -g 2 option for hist and don't
	let axis try to plot down to 0.)

	Another way to get output from Mathematica is with the Write function.

	Do[ Write[ "exponential2.dat", v[[i]] ], {i,Length[v]} ]

	will do the trick.  The first argument of the Write function is
	the name of the file.  It must be included in double quotes.
	You may have additional arguments to specify what to write.  Here
	we only want a single element of the vector v.  Note the use of
	the function Length to determine the length of the vector v.  We
	thus don't need to explicitly remember the length of v in the
	iterator expression.

	We can use hist to make a histogram of the numbers.  If we use the
	-g 2 option to hist, and the logarithmic option to axis, we find
	a plot that looks rather linear, as we would expect with an
	exponential distribution.

	hist -n 20 -g 2 < exponential2.dat |axis y l 1 500 |xplot

	Box-Muller Method

	The Box-Muller method is a clever way to generate random
	numbers with a Gaussian Distribution.  It is described on page
	360 of CSM.  For homework you will get to try out this method and
	verify that the distribution is Gaussian.  You can also find
	this method described in Numerical Recipes.  We derived the Box-Muller
	method.

	Making Gaussian Random Numbers in Mathematica

	Here is an implementation of the Box-Muller method with sigma=1:

	gaus[] := Module[ {x,y,rho}, x=Random[];y=Random[]; rho = -Log[x];
			N[{Sqrt[2 rho]Cos[2.*Pi*y],Sqrt[2 rho]Sin[2.*Pi*y]}] ]

	gaus is now a function with no argument.  With the above definition, 
	you can just type
	gaus[]
	and you will get a list with two Gaussian random numbers in it.

	Of course, we want a lot to test the distribution, so try
	v = Table[ gaus[], {10000} ];

	Don't forget the semicolon at the end of the command!  Otherwise
	you are going to be looking at a lot of output.

	Check some values:
	v[[4]]

	Now export v to a file:
	Export["myfile", v, "Table"]
	This will result in a pair of numbers on each line of the file myfile.

	axis m 0 <myfile |xplot
	will produce a nice scatter plot.  To make the histogram, type

	~sg/bin/hist -n 25 -g <myfile |axis |xplot

	It is not too hard to use awk to continue to clean up the file.
	For instance, we can easily get awk to print each number on a
	separate line.  We can then use another awk command to print
	two numbers on a line.  The keys to this are the variable NR for
	the record or line number and NF for the number of fields, or words
	on each line.  Here is the awk command to make only one number on
	each line

	awk '{ for(i=1; i <= NF; ++i) print $i}' myfile >myfile.2

	To take a file with one number per line and then print pairs of numbers
	is also pretty easy.

	awk '{ if( NR % 2 == 1) x=$i; else print x,$1}' myfile.2 >myfile.3

	In Mathematica we could instead have used

	Do[Write["myfile", v[[i,1]], "  ", v[[i,2]]], {i, Length[v]}]

	to print two words per line.  Note that we really do need the
	string of blanks in double quotes or the two numbers will have
	no space between them.  How annoying!  It is good to learn the options
	of the Export function.

	The file that has two words per line will work for either the
	scatter plot or histogram.