A Simple Monte Carlo Program
Let’s start with one of the simplest Monte Carlo (MC) programs. MC programs give a statistical estimate of an answer, and this estimate gets more and more accurate the longer you run it. This basic characteristic of simple programs producing noisy but ever- better answers is what MC is all about, and it is especially good for applications like graphics where great accuracy is not needed.
Estimating Pi
As an example, let’s estimate π. There are many ways to do this, with the Buffon Needle problem being a classic case study. We’ll do a variation inspired by that. Suppose you have a circle inscribed inside a square:
Figure 1: Estimating π with a circle inside a square
Now, suppose you pick random points inside the square. The fraction of those random points that end up inside the circle should be proportional to the area of the circle. The exact fraction should in fact be the ratio of the circle area to the square area. Fraction:
πr2(2r)2=π4
Since the rr cancels out, we can pick whatever is computationally convenient. Let’s go with r=1r=1, centered at the origin:
#include “rtweekend.h”
#include <iostream>
#include <iomanip>
#include <math.h>
#include <stdlib.h>
int main() {
int N = 1000;
int inside_circle = 0;
for (int i = 0; i < N; i++) {
auto x = random_double(-1,1);
auto y = random_double(-1,1);
if (x*x + y*y < 1)
inside_circle++;
}
std::cout << std::fixed << std::setprecision(12);
std::cout << “Estimate of Pi = ” << 4*double(inside_circle) / N << ‘\n’;
}
Listing 1: [pi.cc] Estimating π
The answer of ππ found will vary from computer to computer based on the initial random seed. On my computer, this gives me the answer Estimate of Pi = 3.0880000000