Lab Assignment: Monte Carlo Techniques

This lab is to be completed mostly on a computer. You can use the python tools you installed on your PC or Mac. In this lab there are four related parts all designed to introduce you to Monte Carlo techniques, essential tools in the analysis of data in all modern scientific endeavors, build familiarity with the three basic statistical distributions important in data analysis, and to further your computing skills.

Part A. Calculate the value of π (pi) using Monte Carlo methods

In this exercise, you will write a program that will be the basis of the rest of the lab. In it, you will use the acceptance/rejection Monte Carlo method to calculate the value of π and its uncertainty. The number of iterations should be sufficient to provide an accuracy (relative error dpi/pi) of π o at least 0.1%. You will need to use the appropriate number of interactions to achieve this precision.

• You can create your program either on your PC or Mac.
• The computer program can be written in any language you want as long as it is FORTRAN, C, C++, Python, or MATLAB.
• If you write a C, C++, or FORTRAN program you will need to compile your code. You do this by issuing the appropriate compiler command "g++" for a C or C++ program, "gfortran" for FORTRAN.
• An outline of the procedure is described below. The exact prescription is detailed in my places, online as well as in Bevington & Robison. Basically, the idea is to compare the area of a circle to the area of a square. This ratio is proportional to π and we use this fact to get π from an MC estimate of the ratio. The pseudo-code that will accomplish this is:
  1. Throw (generate) two random numbers "x" and "y" with a random number generator. In c++ for example this can be done by calling the function rand(). It is available in the standard C++ library stndlib. Make sure the number returned in the interval of [0,1].
  2. Next check if x,y point falls within a circle by the condition if (x**2+y**2 < 1) then increment a counter (integer) by one.
  3. Compute π by taking the ratio of your counter to the number of iterations. As this is proportional to π you have completed the computation of the central value of π.
  4. Finally, calculate the statistical uncertainty of π given that the process used to determine π is a binomial one the uncertainty will thus be proportional but not equal to the binomial uncertainty.

• Also part of the assignment is to find the number of events needed to obtain statistical uncertainties of 1%, 0.01%, and 0.00001%. One way to do this is the brute force way and use the program changing the number of iterations until you get the desired result. This could take a while at least for the smallest uncertainties. There is an analytic method that involves solving for N using the value of pi from your calculator. Make sure you include these estimates in your paper.

Part B. Calculate the Volume of N-dimensional spheres

Modify your program, or write a new one, to determine the volumes of N-dimensional spheres (N=1,2,3,4,5 but you've already done the 2D one) with radius=1.0. Use 10,000 tries and make sure you quote your uncertainties. You should also compare your results with analytic values that you can compute by some other means. Please include the means used to compute the analytic values.

Part C: Generate Poisson deviates, histogram, and plot the distribution

In Part C the assignment is to generate a sample of Poisson deviates using Monte Carlo methods, specifically the transformation method. The method depends on the conservation of probabilities which stipulates that a draw from a uniform random distribution is equivalent to a draw from any other distribution, including the Poisson. Symbolically,

where is the uniform distribution defined on the interval and is, in this case, the Poisson distribution function, a discrete distribution,

with x an integer and mean μ.

Here is how to do it:

1. draw a random number r from a uniform distribution by calling the random() function
2. now compute the Poisson probability, the argument of the summation above, for integers x starting at 0 and continuing on to some n. Each time add terms to form the summation as required by the equality above.
3. Compare the summation as the integers increase and compare to the random number r drawn from the uniform distribution above.
4. Once you have a summation that is greater than the r, which corresponds to the LHS, write the integer in the previous iteration x into a table and that will be your Poisson deviate that corresponds to r,

Repeat this process N = 10,000 times to obtain 10,000 Poisson deviates.

You will need as an input parameter the mean of the distribution. You will then histogram the deviates and plot the histograms, include these in your report. Remember that histogram binning is important. If you don't get this right the plots will look funny, either there will be spaces between bins of some bins will be much larger than their neighbors, indicating overflow. Given that you are plotting Poisson distributions that are discrete the bin width is = 1. Please repeat the generation of Poisson deviates 3 times: one with μ = 1.0, μ = 10.3 and μ = 102.1.

Part D. Generate and Fit a Gaussian

In this exercise, you employ the Monte Carlo technique to generate data distributed as a normal or Gaussian distribution. Then you will fit the generated data sample to a Gaussian Probability Distribution Function (PDF), determine the fit parameters and their errors, and discuss the goodness of the fit via the reduced chi2.

There are many ways to generate a normally distributed set of data. One way is to start off with a uniform (flat) distribution of data points, generated between 0 to 1, and transform them to a normally distributed set of data via the Box-Muller 2D Gaussian transformation described in Bevington and Robinson or here<a href:http://pdg.lbl.gov/2011/reviews/rpp2011-rev-monte-carlo-techniques.pdf > http://pdg.lbl.gov/2011/reviews/rpp2011-rev-monte-carlo-techniques.pdf </a> This clever algorithm is both efficient and easy to implement. There are many sources you can read about it, there is even a Wikipedia article on it. There is an easier way but one I do not want you to use. It turns out that given the usefulness of having random numbers drawn from a Gaussian area, every programming language has normal or Gaussian deviate generators in standard libraries. For example, numerical Python has a built-in normal random number generator in NumPy and probably all other math libraries. Again please do not use a built-in method or any functions that invoke a normal random number generator. The assignment calls for you "coding up" the Box-Muller procedure yourself.

Here is how to do it:

• Use uniform random number generators to generate random numbers drawn from the Gaussian centered at 0 with a variance of 1.0 using the Box-Muller transform as described in the references above.
• Bin the transformed random number data into 40 bins from -4.0 to + 4.0. 1.0.
• Fit the histogram generated above with a Gaussian function, or Gaussian PDF.

• Find the parameters of the fit and quote their uncertainties. Do you know how these were determined? Feel free to explore how fitting package work in general, and how they find fit parameters. Quote the parameters and their uncertainties prominently in your lab report. You should discuss how these were determined by the fitter. You should look through Bevington and Robinson for background and feel free to explore how the functions used in Python actually do the fits (this may take a bit of digging around Dr. Boeglin's LT.box code since those tools wrap standard Python curve fitting libraries found in scipy and maybe matplotlib. Learning about this is something everyone do even if it is done at a superficial level.

• Determine the goodness of fit. You can use the reduced chi2 if your fitter uses the least-squares method to determine the fit parameters. Also, discuss in your lab report the use of the reduced chi2, what it means and why can we use it as a means by which to estimate the "goodness of fit".

GradingRubric