I'd like to show an example of how Python is useful for simulations. I saw it in a terrific book called Probability by Grinstead & Snell. You can get a free pdf of the book online, but after reading that for a while I went out and bought my own copy. They illustrate use of the Monte Carlo method for integrating a function. I always thought it was Johnny von Neumann who coined the name, but the Wikipedia article attributes it to Stanley Ulam and his gambling uncle who was partial to a casino at Monte Carlo.

In any event, it works just like Battleship (png).

We throw darts at a section of the 2D plane, let's say a square containing points 0 < x <= 1 and 0 < y <= 1, by getting two random numbers, x and y. We then evaluate the function f(x) and if y <= f(x), we call that a "hit." We calculate the fraction of darts that hit the target.

Below is a Python program that integrates y = x

^{2}. We can see at a glance that the shaded area is somewhat less than 1/2, but how much less? If you remember your calculus, you'll know our function integrates to 1/3 x

^{3}, evaluated between 0 and 1, which equals 1/3. Python gives us 3 significant digits from a million trials. Pretty slick. It can also be quite inefficient. I should be able to use my recent study of elementary statistics to calculate the expected accuracy for a given number of trials, but it eludes me right now. I do know that you'll use up a lot of darts trying to get an approximation for pi by using a circle as the target.

Here is the R code for the plot:

f <- function(x) { return (x*x) } |