Introduction to Monte Carlo Simulation

1. Introduction to Monte Carlo Simulation Robert L. Harrison University of Washington Medical Center Seattle, Washington, USA Supported in part by PHS grants CA42593 and CA126593

2. What is Monte Carlo? • Monte Carlo methods are a class of algorithms that rely on repeated random sampling to compute their results. (Wikipedia)

3. An example: estimating an integral • Randomly generate points in square. • Count the points under the curve. • Divide the number under curve by total number of points. • 100 points. • 69 under curve. • Estimate = 0.69. • 400 points. • 278 under curve. • Estimate = 0.695. • 1000 points. • 717 under curve. • Estimate = 0.717.

4. A little history… • Buffon’s needle. • Student / William Gosset. • Manhattan project.

5. Buffon’s needle • Georges-Louis Leclerc, Comte de Buffon (1707 – 1788). • A party game using a lined piece of paper and a pin. • Estimates the value of   = 2*drops/hits. http://www.shodor.org/interactivate/activities/buffon

6. William Gosset (1876-1937) • Supervised the experimental brewery at Guiness. • Developed methods to draw conclusions from small sample size experiments. • Published under the penname A. Student. • Student t-test.

7. Manhattan project • Development of first nuclear weapons. • Needed to extrapolate knowledge of a single neutron’s behavior to the behavior of a system of neutrons. • Experimentation was to dangerous, time-consuming, costly.

8. Manhattan project… • John von Neumann (1903-1957) and Stanislaus Ulam (1909-1984). • Developed particle tracking Monte Carlo simulation.

9. Why do people simulate?

10. Why do people simulate? • Low cost way to experiment. • Easy to repeat/modify. • Study systems to complex to solve analytically. • Examine hidden quantities.

11. Disadvantages • Time consuming. • Doesn’t give exact solution. • Is the model correct? • Bugs, bugs, bugs.

12. Alternatives • Analytic solution. • Experiment.

13. Mechanics of simulation • General method. • Random number generator. • Sampling from a distribution. • Accelerating a simulation.

14. General Monte Carlo Method • Model a physical system as a (series of) probability density function(s). • Repeatedly sample from the pdf(s). • Tally the statistics of interest. • (There are many Monte Carlo simulations that don’t fit the above exactly.)

15. Random numbers • Generally we sample from a uniform distribution of integers over [0,n]. • Convert to a uniform distribution on the interval (0,0), [0,1) or [0,1]. • (e.g., divide integer value by n or n+1.) • Lists of truly random numbers. • E.g., time between radioactive decays. • Pseudo-random number generators. • Calculate a randoms number!?

16. Pseudo-random number generators • Usually called random number generators (RNG). • (Early - bad) Middle squares. • (Still around - bad) Linear congruential generator. • (Good) Mersenne twister.

17. Middle square • Von Neumann, 1946. To generate a sequence of 10 digit integers, start with one, and square it and then take the middle 10 digits from the result as the next number in the sequence • e.g. 57721566492 = 33317792380594909291 the next number is • Not a recommended method!

18. Linear congruential generator • Lehmer, 1948. In+1 = (aIn + c) mod m • I0 = Starting value (seed) • a, c, and m are specially chosen constants.

19. Linear congruential generator… • Values used by some built-in random functions are very poor. • Good values for a, c, and m on-line. • Even with good values, there can be problems. Three dimensional sampling from a linear congruential generator with poorly chosen values, RANDU (Wikipedia).

20. Mersenne twister • Used in MatLab, Maple. • Does very well on most tests of randomness. • Ridiculously long period. • Reasonably fast.

21. Mersenne twister • Used in MatLab, Maple. • Does very well on most tests of randomness. • Ridiculously long period. • Reasonably fast. • Quite complicated to implement/describe.

22. Sampling from a distribution • We have a uniform distribution, how do we sample from • An exponential distribution? • Other continuous distributions? • A Poisson distribution?

23. Exponential distribution • Time between radioactive decays. • Distance a photon travels before an interaction. • Mean 1/, variance 1/ Probability density function (pdf) (Wikipedia) Cumulative distribution function (cdf)

24. Sampling from the exponential distribution • Sample u from the uniform distribution on [0,1).

25. Sampling from the exponential distribution • Sample u from [0,1). • Locate u on the y-axis of the exponential distribution cdf. u

26. Sampling from the exponential distribution • Sample u from [0,1). • Locate u on y-axis. • Find the x value that corresponds to u. u x

27. Sampling from the exponential distribution • Sample u from [0,1). • Locate u on y-axis. • Find x = F-1(u). u x

28. Sampling from an invertible continuous distribution • Sample u from the unit interval. • Find x = F-1(u).

29. Other sampling methods for continuous distributions • Approximate F-1(u), e.g., as a piecewise linear function. • Table lookups are fast! • Approximate => bias. • Acceptance-rejection method. • AKA rejection method, accept-reject method, rejection sampling.

30. Acceptance-rejection • To sample from a pdf f(x): • Choose pdf g(x) and constant c such thatc*g(x) ≥ f(x) for all x www.mathworks.com

31. Acceptance-rejection • To sample from f(x): • Choose g(x) and c: c*g(x) ≥ f(x) for all x. • Generate a random number v from g(x). • Generate a random number u from the uniform distribution on (0,1). v

32. Acceptance-rejection • To sample from f(x): • Choose g(x) and c: c*g(x) ≥ f(x) for all x. • Generate v from g(x). • Generate u from uniform distribution. • If c*u ≤ f(v)/g(v) accept v • else reject v and go to 2.

33. Acceptance-rejection • To sample from f(x): • Choose g(x) and c: c*g(x) ≥ f(x) for all x. • Generate v from g(x). • Generate uniform u. • If c*u ≤ f(v)/g(v) accept v • else reject, go to 2. • The average number of iterations is c.

34. Acceptance-rejection • Sampling from a normal distribution. • While the cdf is invertible, the resulting formula is computationally inefficient. • Generally sampling is done using the Box-Muller or Ziggurat algorithms. • Use a library function! • Acceptance-rejection also be used to remove the bias from an approximation/table lookup sampling method.

35. Sampling from a normal distribution • Gaussian, bell curve. • Everything is normal. • Central limit theorem: the sum of a large number of independent random variables is approximately normally distributed. • Generally sampled using the acceptance-rejection method. • Use a library function! pdf cdf (Wikipedia)

36. Poisson distribution • Discrete probability distribution. • {0,1,2…} (non-negative integer) result. • Mean = variance =  Probability mass function (Wikipedia) Cumulative density function

37. Poisson distribution • Distribution for: how many light bulbs will burn out in a month? • Large building. • Maintenance department immediately replaces burned out bulb. • Choose  = the average number that burn out in a month How many mathematicians does it take to change a light bulb? In earlier work, Haynor[1] has shown that one mathematician can change a light bulb. (www.clipartof.com)

38. Poisson distribution • Other examples: • Number of radioactive decays. • Number of bugs discovered in a program. • Number of calls to a call center. • Number of scintillation photons. How many mathematicians does it take to change a light bulb? None. It's left to the reader as an exercise. (www.clipartof.com)

39. Sampling from Poisson distributions • Use a library function! • Use the exponential distribution. • Let T be the time period, set t = 0, k = 0. • Sample tk from the exponential distribution. • t = t + tk, k = k + 1; • If t < T, go to 2. • Return (k - 1). • For  large approximate using a normal distribution. How many quantum physicistsdoes it take to change a light bulb? They can't. If they know where the socket is, they cannot locate the new bulb. (www.clipartof.com)

40. Keeping track of results • Tally (score, bin, histogram) the quantities we are interested in, e.g.: • Number of light bulbs replaced by month; • Time to failure for a product.

41. Building a simulation • Often a single ‘event’ in a simulation will require sampling from many different pdfs, some many times, and be tallied in many ways.

42. Building a simulation • For example, a nuclear physics simulation might model • the time between decays • the decay process • the paths of photons and particles • the energy deposited during tracking and any secondary particles created thereby • the detection electronics… • and tally • energy, position, and type of detected particles.

43. Take away • Simulation a key tool for analytically intractable problems. • Currently the Mersenne Twister is a good choice of random number generator. • Random sampling for functions: • cdfs are analytically invertible; • use library functions if available; • acceptance-rejection method; • piecewise approximation. • Four key distributions: uniform, exponential, normal, Poisson.

45. References J.T. Bushberg, The essential physics of medical imaging, Lippincott Williams & Wilkins, 2002. K.P. George et al, Brain Imaging in Neurocommunicative Disorders, in Medical speech-language pathology: a practitioner's guide, ed. A.F. Johnson, Thieme, 1998. D.E. Heron et al, FDG-PET and PET/CT in Radiation Therapy Simulation and Management of Patients Who Have Primary and Recurrent Breast Cancer, PET Clin, 1:39–49, 2006. E.G.A. Aird and J. Conway, CT simulation for radiotherapy treatment planning, British J Radiology, 75:937-949, 2002. R. McGarry and A.T. Turrisi, Lung Cancer, in Handbook of Radiation Oncology: Basic Principles and Clinical Protocols, ed. B.G. Haffty and L.D. Wilson, Jones & Bartlett Publishers, 2008. R. Schmitz et al, The Physics of PET/CT Scanners, in PET and PET/CT: a clinical guide, ed. E. Lin and A. Alavi, Thieme, 2005. W.P. Segars and B.M.W. Tsui, Study of the efficacy of respiratory gating in myocardial SPECT using the new 4-D NCAT phantom, IEEE Transactions on Nuclear Science, 49(3):675-679, 2002. V.W. Stieber et al, Central Nervous System Tumors, in Technical Basis of Radiation Therapy: Practical Clinical Applications, ed. S.H. Levitt et al, Springer, 2008. P. Suetens, Fundamentals of medical imaging, Cambridge University Press, 2002. www.wofford.org/ecs/ScientificProgramming/MonteCarlo/index.htm www.shodor.org/interactivate/activities/buffon www.impactscan.org/slides/impactcourse/introduction_to_ct_in_radiotherapy