Basic notions about Monte Carlo method

---

Index:

• Introduction
• Intagration via MC
• Sampling methods
• Random Walk

---

Introduction

With Monte Carlo method are indicated a bunch of stochastic procedures which exploit large sample of random numbers.

A classic example is the evaluation of \pi : Imagine a circle with unitary radius R=1 inside a square of length 2R. Now let’s sample uniformly a certain number of points. These points will fall inside the circle with a probability which is proportional to the fraction of the two areas: P =\frac{\pi R^2}{(2R)^2} . In this case the meaning of probability is: Number of points in circle (Nc) over the total number of points (Nt), then \frac{Nc}{Nt} = \frac{\pi R^2}{(2R)^2} = \frac{\pi}{2}.

So the MC procedure is:

• Sample (uniformly) a certain number or position inside the square;
• Count the points inside the circle;
• Evaluate the fraction Nc/Nt

• 

• 

• 

It’s obvius that more points mean higher accuracy, in particular the results from these three tests are:

Number of sample	\pi
1e3	3.1120
1e4	3.1456
1e5	3.1458
5e5	3.1404

Code in Matlab:

close all
clear

x = 2*rand(500000,1)-1;
y = 2*rand(500000,1)-1;

r = sqrt(x.^2 + y.^2);

ID_in = find(abs(r)<1);
ID_out = find(abs(r)>=1);

x_in = x(ID_in);
y_in = y(ID_in);

x_out = x(ID_out);
y_out = y(ID_out);

plot(x_in,y_in,'r.')
hold on
plot(x_out,y_out,'b.')

PI_greek = length(x_in)/length(x)*4

Integration via MC

Through the sampling of random numbers is possible to evaluate complex integrals. let’s consider a generic function g(x) with x \in (a,b).

CASE 1: The integral we want to evaluate is G = \int g(x)dx. It is always possible to apply the mean value theorem: G = \bar g \int dx = \bar g(b-a). Now we can use MC method, since we can evaluate the average value of g by sampling from an uniform distribution x_i in the domain (a,b): \bar g = \frac{1}{N} \sum_i g(x_i)

CASE 2: Now let’s consider a sligthly different situation: our integral now is: G = \int g(x)f(x)dx Similarly as before we can say G = \bar g\int f(x)dx and so the average value of g is \bar g =\frac{ \int g(x)f(x)dx}{\int f(x)dx} where \int f(x)dx is the normalization factor. From this we can understand that f(x) can be seen as a probability distribution from which we can extract our samples and evaluate the integral. (case 1 is the particular case in which f(x) is uniform and equal to 1)

Sampling methods

The sampling from a distribution is the basic of Monte Carlo methos, several methods exist, here we will see two of them.

Inverse tranform method:

Let’s call f(x) our normalized distribution function and CDF(x) the cumulative distribution fuction defined as CDF(x) = \int_{-\infty}^x f(x')dx'.

Inverse transform method consist in sample F from an uniform distribution (between 0 and 1), calling the sample Rand we have Rand = CDF(x) = \int_{-\infty}^x f(x')dx', then inverting the formula: CDF^{-1}(Rand) =x.

Example: sampling from f(x) = x with x \in (0,1)

CDF(x) = \int_{-\infty}^x x' dx' = \frac{x^2}{2}

Now we sample a random number from an uniform distribution, let’s call it Rand

Rand = \frac{x^2}{2}

Finally our sample for x is:

x = \sqrt{2Rand}

This example is the same as “sampling uniformly inside a circle”. In fact let’s imagine to put a point randomly in the circle, the probability of having the radius in the range [r, r+dr] is rapresented by the fraction of the areas: \frac{2\pi rdr}{\pi R^2} this mean that the distribution we have to adopt in order to sample a point uniformly in a circle is: f(r) = \frac{2 r}{ R^2} quite similar to the one in the example. In polar coordinate the other quantity we have to evaluate is the angle \theta which can be sampled uniformly from zero to 2\pi.

Matlab script example:

Rejection method:

If our probability distribution function is particularly complex to integrate we can adopt the Rejection method. As usual our p.d.f. is named f(x) (blue line), let’s consider an uniform distribution with value H higher than the maximum value of f(x). The procedure consist in sampling a tentative value R uniformly in (a,b), then the acceptance probability is p_{acc} = \frac{f(R)}{H}, in the figure is rapresented by the green line. Now we sample another random number (between 0 and 1) and we compare it with the acceptance probability, if it is smaller than p_{acc} , tentative value is accepted, if not is discarded and another tentative value is sampled.

Closest is H to the maximum value of the p.d.f. higher is the efficiency of the method.

Matlab example:

Random walk

MC method is often used to track the trajectory of particles, imagine a neutral particle moving in a material, it will move along a straight line until it collide changing its direction. This motion, if we consider a large number of particles, can be seen as a sequence of random numbers (Sampled from an appropriate distribution) and it’s called random walk.

Let’s align our reference frame with the free flight direction of the particle, calling \Sigma_t the total macroscopic cross section (which is a probability of interaction per unit length) we have that the probability of non-interaction in a length x is e^{-\Sigma_t x}, instead the probability of interaction in dx is \Sigma_t dx, then the probability of interaction between x and x+dx is obtained by the probability of non interaction before x multiplied by the probability of interaction in dx: e^{-\Sigma_t x}\Sigma dx. So finally our pdf for the sampling of the free flight is:

f(x) = \Sigma_t e^{-\Sigma_t x}.

Free flight (let’s call it d) can be easly sampled using the Inverse transform method:

d = -\frac{1}{\Sigma_t} ln(1-Rand)

After the collision, the particle will have a different direction, in polar/spherical coordinate, a different angle. The distribution from which the angle must be sampled depend on the particular situation, the simplest approach is to consider an uniform distribution for the solid angle. This mean that:

Cos\theta = 2Rand - 1, \phi = 2\pi Rand

More on classical collisions