Approximating π using Monte Carlo properly

Approximating \(\pi\) using Monte Carlo

A Monte Carlo algorithm uses random variables to approximate some solution to which the solution might be unfeasable to calculate directly. The example of choice of every professor I know of is to approximate \( \pi = 3.14\dots \).

The procedure is at first glance simple enough:

1. Generate \( n \) samples uniformly from the unit square,\(X_i= (x_i, y_i) \sim Un([0, 1]^2) \)
2. If \( x_i ^2 + y_i^2 \leq 1 \) then the point in inside the unit circle. Count them as \( n_+\)
3. Since the area of the unit circle is known to be \(\pi\) and we generate points over only a forth of it. Therefore \[ \pi = 4\, p((x_i, y_i)\in B) \approx 4\, \frac{n_+}{n} \]

But when it comes to the error estimate the method breaks down. Some suggest to run the simulation again and again, but that's more than wasteful. That is because the method is using statistics without really understanding it.

... properly

Actually we throw away much information if we just take the mean to approximate \(p\). The decision if a point that is randomly generated is in the circle is actually distributed as Bernoulli (1 with probability \(p\) and 0 with probability \(1-p\) ).

Now we can estimate \(p\) smarter if we recognize it to be Beta distributed -- or its estimate to be Beta distributed if you want to be frequentist. The Bayes estimator (assuming a uniform prior) is \[ \hat p_\sim Beta(1+n_+, 1 + n - n_+)\] with the expected value for \(p\) being \(\mathbb{E}(\hat p) = \frac{1+n_+}{1+n}\) and thus the same as the naive approach except for some smoothing. The variance is \[ var(\hat p) = \frac{(1+n_+)\, (1+n-n_+))}{(2+n)^2\, (n+3)} \] but since you have the full posterior probability distribution, you can be much more precise if you'd like.