Monte Carlo Simulation (Empirical Probabilities)

Context: FIT2086_MOC · W10 opens the simulation half of LO4 — replace mathematical ingenuity with computational power · convergence is the WLLN from Expectations and Covariance (FIT2086) · draws come from a PRNG via R Simulation and Random Sampling · applied to estimators in Monte Carlo Estimator Comparison, to the data itself in Bootstrap and Permutation Tests Parent Framework: Statistical Modelling and Inference

Quick Revision

• 🎯 Objective: any $\mathbb{P}((X_1,\dots ,X_p)\in Z)$ ➔ simulate $m$ draws ➔ average the indicators ➔ WLLN guarantees convergence ➔ answers problems with no closed form.
• 📦 Core Components: generate ➔ check $I(\cdot \in Z)$ ➔ average ➔ converge | resolution $1/m$ ➔ need $m\gg 1/p$.
• ⚡ Key Constraint: an event rarer than $1/m$ is usually never observed ⟹ estimate $=0$; size $m$ so the event is expected many times.

📝 How It Works

1. The General Problem

• Target ➔ RVs $X_1,\dots ,X_p\in \mathcal{X}$ and a set $Z\subseteq {\mathcal{X}}^p$ ➔ find $\mathbb{P}((X_1,\dots ,X_p)\in Z)$.
• Hard by hand ➔ $X_1,X_2\sim N(0,1)$: $\mathbb{P}(\sqrt{|X_1|}+\sqrt{|X_2|}>2)$ · $X_1\sim \text{Poi}({\lambda }_1)$, $X_2\sim \text{Poi}({\lambda }_2)$: $\mathbb{P}(X_1-X_2^{3/2}>0)$ ➔ answerable analytically but hard; trivial to pose ones that are impossible. Simulation always gives an answer.
• Buffon’s needle (18th c.) ➔ $\mathbb{P}(\text{needle crosses a strip line})$ needs integral geometry; simulation needs nothing: drop, record, repeat, estimate $\frac{\#\text{crossings}}{\#\text{drops}}$ (e.g. $7/20=0.35$).

2. The Algorithm

1. Generate ➔ $(X_1^{(i)},\dots ,X_p^{(i)})$ from their joint distribution, $i=1,\dots ,m$.
2. Check ➔ $I((X_1^{(i)},\dots ,X_p^{(i)})\in Z)=1$ if the condition holds, $0$ otherwise.
3. Average ➔ ${\hat{\theta }}_m=\frac{1}{m}{\sum }_{i=1}^{m}I((X_1^{(i)},\dots ,X_p^{(i)})\in Z)$ — the empirical probability.
4. Converge ➔ larger $m$ ⟹ more accurate.

3. Resolution and Rare Events

• Resolution ➔ $1/m$ = the smallest non-zero probability the estimate can take; the estimate moves in steps of $1/m$.
• Rare events ➔ $p<1/m$ cannot be estimated well ⟹ choose $m\gg 1/p$.
• Break Q1 ➔ event seen $23$ times in $m=1000$ ⟹ $\hat{p}=0.023$; $\hat{p}\to p$ as $m\to \infty$.

4. Beyond Probabilities

• Any distribution summary ➔ x $=m$ draws from $p(x)$: $\mathbb{E}[X]\approx$ mean(x) · $\mathbb{V}[X]\approx$ var(x) · $Q(0.2)\approx$ quantile(x, 0.2) · $p(x)\approx$ hist(x, probability = TRUE) ➔ all sharpen as $m$ grows.

🧮 Proof Blueprint

• Theorem ➔ ${\hat{\theta }}_m=\frac{1}{m}{\sum }_{i=1}^{m}I({\mathbf{X}}^{(i)}\in Z)\to \theta =\mathbb{P}(\mathbf{X}\in Z)$ as $m\to \infty$.
• Strategy ➔ each indicator is a Bernoulli RV ⟹ ${\hat{\theta }}_m$ is a sample mean ⟹ WLLN.
• Derivation Steps:

$$\begin{array}{rl}{B}_i& =I({\mathbf{X}}^{(i)}\in Z)\in \{0,1\},\mathbb{P}(B_i=1)=\mathbb{P}(\mathbf{X}\in Z)=\theta \Rightarrow B_i\stackrel{iid}{\sim }Be(\theta )\\ {\hat{\theta }}_m& =\frac{1}{m}\sum _{i=1}^m{B}_i\text{(the mean of }m\text{ Bernoulli RVs)}\\ {\hat{\theta }}_m& \to \mathbb{E}[B_i]=\theta \text{as }m\to \infty \text{(weak law of large numbers)}\end{array}$$

• Q.E.D. ➔ the empirical probability is consistent for $\theta$; being a Bernoulli mean, its variance is $\theta (1-\theta )/m$ (Studio 4 drill) — why small-$m$ runs wobble.

⚙️ Core Implementation

🔹 Monte Carlo in two lines of R

Code / Layout Details

X = rnorm(m, 0, 1)          # m draws of X ~ N(0,1)
mean(sqrt(abs(X)) > 1)      # logical vector -> mean = proportion TRUE = empirical probability

💡 Common Mistake: hist(x, Probability = T) as printed on the slide ➔ R’s argument is lowercase probability (or freq = FALSE); the capitalised name is not matched and the y-axis stays as counts.

📊 Exam Execution Trace & Applied Exercises

Manual Execution Trace

$X\sim N(0,1)$, target $\mathbb{P}(\sqrt{|X|}>1)$:

$m$	Resolution $1/m$	Estimate
$100$	$0.01$	$0.3300$
$1,000$	$0.001$	$0.3310$
$10,000$	$0.0001$	$0.3179$
$100,000$	$0.00001$	$0.3176$
$10^9$	$10^{-9}$	$0.3172$

• Reading ➔ convergence is not monotone ($m=1000$ is no closer than $m=100$) — each run is random; only the typical error shrinks.

Applied Exercise

Problem: derive the exact value the simulation chases: $X\sim N(0,1)$, $Q=\sqrt{|X|}$, find $\mathbb{P}(Q>1)$ — the slide’s CDF route, then the shortcut.

$$\begin{array}{rl}{F}_Q(q)& =\mathbb{P}(|X|\le q^2)=2\Phi (q^2)-1\\ p(q)& =2\phi (q^2)\cdot 2q={\left(\frac{2^3}{\pi }\right)}^{1/2}q\exp \left(-\frac{q^4}{2}\right),q\ge 0\\ \mathbb{P}(Q>1)& ={\int }_1^{\infty }p(q)dq=1-\mathrm{erf}(1/\sqrt{2})\\ \text{shortcut: }\mathbb{P}(Q>1)& =\mathbb{P}(|X|>1)=2(1-\Phi (1))=2(1-0.8413)=0.3173\end{array}$$

Final Extracted Output: $\mathbb{P}(Q>1)\approx 0.317$. The slides print $0.3171$; $1-\mathrm{erf}(1/\sqrt{2})=0.31731$ — a rounding slip on the slide, and the $10^9$ run ($0.3172$) agrees with $0.3173$ to within simulation error.

⚠️ Common Mistakes

• 💡 $m$ too small for a rare event ➔ $\hat{p}=0$ is not evidence that $p=0$; it means $p\lesssim 1/m$.
• 💡 Sampling each $X_j$ from the wrong law ➔ draw from the joint distribution; independent marginals are only valid when the $X_j$ are independent.
• 💡 No seed ➔ the estimate changes every run; set.seed() first when the number must be reproduced (Pseudo-Random Number Generators).

🧠 Active Recall

Why does the proportion of simulated draws landing in $Z$ converge to $\mathbb{P}(\mathbf{X}\in Z)$?

Answer

• Short answer: each indicator is a $Be(\theta )$ RV with $\theta =\mathbb{P}(\mathbf{X}\in Z)$, so the proportion is a sample mean, which the WLLN sends to $\theta$.
• Why: Indicator ➔ Bernoulli ➔ WLLN ➔ $\frac{1}{m}\sum I({\mathbf{X}}^{(i)}\in Z)\to \mathbb{E}[I]=\theta$ as $m\to \infty$.

You need $\mathbb{P}(A)\approx 10^{-5}$ and simulate $m=10^4$ times. What goes wrong, and what $m$ do you need?

Answer

• Short answer: the expected count is $mp=0.1$, so the estimate is almost always $0$ — below the resolution $1/m=10^{-4}$; take $m\gg 10^5$ (e.g. $10^7$, expecting $\approx 100$ hits).
• Why: Resolution ➔ the estimate lives on the grid $\{0,\frac{1}{m},\frac{2}{m},\dots \}$; it can only resolve $p$ once the event is observed many times.