Unit 11, lesson 72

Monte Carlo Probability and Expectation

How the Monte Carlo method estimates an expectation or a probability by averaging simulated draws, and what its standard error, the square-root cost of precision and an approximate interval mean.

  • Video 72
  • Exact, Approximate and Simulated Answers

What is E[U2]E[U^2] when U∼Uniform⁡(0,1)U \sim \operatorname{Uniform}(0, 1)? An integral gives 13\tfrac13 at once, but a computer can also estimate it by drawing many values of UU, squaring each one and averaging the squares. The same idea works for models where no formula is available, and the estimate becomes useful once we can say how far a simulated average is likely to land from its target. This lesson sets up the Monte Carlo estimator for expectations and probabilities, computes its standard error, and separates simulation error, which more draws reduce, from model error, which they do not.

Key terms

Monte Carlo method

Definition The Monte Carlo method estimates a probability or an expectation under a probability model by simulating draws from the model and averaging a summary of them. A Monte Carlo calculation has four parts, namely the target quantity, the probability model, the rule that generates draws and the rule that summarizes them. Computers produce pseudorandom numbers, a deterministic sequence fixed by an algorithm and a starting value called the seed, and use them in place of independent uniform draws. A run is reproducible when its generator and software version, seed, sample size, order of draws and the steps that turn draws into values are recorded.

In plain words Instead of computing an expectation with a sum or an integral, let the computer play out the random experiment many times and average the results. The law of large numbers is the reason this works. A simulation answers a question about the model. More draws make the estimate of the model’s value more precise, but they cannot show that the model fits the real process. A wrong rate, a wrong support or ignored dependence is a model error, which is separate from simulation error.

Example To estimate θ=E[U2]\theta = E[U^2] for U∼Uniform⁡(0,1)U \sim \operatorname{Uniform}(0, 1), the target is θ\theta, the model is the uniform law, each draw is a pseudorandom number UiU_i, and the summary is the average of the squares Ui2U_i^2. One recorded run in Python 3.12 used the Mersenne Twister generator random.Random with seed 5001250012 and took 10,00010{,}000 draws, one call of random() for each draw. Repeating the same steps with the same generator and seed reproduces its results exactly.

Monte Carlo estimator

Definition Let XX be a random variable from the model and gg a function with E[∣g(X)∣]<∞E[\lvert g(X) \rvert] < \infty, and let θ=E[g(X)]\theta = E[g(X)]. For independent draws X1,…,XnX_1, \dots, X_n, each with the distribution of XX, the Monte Carlo estimator is

θ^n=1n∑i=1ng(Xi).\hat{\theta}_n = \frac{1}{n} \sum_{i=1}^{n} g(X_i).

It is unbiased, E[θ^n]=θE[\hat{\theta}_n] = \theta, and by the weak law of large numbers P(∣θ^n−θ∣≥ε)→0P(\lvert \hat{\theta}_n - \theta \rvert \ge \varepsilon) \to 0 as n→∞n \to \infty for every ε>0\varepsilon > 0. For a fixed nn it is a random variable, and each run produces one value of it.

In plain words Apply the function to each draw first and average afterward. The order matters. Averaging the values Ui2U_i^2 estimates E[U2]=13E[U^2] = \tfrac13, while squaring the average of the UiU_i estimates (E[U])2=14(E[U])^2 = \tfrac14, which is a different number.

Example For U∼Uniform⁡(0,1)U \sim \operatorname{Uniform}(0, 1) and g(u)=u2g(u) = u^2 the target is

θ=E[U2]=∫01u2 du=13.\begin{aligned} \theta = E[U^2] &= \int_0^1 u^2 \, du \\[0.4em] &= \frac13 . \end{aligned}

The recorded run of 10,00010{,}000 draws gives θ^≈0.329583\hat{\theta} \approx 0.329583, so the error of this run is θ^−θ≈−0.003750\hat{\theta} - \theta \approx -0.003750. The same draws give Uˉ2≈0.245527\bar{U}^{2} \approx 0.245527, which is close to 14\tfrac14, not to 13\tfrac13.

Probability as an expected indicator

Definition For an event AA, the indicator IAI_A equals 11 when AA occurs and 00 otherwise, so E[IA]=P(A)E[I_A] = P(A). A probability is therefore an expectation, and its Monte Carlo estimator is the fraction of draws in which the event occurs. For the event that XX falls in a set BB, with IiI_i the indicator of Xi∈BX_i \in B,

p^n=1n∑i=1nIi\hat{p}_n = \frac{1}{n} \sum_{i=1}^{n} I_i

estimates p=P(X∈B)p = P(X \in B).

In plain words Score each draw 11 if the event happens and 00 if it does not, and average the scores. The result is a relative frequency in a finite simulation, an estimate of the model probability rather than the probability itself. A run with no occurrences gives the estimate 00, which does not show that the event has probability 00.

Example Under the uniform model the event U>0.7U > 0.7 has probability P(U>0.7)=0.3P(U > 0.7) = 0.3. Of the same 10,00010{,}000 recorded draws, 29812981 exceed 0.70.7, so p^=2981/10000=0.2981\hat{p} = 2981/10000 = 0.2981. The squares and the indicators come from the same draws, so the two estimates from this run are not independent of each other.

Monte Carlo standard error

Definition Suppose the draws are independent and σ2=Var⁡(g(X))\sigma^2 = \operatorname{Var}\bigl(g(X)\bigr) is finite. Then

Var⁡(θ^n)=σ2n,SE⁡(θ^n)=σn.\begin{gathered} \operatorname{Var}(\hat{\theta}_n) = \frac{\sigma^2}{n}, \\[0.4em] \operatorname{SE}(\hat{\theta}_n) = \frac{\sigma}{\sqrt{n}}. \end{gathered}

The Monte Carlo standard error is this standard deviation of the estimator across repeated independent runs of the same size. For an indicator σ2=p(1−p)\sigma^2 = p(1 - p). When σ\sigma is unknown it is estimated, by the sample standard deviation of the values g(Xi)g(X_i), or by p^(1−p^)\sqrt{\hat{p}(1 - \hat{p})} for an indicator, and the estimated standard error is itself random.

In plain words The standard error tells how much the estimate would vary from run to run. It is not a bound on the error of a particular run, which can exceed one standard error and now and then two. The formula needs independent draws and a finite variance, and it measures sampling variability within the model, not whether the model is right.

Example For g(u)=u2g(u) = u^2 the fourth moment E[U4]=∫01u4 du=15E[U^4] = \int_0^1 u^4 \, du = \tfrac15 gives

Var⁡(U2)=15−19=445.\operatorname{Var}(U^2) = \frac15 - \frac19 = \frac{4}{45}.

With n=10,000n = 10{,}000 the standard error is 4/450000≈0.002981\sqrt{4/450000} \approx 0.002981, and the run’s error of −0.003750-0.003750 is about −1.26-1.26 standard errors, which is compatible with ordinary run-to-run variation. For the event U>0.7U > 0.7 the standard error is 0.3⋅0.7/10000≈0.004583\sqrt{0.3 \cdot 0.7/10000} \approx 0.004583, and the estimate 0.29810.2981 is about 0.410.41 standard errors below 0.30.3. The estimated version 0.2981⋅0.7019/10000≈0.004574\sqrt{0.2981 \cdot 0.7019/10000} \approx 0.004574 is close to the exact one.

Square-root cost of precision

Definition Suppose 0<σ2<∞0 < \sigma^2 < \infty. Because SE⁡(θ^n)=σ/n\operatorname{SE}(\hat{\theta}_n) = \sigma/\sqrt{n}, multiplying the number of independent draws by a factor cc divides the standard error by c\sqrt{c}. Halving the standard error takes four times as many draws, and a tenfold reduction takes one hundred times as many. A target standard error s>0s > 0 is reached exactly when n≥σ2/s2n \ge \sigma^2/s^2. If σ=0\sigma = 0, every value g(Xi)g(X_i) equals θ\theta with probability 11 and the standard error is 00 for every nn.

In plain words Precision gets expensive. Doubling the draws shrinks the standard error only by the factor 2≈1.41\sqrt{2} \approx 1.41. The rule describes the spread of the estimator, so a particular run with four times as many draws need not have half the error of a particular smaller run.

Example For E[U2]E[U^2] the standard error is about 0.0029810.002981 with 10,00010{,}000 draws, 0.0021080.002108 with 20,00020{,}000 draws and 0.0014910.001491 with 40,00040{,}000 draws. For the event U>0.7U > 0.7, a standard error of 0.0010.001 needs

n≥0.3⋅0.70.0012=210,000n \ge \frac{0.3 \cdot 0.7}{0.001^2} = 210{,}000

draws.

Approximate interval for a Monte Carlo estimate

Definition If the draws are independent and 0<σ2<∞0 < \sigma^2 < \infty, the central limit theorem says that (θ^n−θ)/(σ/n)(\hat{\theta}_n - \theta)/(\sigma/\sqrt{n}) converges in distribution to N(0,1)N(0, 1). For large nn this motivates the approximate 9595 percent interval

θ^n±1.96 σn,\hat{\theta}_n \pm 1.96 \, \frac{\sigma}{\sqrt{n}},

with σ\sigma replaced by an estimate when it is unknown. Its coverage over repeated runs is close to 0.950.95 when nn is large, but it is not guaranteed at a fixed nn. With the exact standard error, Chebyshev’s inequality gives a guarantee for every nn and every k>0k > 0,

P(∣θ^n−θ∣≥kSE⁡(θ^n))≤1k2.\begin{aligned} &P\bigl(\lvert \hat{\theta}_n - \theta \rvert \ge k \operatorname{SE}(\hat{\theta}_n)\bigr) \\[0.4em] &\quad \le \frac{1}{k^2}. \end{aligned}

In plain words When the normal approximation is good, about 9595 percent of runs produce an interval that covers the target. A particular interval either contains θ\theta or it does not, and the 9595 percent describes the procedure. The interval accounts for simulation error only, so a precise interval for a wrong model is still wrong about the real process.

Example For the recorded run, 0.329583±1.96⋅0.0029810.329583 \pm 1.96 \cdot 0.002981 runs from about 0.3237390.323739 to 0.3354270.335427 and contains 13\tfrac13. For the event, 0.2981±1.96⋅0.0045830.2981 \pm 1.96 \cdot 0.004583 runs from about 0.2891180.289118 to 0.3070820.307082 and contains 0.30.3. Chebyshev’s inequality with k=20≈4.47k = \sqrt{20} \approx 4.47 bounds the miss probability by 120=0.05\tfrac{1}{20} = 0.05, and for the squared-uniform estimate it gives the wider interval 0.329583±1750.329583 \pm \tfrac{1}{75}, from about 0.3162500.316250 to 0.3429160.342916.

Check yourself

  1. A study estimates P(U<0.2)P(U < 0.2) for U∼Uniform⁡(0,1)U \sim \operatorname{Uniform}(0, 1) from nn independent draws. Find the standard error of the estimate for n=2500n = 2500, and the number of draws that gives a standard error of 0.0040.004.
Show the answer

The indicator of U<0.2U < 0.2 has variance 0.2⋅0.8=0.160.2 \cdot 0.8 = 0.16, so for n=2500n = 2500 the standard error is

0.162500=0.450=0.008.\sqrt{\frac{0.16}{2500}} = \frac{0.4}{50} = 0.008.

A standard error of 0.0040.004 is half of that, so it takes four times as many draws, n=10,000n = 10{,}000. The condition n≥0.16/0.0042=10,000n \ge 0.16/0.004^2 = 10{,}000 gives the same answer.

  1. A Monte Carlo study estimates E[U3]E[U^3] for U∼Uniform⁡(0,1)U \sim \operatorname{Uniform}(0, 1) with independent draws. Find the target, the variance of one value U3U^3, and the standard error for 90009000 draws. How many draws halve that standard error?
Show the answer

The target is E[U3]=∫01u3 du=14E[U^3] = \int_0^1 u^3 \, du = \tfrac14. The second moment of U3U^3 is the sixth moment of UU, E[U6]=17E[U^6] = \tfrac17, so

Var⁡(U3)=17−116=9112.\operatorname{Var}(U^3) = \frac17 - \frac{1}{16} = \frac{9}{112}.

For 90009000 draws the estimator has variance 9/(112⋅9000)=1/1120009/(112 \cdot 9000) = 1/112000, so the standard error is 1/112000≈0.002988\sqrt{1/112000} \approx 0.002988. Halving it takes four times as many draws, 36,00036{,}000, which give a standard error of about 0.0014940.001494.

  1. A simulation of 10001000 independent draws sees no occurrence of an event, and a report states that the event has probability 00 with standard error 0⋅1/1000=0\sqrt{0 \cdot 1/1000} = 0. What is wrong with this report?
Show the answer

The estimate 0/10000/1000 is a relative frequency, not the model probability, and the plug-in standard error is 00 only because it replaces the unknown pp by 00. Small positive probabilities easily produce a run with no occurrences. For p=0.001p = 0.001 the chance of seeing none in 10001000 independent draws is 0.9991000≈0.3680.999^{1000} \approx 0.368, and every pp up to about 0.0030.003 gives that outcome a chance of at least 55 percent. The report should state the count, 00 occurrences in 10001000 draws, and should not claim that the probability is 00 or that the estimate is exact.

How this lesson connects