PHYCS 6730 Lab Exercise: Monte Carlo Simulation

These exercises play with elementary simulations that lead to some frequently encountered probability distributions. The example files are in ~p6730/exercises/monte.

Exercise 1. Poisson Distribution

The Poisson distribution arises in many applications in science and engineering. For example in radioactive decay in which the sample is not significantly depleted so the average decay intensity is constant, the number of decays occuring over a fixed time interval follows a Poisson distribution. In any similar process in which the event rate is constant, but the events are random and uncorrelated, the number of events in a given time interval follows a Poisson distribution.

Our goal in this exercise is to simulate a radioactive decay process, using the standard library random number generator rand(), which gives numbers in the range [0,RAND_MAX] where the upper bound is defined in stdlib.h. The radioactive decay formula gives

   dN/dt = -k N   or   p = dN/N = -k dt

where N is the number of atoms at a given time t, k is the decay constant, and p is the probability for decay over the time interval dt. Since we are assuming the intensity is constant, p should be small (so only a tiny fraction decays) to get a good approximation to a Poisson process. (If it is too large, we will get the binomial distribution instead.)

You are given the decay probability p and the number of atoms N. Think about how you would construct the simulation.

We can simulate the decay process by examining each of the N atoms at a time and asking whether it decays. How could we use the random number generator to decide whether a particular atom decays? How would we determine the number of decays of N atoms in a given simulation?

To see the Poisson distribution, we must run the simulation lots of times (say 1000) to generate enough results to construct a histogram.

The code genpoi.cc does the simulation. Run 1000 simulations for p = 0.1, N = 50 and construct a histogram of counts. Compare your results with the expected Poisson distribution (mean value pN = 5):

   P(k) = 5^k exp(-5)/k!

Hint: Since the counts are integers, construct the histogram bins so they have width 1 and divisions between the integers.

Report your findings in your answer file. Don't include the 1000 random numbers.

Exercise 2. Normal distribution

Here our task is to generate Gaussian normal deviates z with mean z0 and standard deviation sig. To do this, it suffices to generate x with mean 0 and standard deviation 1 and then do the linear transformation

   z = x * sig + z0

To get x from rand() we use a trick. While the distribution

    P(x)dx = 1/sqrt(2*pi) exp(-x²/2)dx

can not be integrated in terms of elementary functions, the two-dimensional distribution

    P(x,y) dx dy = P(x)P(y) dx dy = 1/(2*pi) exp(-r²/2) r dr dtheta

can be integrated in cylindrical polar coordinates. If we can generate a random r and theta according to the above distribution, we can obtain two random Guassian deviates x = r cos(theta) and y = r sin(theta).

The code gengau.cc does this. Use it to generate 1000 random deviates with mean -1 and standard deviation 2. Then use meansd to check the population mean and standard deviation to see if they agree with your design.

Report your findings in your answer file. Don't include the 1000 random numbers.

Exercise 3. Chi square distribution

The chi square distribution gives the probability for getting

   x² = x₁² + x₂² + ... + x_df²

when x_i are independent Gaussian normal deviates with zero mean and unit standard deviation.

Given a procedure for generating x_i, how would you generate x²?

The code getchi.cc does the job. Generate 1000 values of x² with df = 4 and find the probability of chi square exceeding 6.2. Compare with the value given by the course utility conf_level, which evaluates the probability by deterministic methods (incomplete gamma function).

Hint: To count the number of times chisq = 6.2 is exceeded, redirect the genchi output to a file and try using the "awk" utility to print only the lines with chisq > 6.2. The Unix utility "wc" counts lines.

Report your findings in your answer file. Again, don't include the 1000 random numbers.