Error in calculation of $\pi$ using Monte Carlo method.

Question

So lets say we are trying to calculate value of $\pi$ using MonteCarlo method. By picking random points in a square and measuring their distance from the center and if $k$ points lie inside the circle using the ratio $\frac{k}N$ to calculate the value of $\pi$.

How do I derive the formula for variance of the calculated value of $\pi$?

This is what I have got so far.

1. Say the probability that the point lies inside the circle is $p$.

2. Say we have N random points in the square, of which k lie inside the circle.

3. The probability that k points lie inside the circle is $\binom{N}{k}p^k(1-p)^{N-k}$.

4. The error in estimate of $\pi$ is $e=||c\pi-\frac{N}k||$.

5. So the expected value of error $E(e) = \sum_{k=0}^{k=N}\left\|c\pi-\frac{N}k\right\|\binom{N}{k}p^k(1-p)^{N-k}$?

Is the expectation value of error that I have written correct, how do I evaluate it?

Answer 1

For very large $n$, as in a simulation, a reasonable 95% confidence interval for the estimate $\hat p = k/n$ of the proportion of the points inside the circle is of the form

$$\hat p \pm 1.96\sqrt{\frac{\hat p(1-\hat p)}{n}}.$$

So the approximate margin of simulation error is $1.96\sqrt{\frac{\hat p(1-\hat p)}{n}}.$

We have $Var(\hat p) = Var(k/n) = (1/n^2)Var(k) = p(1-p)/n,$ so $SD(\hat p) = \sqrt{\frac{p(1-p)}{n}},$ where $p$ is the unknown theoretical fraction. We use the observed value $\hat p$ to get $\widehat{SD}(\hat p),$ the estimated standard deviation of $\hat p.$

---

But be careful! You are not specific about 'a square' and 'the circle', but the usual simulation is to use a square of area $4$ with vertices at $(-1,-1)$ and $(1,1)$ and the inscribed circle of radius $1$ and area $\pi.$ Thus the ratio of 'accepted' points in the circle to the total points in the square estimates $\pi/4.$ If you multiply the estimate by $4$ to get an estimate of $\pi,$ then the margin of error is also multiplied by $4.$

Let's try the actual simulation experiment in R statistical software:

> set.seed(1066)  # uae this seed to replicate this simulation
> n = 10^6;  x = runif(n, -1, 1);  y = runif(n, -1, 1)
> r = sqrt(x^2 + y^2)
> p = sum(r < 1)/n;  pi.est = 4*p;  pi.est
[1] 3.140892
> abs(pi.est - pi)  # observed simulation error
[1] 0.0007006536
> me.p = 2*sqrt(p*(1-p)/n); me.p
[1] 0.0008213351
> me.pi = 4*me.p;  me.pi  # 95% est. marg of sim err, using sim results
[1] 0.00328534
q = pi/4
4*1.96*sqrt(q*(1-q)/m)    # exact 95% marg of sim err, using pi=3.14.59...
[1] 0.003218679

The observed simulation error (about 0.001) for estimating $\pi$ is within the 95% simulation margin of error (about 0.003) for the estimate of $\pi$, but not within the 95% margin of simulation error (about 0.0001) for estimating $\pi/4.$

Here is a plot of 10,000 of the million points.