Skip to content

5. Monte Carlo Methods

5.1 Random numbers

Consider the output of the following code:

random module
import random            

random.seed(1000)        
print (random.random(), random.random())

The output will be two numbers printed on the screen. The random.random() command gives you a random number in the range [0.0, 1.0). The likelihood of getting a particular value between 0 and 1 is the same. In statistical language, such a likelihood/probability is called a uniform distribution. If you want to get a random number in a different range, say [a, b], you can either use the result of random.random() and rescale it appropriately, or directly use the function random.uniform(a, b).

Use of random.seed(): If you use the same seed number (argument) in calling random.seed(seed) before calling random.random(), the output is identical for the same seed. Because the returned random numbers depend on the seed, the number generator used by the random module is called a pseudo-random number generator. Such a generator is useful if you have to reproduce results (e.g. simulation) at a later time. If you do not provide an argument to random.seed or set random.seed(None), the current system time is used as a proxy to a random seed and you will get different results each time you call random.random().

5.1.1 Estimating area using random numbers

We can make use of random numbers to estimate areas. Suppose we want to estimate the area of the unit circle.

Area of a unit circle
\[ A = \pi r ^2 = \pi (1)^2 = \pi \]

In this case, we know the answer is \(\pi\). Therefore, we can frame this exercise as a method to estimate the numerical value of \(\pi\) as well.

Consider a unit circle embedded in a square. Now, in each step, we will generate two random numbers \((x, y)\) inside the embedding square, and check whether this randomly generated point lies inside of the unit circle or outside.

Then an estimate of the area of the unit circle can be made:

\[ A = \left( \frac{\rm points\ inside\ the\ circle}{\rm total\ points} \right) \times {\rm area\ of\ embedding\ square}\]



Total Points:

Points Inside:

Area Estimate:

Refresh the page to estimate using a new set of random points.

Python code to estimate \(\pi\)

Estimate area of a unit circle
def pi_estimate(total=1000):
    pts_in = 0

    for i in range(total):
        x = random.uniform(-1, 1)
        y = random.uniform(-1, 1)

        if (x**2.0+y**2.0 <= 1.0):
            pts_in = pts_in + 1.0

    return 4.0*pts_in/total

print (pi_estimate())

Increasing the number of points used through the total argument in the pi_estimate function improves the estimate. If \(N\) is the number of points used, the standard error scales as \(1/\sqrt{N}\).

5.2 Monte Carlo Integration

5.3 Markov Chain Monte Carlo (MCMC)

5.4 Examples