Monte Carlo Estimation of Pi and e

Worked Solution

How to Think About It: Monte Carlo is all about finding a random experiment whose expected value equals the constant you want, then averaging many independent copies. For $\pi$, the classic trick is that the area of a quarter circle inscribed in the unit square is $\pi/4$ -- throw darts at the unit square and the fraction landing inside the circle converges to $\pi/4$. For $e$, you need a different experiment: it turns out the expected number of independent $U(0,1)$ draws needed for the running sum to first exceed $1$ is exactly $e$.

Quick Estimate: With $N = 10{,}000$ samples, the $\pi$ estimator $\hat\pi = 4H/N$ has standard error $4\sqrt{\tfrac{(\pi/4)(1-\pi/4)}{N}} = \sqrt{\tfrac{\pi(4-\pi)}{N}} \approx \sqrt{2.70/10000} \approx 0.016$, so about $2$ decimal places. For $e$, the stopping count $N$ has variance $\operatorname{Var}(N) = 3e - e^2 \approx 0.765$ (derived below), so with $M = 10{,}000$ trials the standard error is $\sqrt{(3e-e^2)/M} \approx \sqrt{0.765/10000} \approx 0.0088$.

Approach: For each constant, identify an unbiased random observable (an indicator for $\pi$, a stopping count for $e$), give the sampling algorithm, justify it via the law of large numbers, and read off the $O(1/\sqrt{N})$ Monte Carlo convergence rate from the CLT.

Formal Solution:

Part (a) -- Estimating $\pi$. Consider the unit square $[0,1]^2$. The quarter disk $\{(x,y): x^2 + y^2 \le 1\}$ has area $\pi/4$.

Algorithm: 1. Sample $N$ points $(x_i, y_i)$ uniformly from $[0,1]^2$. 2. Count hits $H = \sum_{i=1}^{N} \mathbf{1}[x_i^2 + y_i^2 \le 1]$. 3. Estimate $\hat{\pi} = 4H/N$.

Why it works: $E\!\left[\mathbf{1}[x^2+y^2\le 1]\right] = \pi/4$, so by the law of large numbers $H/N \to \pi/4$ and $\hat\pi \to \pi$. Since $H \sim \text{Bin}(N, \pi/4)$, $\operatorname{Var}(\hat\pi) = 16\cdot\tfrac{(\pi/4)(1-\pi/4)}{N} = \tfrac{\pi(4-\pi)}{N}$.

Part (b) -- Estimating $e$. Let $U_1, U_2, \dots \sim U(0,1)$ be i.i.d. and define the stopping count $$N = \min\Bigl\{k : U_1 + U_2 + \cdots + U_k > 1\Bigr\}.$$

Algorithm: 1. Repeat $M$ independent trials. 2. In each trial, draw $U(0,1)$ values and accumulate the sum until it first exceeds $1$; record the count $N_j$. 3. Estimate $\hat{e} = \frac{1}{M}\sum_{j=1}^{M} N_j$.

Why it works: the event $\{N > n\}$ means the first $n$ uniforms sum to at most $1$, i.e. $(U_1,\dots,U_n)$ lies in the $n$-simplex $\{u_i \ge 0,\ \sum u_i \le 1\}$, whose volume is $1/n!$. So $P(N > n) = 1/n!$ for $n \ge 0$, and by the tail-sum formula for a nonnegative integer random variable, $$E[N] = \sum_{n=0}^{\infty} P(N > n) = \sum_{n=0}^{\infty} \frac{1}{n!} = e.$$

*Variance of the stopping count.* Use $E[N(N-1)] = \sum_{n\ge 0} 2n\,P(N>n)$ (or the pmf $P(N=n) = \tfrac{1}{(n-1)!} - \tfrac{1}{n!}$ for $n\ge 1$). The second factorial moment is $$E[N(N-1)] = \sum_{n=0}^{\infty} 2n\cdot\frac{1}{n!} = 2\sum_{n=1}^{\infty}\frac{1}{(n-1)!} = 2e,$$ so $E[N^2] = E[N(N-1)] + E[N] = 2e + e = 3e$. Therefore $$\operatorname{Var}(N) = E[N^2] - (E[N])^2 = 3e - e^2 \approx 3(2.71828) - (2.71828)^2 \approx 0.7658.$$ Hence the Monte Carlo standard error of $\hat e$ over $M$ trials is $\sqrt{(3e - e^2)/M}$.

Code: ```python import random

def estimate_pi(N=1_000_000): hits = sum(1 for _ in range(N) if random.random()2 + random.random()2 <= 1) return 4.0 * hits / N

def estimate_e(M=1_000_000): total = 0 for _ in range(M): s, k = 0.0, 0 while s <= 1.0: s += random.random() k += 1 total += k return total / M ```

Convergence: Both estimators are sample averages of i.i.d. observables with finite variance, so by the CLT they converge at rate $O(1/\sqrt{N})$ -- the standard Monte Carlo rate. To gain one extra decimal place of accuracy you need about $100\times$ more samples. Variance-reduction techniques (antithetic variates, importance sampling, control variates) improve the constant but not the $O(1/\sqrt{N})$ rate.

Answer: For $\pi$: sample uniform points in the unit square and estimate $\hat\pi = 4\times(\text{fraction inside the quarter disk})$, since that fraction has mean $\pi/4$ (estimator variance $\pi(4-\pi)/N$). For $e$: count how many $U(0,1)$ draws are needed for the running sum to first exceed $1$, and average over many trials, because $P(N>n) = 1/n!$ gives $E[N] = \sum 1/n! = e$; the stopping count has $\operatorname{Var}(N) = 3e - e^2 \approx 0.765$, so the standard error over $M$ trials is $\sqrt{(3e-e^2)/M}$. Both estimators converge at the Monte Carlo rate $O(1/\sqrt{N})$.