Worked examples — Monte Carlo simulation — law of large numbers basis

Before symbols fly, one reminder of the three characters from the parent note:

---

The scenario matrix

Every Monte Carlo problem falls into one of these case classes. The table is the map; the examples below fill every row.

Cell	Case class	What makes it different	Example
A	Bounded integrand on $[0,1]$	Textbook case, finite variance	Ex 1
B	Integral on $[a,b]\ne [0,1]$	Must rescale to a uniform	Ex 2
C	Improper / infinite range $[0,\infty )$	Uniform impossible — change the sampling density	Ex 3
D	Probability = indicator average (0/1 output)	$g$ returns only 0 or 1 (Bernoulli)	Ex 4
E	Degenerate: $g$ constant → zero variance	$\text{SE}=0$, converges in one sample	Ex 5
F	Limiting behaviour: how does error scale?	Plan $n$ for a target accuracy	Ex 6
G	Heavy tail: variance infinite → LLN breaks	The one case that does NOT converge	Ex 7
H	Real-world word problem	Translate story → expectation	Ex 8
I	Exam twist: variance reduction	Same $\theta$, smaller $\text{SE}$	Ex 9

Prerequisite links you may want open: Bernoulli Distribution, Variance and Covariance, Central Limit Theorem, Chebyshev's Inequality, Numerical Integration, Importance Sampling, Strong Law of Large Numbers. This page is a child of Monte Carlo simulation — law of large numbers basis.

---

Ex 1 — Cell A: bounded integrand on $[0,1]$

Active symbol: the random input is $U\sim \text{Unif}(0,1)$ (uniform).

Forecast: guess the true value and guess whether $n=1000$ samples get you 2 decimals right. (Write it down before reading on.)

Step 1 — Turn the integral into an expectation. ${\int }_0^1{x}^{2}dx=\mathbb{E}[U^2],U\sim \text{Unif}(0,1).$ Why this step? The uniform density on $[0,1]$ is the constant $1$, so ${\int }_0^{1}g(x)\cdot 1dx$ is by definition the average of $g(U)$. Here $g(x)=x^2$.

Step 2 — Write the estimator. ${\hat{\theta }}_n=\frac{1}{n}{\sum }_{i=1}^n{U}_i^2.$ Why? Replace the true average by the sample average — the LLN promises it converges.

Step 3 — True value & true variance (so we know what to expect). $\theta =[\frac{x^3}{3}{]}_0^1=\frac{1}{3},\mathbb{E}[U^4]=\frac{1}{5},{\sigma }_g^2=\mathbb{E}[U^4]-{\theta }^2=\frac{1}{5}-\frac{1}{9}=\frac{4}{45}.$ Why compute ${\sigma }_g^2$? Because $\text{SE}={\sigma }_g/\sqrt{n}$ tells us the error size before running anything.

Step 4 — Predicted error at $n=1000$. $\text{SE}=\sqrt{\frac{4}{45}}/\sqrt{1000}=\frac{0.2981}{31.62}\approx \mathrm{0.00943.}$ Why plug in $n=1000$? We picked the sample size the forecast asked about, so we can decide in advance whether it buys 2 decimals — no simulation needed.

Verify: $\theta =0.3333\dots$, and $\text{SE}\approx 0.0094<0.01$, so a typical run lands within about $\pm 0.019$ (two SEs) of $\frac{1}{3}$ — good to roughly 2 decimals. Forecast checks out. Units: dimensionless (pure number).

---

Ex 2 — Cell B: integral on a general interval $[a,b]$

Active symbol: the random input is $U\sim \text{Unif}(2,5)$ (uniform, but now on $[2,5]$).

Forecast: the interval is length 3, not 1. Guess: does the "$\times$ length" factor go inside or outside the average?

Step 1 — Sample uniformly on $[2,5]$. A point $U\sim \text{Unif}(2,5)$ has density $\frac{1}{b-a}=\frac{1}{3}$ (constant, but now $\frac{1}{3}$ not $1$). Why this step? On a longer interval, "evenly spread" means each unit of length gets probability $\frac{1}{b-a}$, so the density shrinks.

Step 2 — Recover the integral from the expectation. $\mathbb{E}[\sqrt{U}]={\int }_2^5\sqrt{x}\cdot \frac{1}{3}dx=\frac{1}{3}\theta \Rightarrow \theta =(b-a)\mathbb{E}[\sqrt{U}]=3\mathbb{E}[\sqrt{U}].$ Why? The density $\frac{1}{3}$ divides the integral; to undo it we multiply the average by $(b-a)=3$. The length factor sits OUTSIDE the average.

Step 3 — Estimator. ${\hat{\theta }}_n=(b-a)\cdot \frac{1}{n}{\sum }_{i=1}^n\sqrt{U_i}=3\cdot \frac{1}{n}\sum \sqrt{U_i}.$ Why this exact form? Step 2 proved $\theta =3\mathbb{E}[\sqrt{U}]$; we simply replace the true average $\mathbb{E}[\sqrt{U}]$ by the sample average $\frac{1}{n}\sum \sqrt{U_i}$ (LLN) and keep the constant $3$ outside — that is the estimator we will actually run.

Step 4 — True value for checking. ${\int }_2^5{x}^{1/2}dx=[\frac{2}{3}{x}^{3/2}{]}_2^5=\frac{2}{3}(5^{1.5}-2^{1.5})=\frac{2}{3}(11.1803-2.8284)=\mathrm{5.5679.}$ Why compute the exact integral here? This case has an elementary antiderivative $\frac{2}{3}{x}^{3/2}$, so we can get the true $\theta$ by hand and use it as the yardstick the Monte Carlo estimate must match — validation, not the method itself.

Verify: $\theta \approx 5.5679$. Sanity: the interval is width 3 and $\sqrt{x}$ ranges from $\sqrt{2}\approx 1.41$ to $\sqrt{5}\approx 2.24$, so the integral must lie between $3\times 1.41=4.24$ and $3\times 2.24=6.71$. Our $5.57$ sits comfortably inside. ✔

---

Ex 3 — Cell C: improper integral on $[0,\infty )$

Active symbol: the random input is $X\sim \text{Exp}(1)$ (exponential), built from a uniform $U$ via $X=-\ln U$.

Forecast: you cannot sample "uniformly on $[0,\infty )$" — an infinite line has no even spread. Guess how we dodge this.

Figure s01 — what you're looking at: the horizontal axis is $x$ running from $0$ toward infinity; the vertical axis is probability density. The red curve is the sampling density $p(x)=e^{-x}$. The black tick marks along the bottom are 40 actual samples drawn from it — notice they crowd near $x=0$ (where the red curve is tall) and thin out to the right (where it flattens). This is the picture behind Step 1: we sample where the density says the mass is.

Step 1 — Recognise a built-in density. The function $p(x)=e^{-x}$ for $x\ge 0$ is a genuine probability density: it is positive and ${\int }_0^{\infty }{e}^{-x}dx=1$. This is the Exponential(1) distribution. Why this step? Since we cannot make a uniform on an infinite line, we instead sample from a density that already lives there. The red curve in the figure packs samples near 0 and thins out toward infinity, exactly where $e^{-x}$ says the "mass" is.

Step 2 — Split the integrand into density $\times$ leftover. $\theta ={\int }_0^{\infty }\underset{\text{density }p(x)}{\underbrace{e^{-x}}}\underset{g(x)}{\underbrace{\cos x}}dx={\mathbb{E}}_{X\sim \text{Exp}(1)}[\cos X].$ Why? Whatever we multiply the density by becomes the function we average. Here the leftover is simply $\cos x$, and the active input is now $X\sim \text{Exp}(1)$ (not a uniform).

Step 3 — Generate exponential samples. Draw $U\sim \text{Unif}(0,1)$ and set $X=-\ln U$. This standard trick (inverse-CDF) turns a uniform into an Exponential(1). Why? We only have uniform RNGs; $-\ln U$ reshapes them into the density we need. This is the seed of Importance Sampling.

Step 4 — Estimator. ${\hat{\theta }}_n=\frac{1}{n}{\sum }_{i=1}^n\cos (X_i),X_i=-\ln U_i.$ Why this form? Step 2 rewrote $\theta$ as $\mathbb{E}[\cos X]$ with $X\sim \text{Exp}(1)$, so the sample mean of $\cos (X_i)$ is the natural estimator — no extra density factor is needed because the $e^{-x}$ is absorbed into how we sample $X$.

Step 5 — True value. ${\int }_0^{\infty }{e}^{-x}\cos xdx=\frac{1}{1^2+1^2}=\frac{1}{2}.$ Why this formula? Standard result ${\int }_0^{\infty }{e}^{-ax}\cos (bx)dx=\frac{a}{a^2+b^2}$ with $a=b=1$.

Verify: $\theta =0.5$. Sanity: $\cos X$ oscillates between $-1$ and $1$ but the samples cluster near small $x$ where $\cos x\approx 1$, so a positive value below 1 is expected. $0.5$ ✔

---

Ex 4 — Cell D: probability as an indicator average

Active symbol: the random input is the pair $(X,Y)$, each coordinate an independent $\text{Unif}(0,1)$.

Forecast: the region above $y=x^2$ inside the unit square — guess its area as a fraction.

Figure s02 — what you're looking at: the axes are the dart's coordinates $X$ (horizontal) and $Y$ (vertical), both in $[0,1]$. The black curve is the parabola $y=x^2$. The red shaded band is the region where $Y>X^2$ — the "success" region. Dots are 120 random darts: red dots landed above the curve (score $1$), black dots below (score $0$). The estimate is just the fraction that are red, which visually is the red area $=2/3$.

Step 1 — Define the indicator. $g(X,Y)=1\{Y>X^2\}=\{\begin{array}{ll}1& \text{point is above the curve}\\ 0& \text{otherwise}\end{array}$ Why this step? A probability is the expectation of a 0/1 indicator: $P(A)=\mathbb{E}[1_A]$. The output of $g$ is a coin flip — this is the Bernoulli Distribution.

Step 2 — Estimator = fraction of "yes" darts. ${\hat{\theta }}_n=\frac{1}{n}{\sum }_{i=1}^{n}1\{Y_i>X_i^2\}.$ Why? Averaging 0s and 1s just counts the fraction of successes — exactly the red-dot fraction in the figure.

Step 3 — True value. The area above the parabola is total square minus area under it: $\theta =1-{\int }_0^1{x}^{2}dx=1-\frac{1}{3}=\frac{2}{3}.$ Why subtract? The whole square has area $1$; the region under the parabola has area ${\int }_0^1{x}^{2}dx=\frac{1}{3}$; everything else is the "above" region, so $1-\frac{1}{3}$.

Step 4 — Its variance (Bernoulli). With $p=\frac{2}{3}$: ${\sigma }_g^2=p(1-p)=\frac{2}{3}\cdot \frac{1}{3}=\frac{2}{9}\approx \mathrm{0.2222.}$ Why care? Indicator outputs always have variance $p(1-p)$, maximised at $p=\frac{1}{2}$. Extreme probabilities (near 0 or 1) are cheaper to estimate.

Verify: $\theta =0.6667$. Sanity: the parabola dips low, so most of the square is above it — a fraction above $\frac{1}{2}$ is expected. $\frac{2}{3}$ ✔

---

Ex 5 — Cell E: degenerate, zero-variance case

Active symbol: the random input is $U\sim \text{Unif}(0,1)$ — though as we'll see, it doesn't matter what we sample.

Forecast: how many samples does this "estimate" need to be exactly right?

Step 1 — Compute the expectation. $\mathbb{E}[7]=7.$ Why? The average of a thing that's always 7 is 7 — no randomness involved.

Step 2 — Compute the variance. ${\sigma }_g^2=\mathbb{E}[(7-7{)}^2]=0.$ Why this matters? $\text{SE}={\sigma }_g/\sqrt{n}=0/\sqrt{n}=0$ for every $n$, even $n=1$.

Step 3 — Interpret the degenerate limit. A single sample gives ${\hat{\theta }}_1=7$ exactly. Adding more samples changes nothing. Why? With zero spread there is no noise to average away — the estimator has already hit the target on sample one.

Verify: $\theta =7$, $\text{SE}=0$. Any run of any length returns exactly 7. ✔

---

Ex 6 — Cell F: limiting behaviour, planning $n$

Active symbol: the random input is $U\sim \text{Unif}(0,1)$, same as Ex 1.

Forecast: Ex 1 hit $\text{SE}\approx 0.0094$ at $n=1000$. To shrink the half-width ~$20\times$, guess: a bit more samples, or a LOT more?

Figure s03 — what you're looking at: both axes are logarithmic. The horizontal axis is the sample size $n$ (from $100$ up to a million); the vertical axis is the $95\%$ half-width $w=1.96{\sigma }_g/\sqrt{n}$. The red line is that half-width — a straight downward line on a log-log plot because $w\propto n^{-1/2}$. The dashed black line marks the target $w=0.001$; the arrow shows the red line crosses it at $n\approx 3.4\times 10^5$. The picture is the $1/\sqrt{n}$ law: to drop $w$ by $10\times$ you slide $100\times$ to the right.

Step 1 — Write the target. The $95\%$ half-width (defined above) is $w=1.96{\sigma }_g/\sqrt{n}$. Demand $1.96\cdot \frac{{\sigma }_g}{\sqrt{n}}\le \mathrm{0.001.}$ Why this step? This half-width is the "$\pm$" of the confidence interval; forcing it below $0.001$ is exactly the accuracy the problem requests, so we solve for the $n$ that achieves it.

Step 2 — Solve for $n$. $\sqrt{n}\ge \frac{1.96{\sigma }_g}{0.001}\Rightarrow n\ge {\left(\frac{1.96{\sigma }_g}{0.001}\right)}^2=\frac{1.96^2\cdot \frac{4}{45}}{0.001^2}.$ Why square? Because error lives under a square root; undoing it squares everything — the source of the brutal $1/\sqrt{n}$ law.

Step 3 — Plug numbers. ${\sigma }_g^2=\frac{4}{45}=0.08889$: $n\ge \frac{3.8416\times 0.08889}{10^{-6}}=\frac{0.34148}{10^{-6}}\approx 341,476.$ So $n\approx 341,477$ samples. Why substitute these particular numbers? $1.96$ is the $95\%$ factor from the Central Limit Theorem, $\frac{4}{45}$ is the variance we already computed in Ex 1, and $10^{-6}=(0.001{)}^2$ is the squared target; plugging them turns the abstract inequality into a concrete sample budget we can actually buy.

Verify: the red line in the figure shows half-width $\propto 1/\sqrt{n}$. Going from half-width $0.0184$ (at $n=1000$, i.e. $1.96\times 0.0094$) down to $0.001$ is an $18.4\times$ reduction, needing $18.4^2\approx 339\times$ more samples: $1000\times 339\approx 3.4\times 10^5$. ✔ Matches. Twenty-fold accuracy costs ~340× the work.

---

Ex 7 — Cell G: heavy tail, LLN breaks

Active symbol: the random input is $X\sim \text{Cauchy}$ (heavy-tailed, not uniform).

Forecast: the density is symmetric about 0, so surely the mean is 0 and the sample mean converges there. Guess whether that reasoning is safe.

Figure s04 — what you're looking at: the horizontal axis is the number of samples $n$ collected so far; the vertical axis is the running sample mean (the average of the first $n$ draws). The black line is a Normal sample: it quickly hugs $0$ and stays there — the LLN working. The red line is a Cauchy sample: it lurches in violent jumps and never settles, because rare gigantic draws from the fat tails keep hijacking the average. This is the visual of "the mean does not exist."

Step 1 — Check the integrability condition. The Weak Law needs $\mathbb{E}|X|<\infty$. Test it: $\mathbb{E}|X|={\int }_{-\infty }^{\infty }\frac{|x|}{\pi (1+x^2)}dx=\infty .$ Why this step? For large $x$ the integrand behaves like $\frac{|x|}{\pi x^2}=\frac{1}{\pi |x|}$, and $\int \frac{dx}{x}$ diverges (logarithmically). The mean does not exist.

Step 2 — Consequence. Because $\mathbb{E}|X|=\infty$, the Strong Law of Large Numbers does not apply. The sample mean ${\bar{X}}_n$ does not settle down — remarkably, ${\bar{X}}_n$ of $n$ Cauchy samples is itself standard Cauchy, no matter how big $n$ is. Why this matters? The red trace in the figure keeps taking wild jumps forever; occasional gigantic samples from the fat tails repeatedly hijack the running average.

Step 3 — The moral. Symmetry suggested "mean $=0$", but the mean is undefined, so there is nothing for the LLN to converge to.

Verify: ${\int }_{-\infty }^{\infty }\frac{|x|}{\pi (1+x^2)}dx$ diverges — we confirm the tail $\frac{1}{\pi x}$ has divergent integral by comparing partial integrals growing without bound. ✔ (See VERIFY.)

---

Ex 8 — Cell H: real-world word problem

Active symbol: the random input per box is a triple of uniforms $U_1,U_2,U_3\sim \text{Unif}(0,1)$; the count of failures $F$ is then a derived random variable.

Forecast: with only a $10\%$ per-bulb failure, guess whether returns are rare ($3\%$) or common ($20\%$).

Step 1 — Model each box. Number of failures $F\sim \text{Binomial}(3,0.1)$ (sum of 3 Bernoulli Distribution flips). Define $g=1\{F\ge 2\}.$ Why this step? "Return probability" is $P(F\ge 2)=\mathbb{E}[1\{F\ge 2\}]$ — again a probability written as an indicator average.

Step 2 — Monte Carlo recipe. Simulate 3 uniforms per box, count how many are $<0.1$ (failures), record whether $\ge 2$, average over many boxes. Why? Each simulated box is one draw of $g$; the LLN converges the fraction returned to the true probability. A uniform below $0.1$ happens exactly $10\%$ of the time, so it faithfully mimics a bulb failing.

Step 3 — Exact value to check against. $P(F\ge 2)=P(F=2)+P(F=3).$ $P(F=2)=\left(\genfrac{}{}{0ex}{}{3}{2}\right)(0.1{)}^2(0.9{)}^1=3\cdot 0.01\cdot 0.9=0.027,$ $P(F=3)=\left(\genfrac{}{}{0ex}{}{3}{3}\right)(0.1{)}^3=\mathrm{0.001.}$ $\theta =0.027+0.001=\mathrm{0.028.}$ Why compute exactly? A word problem with a small discrete model can be checked exactly with the binomial formula $\left(\genfrac{}{}{0ex}{}{n}{k}\right)p^k(1-p{)}^{n-k}$, so we validate the simulator against ground truth instead of trusting it blindly.

Verify: $\theta =0.028$ (about $2.8\%$) — returns are rare, the "rare" forecast wins. Units: a probability, dimensionless, in $[0,1]$. ✔

---

Ex 9 — Cell I: exam twist, variance reduction

Active symbol: the random input is $U\sim \text{Unif}(0,1)$; the antithetic partner $1-U$ is also $\text{Unif}(0,1)$ but built from the same draw.

Forecast: both give the same $\theta$. Guess whether pairing helps, hurts, or does nothing.

Step 1 — True value (for checking). $\theta ={\int }_0^1{e}^{x}dx=e-1\approx \mathrm{1.71828.}$ Why start here? $e^x$ has the elementary antiderivative $e^x$, so we can pin the exact $\theta$ and later confirm both estimators aim at the same target.

Step 2 — Plain estimator variance. $g(x)=e^x$, and ${\sigma }_g^2=\mathbb{E}[e^{2U}]-{\theta }^2=\frac{e^2-1}{2}-(e-1{)}^2.$ Numerically $\frac{e^2-1}{2}=3.19453$ and $(e-1{)}^2=2.95249$, so ${\sigma }_g^2\approx 0.24204$. Why? $\mathbb{E}[e^{2U}]={\int }_0^1{e}^{2x}dx=\frac{e^2-1}{2}$; variance is second moment minus mean squared, the standard $\mathrm{Var}=\mathbb{E}[g^2]-(\mathbb{E}[g]{)}^2$.

Step 3 — Antithetic estimator. For each draw use $h(U)=\frac{1}{2}(e^U+e^{1-U})$, then average the $h$'s. Why this step? $e^U$ and $e^{1-U}$ are negatively correlated (when one is large the other is small), so their average has smaller variance while keeping the same mean $\theta$. This is the workhorse of Importance Sampling-style variance reduction.

Step 4 — Antithetic variance (per pair). $\mathrm{Var}(h)=\frac{1}{4}(\mathrm{Var}(e^U)+\mathrm{Var}(e^{1-U})+2\mathrm{Cov}(e^U,e^{1-U})).$ By symmetry the two variances equal ${\sigma }_g^2$; the covariance is negative: $\mathrm{Cov}(e^U,e^{1-U})=\mathbb{E}[e^U{e}^{1-U}]-{\theta }^2=e-(e-1{)}^2=2.71828-2.95249=-\mathrm{0.23421.}$ So $\mathrm{Var}(h)=\frac{1}{4}(2\times 0.24204+2\times (-0.23421))=\frac{1}{4}(0.01566)=\mathrm{0.003914.}$ Why $\mathbb{E}[e^U{e}^{1-U}]=e$? Because $e^U{e}^{1-U}=e^{U+(1-U)}=e^1$ is constant — the exponents always add to 1! That's why the correlation is so strongly negative and the variance nearly cancels.

Verify: plain per-sample variance $\approx 0.2420$; antithetic per-pair variance $\approx 0.00391$. Ratio $\approx 62\times$ smaller variance — a massive win. Same target $\theta =e-1\approx 1.71828$. ✔

---

Active Recall