Python Monte Carlo Simulation Example
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.
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:
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:
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)
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.
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!