Monte Carlo Method Fundamentals PDF

Part II Fundamentals of Monte Carlo Method Samira Abbaspour, PhD Department of Medical Physics, Tehran University of Medical Sciences, Tehran, Iran Outline ELEMENTARY PROBABILITY THEORY RANDOM NUMBER GENERATOR SAMPLING THEORY PHOTON TRANSPORT I...

Part II Fundamentals of Monte Carlo Method Samira Abbaspour, PhD Department of Medical Physics, Tehran University of Medical Sciences, Tehran, Iran Outline ELEMENTARY PROBABILITY THEORY RANDOM NUMBER GENERATOR SAMPLING THEORY PHOTON TRANSPORT IN MEDIA INTERACTION MODELING Elementary Probability Theory Fundamental to understanding the operation of a Monte Carlo process and interpreting the results of Monte Carlo calculations, is some understanding of elementary probability theory. A probability distribution function (pdf), p(x) is a measure of the likelihood of observing x. For example, x could be the position at which a photon interacts via the Compton interaction.  Elementary Probability Theory In general, p(x) has some special properties that distinguish it from other functions: p(x) ≥ 0 since negative probabilities have no meaning p(x) is normalized in the following fashion, If the variable is discrete −∞ < xmin < xmax < +∞, that is, xmin and xmax can be any real number so long as xmin is less than xmax. If the variable is continuous. It is also possible to create a probability density function from a distribution function x(f) which is not normalized. Elementary Probability Theory Associated with each probability distribution function is its cumulative probability distribution function (cpdf) Cumulative probability distribution functions have the following properties which follow directly from its definition and the properties of probability distribution functions: p(x) and c(x) are related by a derivative, c(x) is zero at the beginning of its range of definition, and unity at the end of its range of definitions, c(x) is a monotonically increasing function of x as a result of p(x) always being positive Elementary Probability Theory pdf cpdf Outline WHAT IS MONTE CARLO ELEMENTARY PROBABILITY THEORY RANDOM NUMBER GENERATOR SAMPLING THEORY PHOTON TRANSPORT IN MEDIA INTERACTION MODELING VARIANCE REDUCTION TECHNIQUES Random Number Generation “Anyone who considers arithmetical methods of producing random digits is, of course, in a state of sin” John von Neumann (1951) The “pseudo” random number generator (RNG) is the “heartbeat" of a Monte Carlo simulation. The operative phrase to be used when considering RNG's is “ use extreme caution ". DO USE an RNG that is known to work well and is widely tested DO NOT FIDDLE with RNG's unless you understand thoroughly the underlying mathematics and have the ability to test the new RNG thoroughly. DO NOT TRUST RNG's that come bundled with standard mathematical packages Random Number Generation Mathematically speaking, the pseudo-random numbers used in Monte Carlo simulation should possess the following properties: Uncorrelated sequences: the sequences of random numbers should be serially uncorrelated. Long period: ideally, the generator should not repeat; practically, the repetition should occur only after the generation of a very large set of random numbers. Uniformity: the sequence of random numbers should be uniform, and unbiased. Reproducibility: when debugging programs, it is necessary to repeat the calculations to find out how the errors occurred. The feature of reproducibility is also helpful while porting the program to a different machine. Speed: It is of course desirable to generate the random numbers fast. Parallelization: The generator used on vector machines should be vectorizable, with low overhead on massively parallel architectures. Test of Uniformity The uniformity of RAND function in MATLAB 50000000 50 500 5000 50000 500000 5000000 Random Number Generators http://random.mat.sbg.ac.at http://www.heparts.com/products Outline WHAT IS MONTE CARLO ELEMENTARY PROBABILITY THEORY RANDOM NUMBER GENERATOR SAMPLING THEORY PHOTON TRANSPORT IN MEDIA INTERACTION MODELING  Sampling Theory Now that we have tackled the essentials of elementary probability theory and random number generation, it is now time to connect the two and demonstrate how random numbers may be employed to sample from probability distributions. We will consider three kinds of sampling techniques: Direct Method Rejection Method Mixed Method Analog Sampling Library function Markov chain  Direct Method Suppose we have a typical pdf as follow: It is define over a range of [a,b] wher neither a nor b Are necessarily finite. Then we construct its cumulative probability function as follow: This method is useful when the cumulative probability function is invertible By its definition, we can map the cumulative probability distribution function onto the range of random variables, r, where 0 ≤ r ≤ 1 and r is distributed uniformly.  Rejection Method While the invertible cumulative probability distribution function method is always possible, at least in principle, it is often impractical to calculate c( ) -1 because it may be exceedingly complicated mathematically or contain mathematical structure that is difficult to control. Another approach is to use the rejection method. In recipe form, the procedure is this: 1. Scale the probability distribution function by its maximum value obtaining a new distribution function, f(x) = p(x)/p(xmax), which has a maximum value of 1 which occurs at x = xmax. Clearly, this method works only if the probability distribution function is not infinite anywhere and if it is not prohibitively difficult to determine the location of the maximum value. If it is not possible to determine the maximum easily, then overestimating it will work as well, but less efficiently.  Rejection Method 2. Choose a random number, r1, uniform in the range [0; 1] and use it to obtain an x which is uniform in the probability distribution function's range [a,b]. (To do this, calculate x = a + (b − a)r1.) (Note: This method is restricted to finite values of a and b. 3. Choose a second random number r2. If r2 < p(x)/p(xmax) (region under p(x)/p(xmax) ) then accept x, else, reject it (shaded region above p(x)/p(xmax)) and go back to step 2. The efficiency of the rejection technique is defined as: This is the ratio of the expected number of random numbers pairs that are accepted to the total number of pairs employed. However, it can save computing time if the c( )−1 is very complicated. One has to “waste" many random numbers to use as much computing time as in the evaluation of a transcendental function! Mixed Method Mixed method is the combination of direct method and rejection method Imagine that the probability distribution function is too difficult to integrate and that it is “spiky", rendering the rejection method inefficient. (Many probability distributions have this objectionable character.) However, imagine that the probability distribution function can be factored as follows: where f(x) is an invertible function that contains most of the “spikiness”, and g(x) is relatively flat but contains most of the mathematical complexity. The recipe is as follows: Library function Method Default functions Questions What distribution function do we want? What parameters does it have? How many random numbers do we want?  Markov Chain Monte Carlo (MCMC) It is suitable for functions that have a very complex distribution. Due to the complexity of these functions, their inverse cannot be drawn. This method is a powerful method in statistics to generate random numbers, but it is not used much in medicine. Thanks For Your Attention