![](https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhfx6fjo3MB1GaJojbzW-xWUZbts1iRlG8XE_7R29LZroq3a5_2FHAOUpL1eNtEoUlwIsQfgwj-se4_WVISsoeT81HQ-h0o7E6YmI7U4M9m1bEFkZaGi797-nk9BYLMoKCpJq_a6y7waevy/s320/Monte_carlo_method.png)
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.
![](https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhz1BNN8tX-s7woYJjFwKjsaaSihs_mlLjz1YHlc_EF6rKLDgCZ2hXv_4YD57vuCDYok0KIftXoSO7CRu-ogiajuJk9kFxfWVyIaVcIWmHYqpk6YbGzMwbTmz-eiDtWzRjGs3NiCbAibKTR/s320/y=f(x).png)
Below is a Python program that integrates y = x2. 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 x3, 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.
![](https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiqXELxkyTGWr-AI4tZBBe0E_sCPcEp0zu5YPBXoWehtMrz9G-TZNR3G6XMPeXDRNrXs_g-U4FLA0lEGcrnFag9Cb6QtIH2uqLL3luDH-hHm-lMakCJa3hEiABluJ56iwWWdLVIpXJG_vpE/s320/monte+carlo.png)
Here is the R code for the plot:
f <- function(x) { return (x*x) } |