**Python Monte Carlo Simulation Example**

**Introduction**

As Python becomes more widely accepted and the Python community broadens, it becomes more difficult for one side of the community to see what the other is doing. In response to comments we received about Python only being useful as a web programming language, we decided it was time to write about the Python random number module, and use it for a couple of simple Monte Carlo simulations.

**The random module**

The random number module may be imported into a Python program by typing “import random.” Once imported, it provides a number of useful functions. The following text is taken from the Python documentation: **Numpy Monte Carlo**

choice (seq): Choose a random element from the non-empty sequence seq and returns it.

randint(a, b): Returns a random integer N such that a <= N <= b.

random(): Returns the next random floating point number in the range [0.0 … 1.0).

uniform(a, b): Returns a random real number N such that a <= N < b.

In this article we’ll mostly be focusing on the random() function. **Scipy Monte Carlo**

When you import the random module, Python uses a Wichmann-Hill pseudo-random number generator. All computer random number generators are generally pseudo-random number generators, because software programs can only generate nearly random sequences. Since this is an article on Python, not on random number generation, we won’t go further into the topic, except to say that the Wichmann-Hill algorithm is thought to be pretty good. People who have high-powered, scientific needs for random numbers might want to check out the RNG random number generator package from Lawrence Livermore National Laboratories available here.

**Calculating pi**

So what, then, can you do with a random number generator? One of the more useful applications of random numbers in scientific computing is something called Monte Carlo simulation. The name refers to the gambling city, and denotes that the technique is based upon random processes, such as dice rolls.

**Monte Carlo Simulation Python Code**

We can write a Python function to simulate six-sided dice as follows:

import random

def six_sided_dice():

return random.randint(1,6)

Fairly simple. The first computer program we ever wrote was a program to roll different sided dice for Dungeons and Dragons, and it is straightforward to extend this function to arbitrary-sided dice, the 4, 8, 12, and 20 sided dice that are used in D&D;, as well as 29-sided-dice, should the mood strike you. Most Monte Carlo processes, however, deal with random floating point numbers rather than random integers, so we’ll find the random() function more useful than randint().

A Monte Carlo algorithm for calculating pi (the ratio of a circle’s circumference to its diameter) comes from the observation that the ratio of the areas of a circle and its circumscribing square is pi/4. The figure below shows this in more detail.

We can make use of the ratio of areas to compute pi using a random number generator. Here is a Python code snippet for this computation:

import random

npoints = 100

total_shots = 0

total_hits = 0

for i in range(npoints):

ax = random.random()

ay = random.random()

a2 = ax*ax + ay*ay

total_shots = total_shots + 1

if a2 < 1:

total_hits = total_hits + 1

pi_approximation = 4.0*total_hits/total_shots

print i, pi_approximation

Given enough time and a good enough random number generator, this algorithm will produce a decent value of pi. After 100 iterations we get a value of ~3.3, after 1000 we get a value of 3.144, and after 10000 iterations we have a value of 3.122. Hmmmm.

This is a really fabulous method of computing pi, except when you compare it to any other technique. A much better way of computing pi is to use Machin’s formula:

A Python code snippet to compute pi using Machin’s formula is

npoints = 10

pi = 0.

for i in range(npoints):

pi = pi + 16*pow(-1,i)*pow(0.2,2*i+1)/(2*i+1) \

– 4*pow(-1,i)*pow(1./239.,2*i+1)/(2*i+1)

print pi

After only 10 iterations with Machin’s technique we obtain a value of 3.14159265359. Here is a comparison of the convergence of the two methods for the first 50 steps:

*Monte Carlo Python pi*

Monte Carlo Methods

This comparison teaches an important lesson about Monte Carlo techniques: Monte Carlo techniques are generally the slowest way possible to solve any given problem, and excel only when there is no other way of solving the problem. It just so happens that there are many problems for which Monte Carlo techniques are the only means of solution.

Multidimensional integration is a good example of such a problem. We can use a Monte Carlo sampling technique to evaluate an integral:

Monte Carlo Integration

In one dimension this isn’t much use; there are much faster ways of evaluating the integral. But when we get much above three dimensions, Monte Carlo methods become generally the only way to evaluate the integrals.

**Python Monte Carlo Integration**

There are many more uses of Monte Carlo techniques, from exact quantum mechanical techniques to economic modeling. We’re going to skip these, since this is a programming article, not a math or science article. But Python is wonderful for all sorts of problems, and perhaps we’ve opened up some new possibilities you haven’t considered. That’s always fun!