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.