Are all Monte Carlo algorithms to approximate $\pi$ equivalent?

935 Views Asked by At

There are several ways one can approximate $\pi$ using a Monte-Carlo type algorithm. For example, one can draw random points in the unit square, and approximate $\pi$ via the ratio of points that fall in the unit disk. Another way would be to use an algorithm based on Buffon's needle problem: from the proportion of needles that fall across the parallel lines, one can approximate $\pi$. One can presumably devise many other random algorithms to estimate $\pi$.

My question is: in terms of speed of convergence, are all these algorithms equivalent? Or are there some that will yield correct decimals of $\pi$ faster? What is the Monte Carlo method to estimate $\pi$ that has the fastest convergence?

2

There are 2 best solutions below

0
On

Trying to keep this simple, answering your original very good question.(+1) Using Riemann approximation by rectangles to get the area under the curve $y = \sqrt{1-x^2},$ for $0 < x < 1,$ and multiplying by 4, you can get pretty close to $\pi$ with only a dozen rectangles. Closer if you use trapezoids. This is not a Monte Carlo method, but the first method below is somewhat analogous (using random points $x,$ instead of a deterministic grid of center points as bases of rectangles).

set.seed(430)  # retain for exactly the same simulation
m = 10^6;  x = runif(m);  h = sqrt(1-x^2)
mean(h)*4
[1] 3.142033           # estimate of pi
2*4*sd(h)/sqrt(m)
[1] 0.001784625        # rough 95% margin of simulation error

So a million random $x$-values gives $3.1420 \pm 0.0018.$

An acceptance-rejection method similar to the one you suggested is somewhat less efficient.

set.seed(430)
m = 10^6;  x = runif(m);  y=runif(m)
acc = y <+ sqrt(1-x^2)
mean(acc)*4
[1] 3.143212           # aprx pi
2*4*sd(acc)/sqrt(m)
[1] 0.003282115        # marg of sim err

Here we get $3.1433 \pm 0.0033,$ which indicates somewhat slower convergence.

Notes: (a) My purpose here is not to give best simulation methods for approximating $\pi,$ but to answer your question by showing two commonly used methods with different rates of convergence. (b) I do not challenge the importance of the Metropolis-Histings algorithm in simulation, but in my experience it is more useful in dimensions higher than one or two. (c) I am not from the 'show-me' state of Missouri, but I have a sister who lives there; if there is a way to get $\pi$ to 47 decimal places by simulating a couple of coin flips, that's news to me and I'd want to see an authoritative reference for that. (d) I have illustrated other Monte Carlo methods on another page; some of them can be adapted to your question.

4
On

All random algorithms that approximate $\pi$ by trying something that succeeds with probability $f(\pi)$ lots of times, getting an estimate of $f(\pi)$, and using that to estimate $\pi$ will be approximately equivalent. (Whether the thing we try is "drop a needle that interects a line" or "pick a point that lands inside a circle" or anything else.)

The reason is that their convergence rate will have nothing to do with the method at all. If the probability we're trying to estimate is $p$, then after $n$ trials the number of successes will follow a binomial distribution, and the typical number of successes will be $np \pm C \sqrt{np(1-p)}$. This gives an estimate of $p$ to within $O(\frac1{\sqrt n})$ error, which means millions of trials before we can estimate $p$ to three decimal places.

I only say "approximately" equivalent because the probability we're estimating can depend on $\pi$ in lots of ways. If the probability of success is $\frac12 + \frac{\pi}{1000}$, then it will take many trials before our estimate of this probability becomes significantly better than "it's about $\frac12$". On the other hand, if the probability of success is $10000\pi - 31415$, then estimating that the probability is approximately $0.9$ is enough to give use $\pi$ to five digits after the decimal: $\pi \approx 3.14159$.