Python for simulations

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². 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³, 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) } plot(f,0,1,cex=2) xvals = seq(0,1,length=100) yvals = f(xvals) x = c(xvals,rev(xvals)) y = c(rep(0,100),rev(yvals)) polygon(x,y,col='gray') text(0.2,0.7,expression(y == x^2),cex=2) lines(c(0,1),c(0,1),lty=2)