November 11, 2005 --- Class 21---Second Semester Outline,
			Inverse Transform Method

Activities:

	Second Semester Outline

	Below is what I planned to do last year.  What I actually did was
	a little less ambitious.  This Spring, assuming enough people sign
	up to offer the class, we can discuss what topics to cover.

Fractals, Kinetic Growth Models & Cellular Automata     CSM Ch 13       3 lec
3 D Graphics with Mathematica                           notes           1 lec
Fields of Static Charges and Currents                   CSM Ch 9        2 lec
Microcanonical Ensemble                                 CSM Ch 15       3 lec
Canonical Ensemble                                      CSM Ch 16       3 lec
Quantum Systems                                         CSM Ch 17       3 lec
Partial Differential Equations                          NR 19           3 lec
Fast Fourier Transform                                  NR 12,13        2 lec
Linear Algebra                                          NR 2,11         3 lec
Performance Analysis (prof, gprof, SvPablo)             notes           2 lec
Introduction to Parallel Computing                      notes           3 lec

	Next homework

	It is convenient use the split utility, or my blockavg program for
	the sample mean problem in the next homework on Prob. 11.6.  
	We discussed how these may be used.

	Inverse Transform Method

	Sometimes we are interested in generating random numbers with a non
	uniform probabilty density.  There are a number of ways of doing
	this and we shall study a few.  This is covered in CSM section
	11.5, starting on page 358.  The first method we will study is
	called inverse transform.  You may be used to the idea of the
	probability density function.  For a real variable, we have
	$\int_{-\Infinity}^{+\Infinity} p(x) dx = 1 $.

	We can define the monotonically increasing function called the
	cumulative probability distribution by
	$ P(x) = \int_{-\Infinity}^{x} p(x') dx' $.

	Of course, $ d P(x)/dx = p(x)$.  P(x) is a number between 0 and 1.
	If we can solve the equation $ P(x)=r$  for $x(r)$, then if we
	generate uniform random numbers r, then the numbers x(r) have the
	distribution we want.  So, the key is we start with p(x), and
	the method only is useful if we can find P(x) and then solve
	P(x)=r for x.  One useful example where this works in a exponential
	distribution.  We were able to show that 

	P(x) = 1 - exp(-x/lambda)
	
	Setting P(x)=r and solving for r, we find

	x = - lambda * ln(1-r) 

	Since 1-r and r are both uniformly distributed between 0 and 1,
	we can simplify the above and replace ln(1-r) by ln(r).
	Thus, to generate numbers with an exponential distribution, we
	first generate a uniform random variable between 0 and 1.  Then
	
	x = -lambda * ln(r)

	has an 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.)

	We can make use of Mathematica Histogram function:

	Needs["Graphics`Graphics"]
	Histogram[v]

	You can click on the graph with the right mouse button and find
	the appropriate item to save it.  Alternatively, use the Display
	command:

	Display[ "file", Histogram[v], "EPS" ]
	
	will save the graph is the file named "file" in Encapsulated
	PostScript Format.  Other formats such as PDF, TIFF and GIF are
	available.  See the documentation for additional options.

	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.