Algorithms Teaching Series · Advanced

蒙地卡羅方法完整教材

Monte Carlo Method

Theory · Algorithms · Applications

用「擲骰子」解決確定性的問題：從大數法則與中央極限定理出發，走過隨機數生成、蒙地卡羅積分、變異數縮減、馬可夫鏈蒙地卡羅（MCMC），到 Miller–Rabin、Freivalds、Karger、模擬退火、蒙地卡羅樹搜尋與選擇權定價—— 十七個可執行的 C++ 實例，每一個都與理論值對照。

配套程式：monte_carlo.cpp（C++17，單檔、固定隨機種子、結果可重現）

2026 年 9 月

你想知道一座不規則形狀的湖有多大，但沒有任何測量工具。有一個笨到極點、卻絕對正確的辦法：站在直升機上，對著包含這座湖的一塊方形土地均勻地亂丟石頭，然後數「掉進湖裡」的石頭比例。方形土地的面積你知道，乘上比例就是湖的面積。丟越多石頭，估得越準；而且不管湖的形狀有多詭異、有多少個島，方法完全不用改。這就是蒙地卡羅方法：用隨機抽樣回答一個確定性的問題。

一句話定義

把想求的量 $I$ 寫成某個隨機變數的期望值 $I=\mathbb{E}[Y]$，然後用樣本平均來估計它： \[\hat I_N = \frac{1}{N}\sum_{i=1}^{N} Y_i,\qquad Y_1,\dots,Y_N \text{ 獨立同分布（i.i.d.）}.\] 大數法則保證 $\hat I_N\to I$；中央極限定理告訴我們誤差以 $O(1/\sqrt N)$ 的速度縮小，且此速率與問題的維度無關。

這個定義看起來平淡，威力卻在三個地方：

不挑問題結構：只要能「抽樣」與「評分」，就能跑——不需要可微、不需要凸、不需要封閉解。
不怕高維度：傳統數值積分的誤差隨維度指數惡化（維度災難），蒙地卡羅的 $1/\sqrt N$ 對 $d=1$ 與 $d=1000$ 一視同仁。
自帶誤差條：每次估計都附贈一個標準誤差與信賴區間，你永遠知道自己「多不準」。

歷史：從投針到曼哈頓計畫

1777 年，Buffon 投針：Georges-Louis Leclerc de Buffon 提出「往有平行線的地板丟針，針壓到線的機率是 $2l/(\pi d)$」——這是第一個用隨機實驗算 $\pi$ 的方法（第 4.2 節會實作）。
1930 年代，Fermi：Enrico Fermi 在研究中子擴散時，用手搖機械計算器做隨機抽樣，據說還以此「猜」實驗結果讓同事大吃一驚。
1946 年，Ulam 與 von Neumann：Stanislaw Ulam 病中玩 Canfield 接龍，想算贏牌機率卻算不出來，靈機一動——「乾脆玩一百次數一數」。他和 John von Neumann 意識到這招可以解決曼哈頓計畫中的中子傳輸問題，並在 ENIAC 上實作。
命名：Nicholas Metropolis 建議以摩納哥的蒙地卡羅賭場命名——Ulam 的叔叔常向親戚借錢去那裡賭博。
1953 年，Metropolis 演算法：Metropolis、Rosenbluth 夫婦、Teller 夫婦發表的論文開啟了馬可夫鏈蒙地卡羅（MCMC），1970 年 Hastings 推廣為 Metropolis–Hastings；1984 年 Geman 兄弟提出 Gibbs 抽樣；1990 年代 MCMC 席捲貝氏統計。
今天：從核反應爐設計、電影渲染（路徑追蹤）、金融衍生品定價、天氣集合預報、粒子物理，到 AlphaGo 的蒙地卡羅樹搜尋，蒙地卡羅是二十世紀最有影響力的演算法之一。

三種用途

蒙地卡羅方法在實務上分成三個家族，本教材依序涵蓋：

Monte Carlo 演算法 vs. Las Vegas 演算法

在演算法理論中，「蒙地卡羅演算法」有更狹義的意思——使用隨機性、執行時間有界、但答案可能錯的演算法。它的對照組是 Las Vegas 演算法：答案一定對、但執行時間是隨機變數（例如隨機化快速排序）。

	Monte Carlo 演算法	Las Vegas 演算法
正確性	有界的錯誤機率（可放大縮小）	永遠正確
執行時間	確定（或有界）	隨機，期望值有界
典型例子	Miller–Rabin、Freivalds、Karger	隨機化快速排序、隨機化選擇
轉換	若能驗證答案 $\Rightarrow$ 反覆執行至正確即成 Las Vegas	限時強制輸出 $\Rightarrow$ 成 Monte Carlo

本教材標題的「蒙地卡羅方法」取廣義：所有以隨機抽樣為核心的數值與演算法技術。第 7 章討論的「Monte Carlo 型隨機化演算法」則是狹義——兩者精神相同：用隨機性換取效率，用重複換取正確性。

何時該用、何時不該用

該用：高維積分（$d\gtrsim 4$）、複雜幾何、只知道未正規化的機率密度（貝氏後驗）、離散結構的計數與抽樣、沒有封閉解的隨機系統（排隊、金融路徑）、對手不確定的決策（博弈）。
不該用：低維光滑函數（高斯求積可達指數收斂，蒙地卡羅的 $1/\sqrt N$ 望塵莫及）、需要精確到 $10^{-12}$ 的量（需要 $N\sim 10^{24}$ 個樣本）、有多項式時間確定性演算法的問題。

數學基礎：為什麼隨機抽樣是對的

本章建立整套方法的理論地基。三個定理支撐一切：大數法則（會收斂）、中央極限定理（誤差多大）、集中不等式（需要多少樣本）。

隨機變數、期望值與變異數

定義 2.1 (期望值與變異數). 設 $Y$ 為隨機變數，其期望值與變異數為 \[\mu=\mathbb{E}[Y]=\int y\,p(y)\,dy,\qquad \sigma^2=\operatorname{Var}(Y)=\mathbb{E}\big[(Y-\mu)^2\big]=\mathbb{E}[Y^2]-\mu^2 .\] 對 i.i.d. 樣本 $Y_1,\dots,Y_N$，樣本平均 $\bar Y_N=\frac1N\sum Y_i$ 滿足 \[\mathbb{E}[\bar Y_N]=\mu,\qquad \operatorname{Var}(\bar Y_N)=\frac{\sigma^2}{N}.\]

第二式是整個蒙地卡羅方法的心臟：平均 $N$ 個獨立樣本，變異數縮小 $N$ 倍，標準差縮小 $\sqrt N$ 倍。

定義 2.2 (估計量的品質). 估計量 $\hat I$ 對目標 $I$：

無偏（unbiased）：$\mathbb{E}[\hat I]=I$。樣本平均是 $\mathbb{E}[Y]$ 的無偏估計量。
一致（consistent）：$\hat I_N\to I$（機率收斂）當 $N\to\infty$。
均方誤差：$\mathrm{MSE}(\hat I)=\mathbb{E}[(\hat I-I)^2]=\underbrace{(\mathbb{E}[\hat I]-I)^2}_{\text{偏差}^2}+\operatorname{Var}(\hat I)$。

無偏估計量的 MSE 就是變異數——所以第 5 章「變異數縮減」等同「提高精度」。

大數法則：樣本平均會收斂

定理 2.1 (弱大數法則). 設 $Y_1,Y_2,\dots$ i.i.d.，$\mathbb{E}[Y_i]=\mu$，$\operatorname{Var}(Y_i)=\sigma^2<\infty$。則對任意 $\varepsilon>0$， \[\mathbb{P}\big(|\bar Y_N-\mu|\ge\varepsilon\big)\le\frac{\sigma^2}{N\varepsilon^2}\xrightarrow[N\to\infty]{}0 .\]

Proof. Chebyshev 不等式：對任何有限變異數的隨機變數 $X$，$\mathbb{P}(|X-\mathbb{E}X|\ge\varepsilon)\le\operatorname{Var}(X)/\varepsilon^2$。代入 $X=\bar Y_N$，$\operatorname{Var}(\bar Y_N)=\sigma^2/N$ 即得。 ◻

定理 2.2 (強大數法則，Kolmogorov). 同上假設（其實只需 $\mathbb{E}|Y_i|<\infty$），則 $\bar Y_N\to\mu$ 幾乎必然（almost surely）： $\mathbb{P}\big(\lim_{N\to\infty}\bar Y_N=\mu\big)=1$。

強大數法則說的是：你跑的那一條樣本路徑（不只是「大多數路徑」）終究會收斂到 $\mu$。這是為什麼我們敢只跑一次程式就相信結果。

中央極限定理：誤差有多大

定理 2.3 (中央極限定理，CLT). 同上假設且 $0<\sigma^2<\infty$，則 \[\frac{\sqrt N\,(\bar Y_N-\mu)}{\sigma}\ \xrightarrow{\ d\ }\ \mathcal N(0,1), \qquad\text{即}\qquad \bar Y_N\ \approx\ \mathcal N\!\Big(\mu,\ \frac{\sigma^2}{N}\Big)\ \text{（$N$ 大時）}.\]

以樣本變異數 $s^2=\frac{1}{N-1}\sum(Y_i-\bar Y_N)^2$ 估計 $\sigma^2$，得標準誤差（standard error） \[\mathrm{s.e.}(\hat I_N)=\frac{s}{\sqrt N}, \qquad 95\%\ \text{信賴區間} \approx \hat I_N\pm 1.96\,\frac{s}{\sqrt N}.\] 誤差 $\propto 1/\sqrt N$：想多一位小數（誤差縮小 10 倍），樣本要多 100 倍。

好消息：速率與維度無關、與被積函數是否光滑無關。壞消息：速率很慢。從 $10^{-3}$ 到 $10^{-6}$ 精度要多花一百萬倍時間。因此蒙地卡羅的實務藝術不是「多丟樣本」，而是縮小 $\sigma$（第 5 章）。

集中不等式：需要多少樣本才夠

CLT 是漸近敘述；若要「保證」（非漸近）的樣本數，用集中不等式。

定理 2.4 (Hoeffding 不等式). 設 $Y_i$ 獨立且 $a\le Y_i\le b$，則 \[\mathbb{P}\big(|\bar Y_N-\mu|\ge\varepsilon\big)\le 2\exp\!\Big(-\frac{2N\varepsilon^2}{(b-a)^2}\Big).\]

推論 2.1 (樣本數公式). 要以至少 $1-\delta$ 的機率保證 $|\bar Y_N-\mu|<\varepsilon$，只需 \[N\ \ge\ \frac{(b-a)^2}{2\varepsilon^2}\,\ln\frac{2}{\delta}.\]

範例 2.1 (估計 $\pi$ 要丟幾個點). 用 hit-or-miss 法（第 4.1 節），$Y_i=4\cdot\mathbf{1}\{\text{落入圓內}\}\in[0,4]$。要 99% 信心下誤差 $<0.01$： $N\ge\frac{16}{2\times 10^{-4}}\ln 200\approx 4.2\times10^{5}$。 Hoeffding 是保守上界；用 CLT 估（$\sigma=4\sqrt{p(1-p)}\approx1.64$）約需 $N\approx(2.58\times1.64/0.01)^2\approx1.8\times10^5$。兩者同一量級，這正是「$\varepsilon$ 縮小 10 倍、$N$ 放大 100 倍」的體現。

與確定性數值積分的比較：維度災難

計算 $I=\int_{[0,1]^d}f(x)\,dx$。用每維 $n$ 個格點的梯形法，共 $N=n^d$ 個函數評估，誤差為 $O(n^{-2})=O(N^{-2/d})$；Simpson 法為 $O(N^{-4/d})$。蒙地卡羅則是 $O(N^{-1/2})$，與 $d$ 無關。

維度 $d$	梯形法	Simpson	蒙地卡羅	勝者
1	$N^{-2}$	$N^{-4}$	$N^{-1/2}$	確定性方法（壓倒性）
2	$N^{-1}$	$N^{-2}$	$N^{-1/2}$	確定性方法
4	$N^{-1/2}$	$N^{-1}$	$N^{-1/2}$	梯形法打平；Simpson 略勝
8	$N^{-1/4}$	$N^{-1/2}$	$N^{-1/2}$	Simpson 打平
20	$N^{-0.1}$	$N^{-0.2}$	$N^{-1/2}$	蒙地卡羅（壓倒性）
100	$N^{-0.02}$	$N^{-0.04}$	$N^{-1/2}$	只有蒙地卡羅可行

$d=100$ 時，梯形法每維只取 2 個點就需要 $2^{100}\approx10^{30}$ 次函數評估——宇宙年齡都不夠。蒙地卡羅用 $10^6$ 個點就給你三位有效數字。這就是維度災難（curse of dimensionality），也是蒙地卡羅存在的根本理由。

蒙地卡羅的誤差是 $\sigma/\sqrt N$，常數 $\sigma$ 是被積函數的標準差，與光滑度無關；確定性方法的常數涉及高階導數，函數不光滑（有跳躍、尖點）時會爆掉。因此即使在低維，不連續的被積函數（如選擇權的 $\max(S-K,0)$）也常用蒙地卡羅。

蒙地卡羅積分的兩個基本估計量

設 $f:[0,1]\to[0,1]$，$I=\int_0^1 f(x)\,dx$。

定義 2.3 (Hit-or-miss 估計量). 在單位正方形內均勻丟 $N$ 個點 $(X_i,Y_i)$， \[\hat I^{\mathrm{HM}}_N=\frac1N\sum_{i=1}^N \mathbf{1}\{Y_i\le f(X_i)\}, \qquad \operatorname{Var}\big(\hat I^{\mathrm{HM}}_N\big)=\frac{I(1-I)}{N}.\]

定義 2.4 (平均值（crude Monte Carlo）估計量). 抽 $X_i\sim U(0,1)$， \[\hat I^{\mathrm{MV}}_N=\frac1N\sum_{i=1}^N f(X_i), \qquad \operatorname{Var}\big(\hat I^{\mathrm{MV}}_N\big)=\frac{1}{N}\Big(\int_0^1 f^2\,dx - I^2\Big).\]

定理 2.5. 對 $0\le f\le1$，恆有 $\operatorname{Var}(\hat I^{\mathrm{MV}}_N)\le\operatorname{Var}(\hat I^{\mathrm{HM}}_N)$。

Proof. 因 $0\le f\le 1$，$f^2\le f$，故 $\int f^2\le\int f=I$，於是 $\int f^2-I^2\le I-I^2=I(1-I)$。 ◻

直觀：hit-or-miss 把每個樣本壓成 0/1 一個位元，丟掉了 $f(X_i)$ 的大小資訊；平均值法保留了它。第 4.3 節的實驗中，兩者變異數差了 4.7 倍。教學上 hit-or-miss 最直觀，實務上幾乎總是用平均值法。

隨機數生成與抽樣技術

蒙地卡羅的一切都建立在「能產生想要的分布的樣本」之上。本章由淺入深：均勻偽隨機數 $\to$ 逆變換 $\to$ 常態分布 $\to$ 拒絕抽樣 $\to$ 準隨機數。

偽隨機數生成器（PRNG）

電腦是確定性機器，「隨機數」其實是偽隨機（pseudo-random）：由種子（seed）出發、確定性遞推、但統計性質像 i.i.d. 均勻分布的數列。

線性同餘法 LCG：$X_{n+1}=(aX_n+c)\bmod m$。簡單、快、但週期 $\le m$ 且高維格點結構明顯（著名的 RANDU 災難）。只適合教學。
Mersenne Twister（MT19937）：週期 $2^{19937}-1$，623 維等分布，C++ 的 std::mt19937 / mt19937_64。科學計算的長期標準，缺點是狀態大（2.5 KB）、非密碼安全。
xorshift / xoshiro / PCG：狀態小、速度極快、統計品質優良，現代模擬與遊戲的主流選擇。
密碼學安全 PRNG（CSPRNG）：ChaCha20、AES-CTR——輸出不可預測，用於安全用途；模擬不需要（且較慢）。

std::mt19937_64 配固定種子 20260912，搭配 std::uniform_real_distribution 與 std::normal_distribution。固定種子讓每次執行結果完全一致——除錯、審稿、教學都需要。正式實驗時再用 std::random_device 播種並記錄下種子。

(1) 部分實作週期只有 $2^{31}$、低位元品質極差；(2) rand() % n 有模餘偏差（modulo bias）——除非 $n$ 整除 RAND_MAX+1，否則小數字被抽到的機率略高；(3) 全域狀態、不可重入。C++11 以後請一律用 <random>。

逆變換法（Inverse Transform Sampling）

定理 3.1 (逆變換). 設 $F$ 為連續分布的累積分布函數（CDF），$F^{-1}(u)=\inf\{x:F(x)\ge u\}$ 為其廣義逆。若 $U\sim U(0,1)$，則 $X=F^{-1}(U)$ 的 CDF 恰為 $F$。

Proof. $\mathbb{P}(X\le x)=\mathbb{P}(F^{-1}(U)\le x)=\mathbb{P}(U\le F(x))=F(x)$，第二個等號用到 $F$ 單調不減。 ◻

範例 3.1 (指數分布). $F(x)=1-e^{-\lambda x}$，$F^{-1}(u)=-\ln(1-u)/\lambda$。因 $1-U$ 亦為 $U(0,1)$，實務寫成 $X=-\ln U/\lambda$。這是排隊模擬、事件間隔時間的基本工具。

範例 3.2 (離散分布). 機率 $p_1,\dots,p_k$，累積 $c_j=\sum_{i\le j}p_i$。抽 $U$，回傳最小的 $j$ 使 $c_j\ge U$（二分搜尋 $O(\log k)$）。這也是遺傳演算法「輪盤選擇」的實作。若要 $O(1)$，用 Walker 的別名法（alias method）。

逆變換法的限制：需要 $F^{-1}$ 可算。常態分布的 $\Phi^{-1}$ 沒有封閉解，因此有下一節。

常態分布：Box–Muller 變換

定理 3.2 (Box–Muller). 設 $U_1,U_2\sim U(0,1)$ 獨立，則 \[Z_0=\sqrt{-2\ln U_1}\,\cos(2\pi U_2),\qquad Z_1=\sqrt{-2\ln U_1}\,\sin(2\pi U_2)\] 為兩個獨立的標準常態隨機變數。

證明的想法：二維標準常態 $(Z_0,Z_1)$ 的密度在極座標 $(R,\Theta)$ 下分解為兩個獨立因子： \[\frac{1}{2\pi}e^{-(z_0^2+z_1^2)/2}\,dz_0\,dz_1 =\underbrace{\tfrac12 e^{-r^2/2}\,d(r^2)}_{R^2\sim\mathrm{Exp}(1/2)}\cdot\underbrace{\tfrac{1}{2\pi}\,d\theta}_{\Theta\sim U(0,2\pi)},\] 而 $R^2=-2\ln U_1$ 正是指數分布的逆變換。C++ 的 std::normal_distribution 內部使用改良版（Marsaglia 極座標法，避開三角函數）。

拒絕抽樣（Rejection Sampling）

當目標密度 $p(x)$ 只能逐點計算（甚至只知道到常數倍）、無法逆變換時，用「包絡分布」$q$ 幫忙。

目標密度 $p$，提議密度 $q$（容易抽樣），常數 $M$ 使 $p(x)\le M\,q(x)$ 對所有 $x$ 成立一個服從 $p$ 的樣本抽 $X\sim q$，抽 $U\sim U(0,1)$ $X$

定理 3.3. 演算法回傳的 $X$ 服從 $p$；每一輪被接受的機率為 $1/M$，故期望迭代次數為 $M$。

Proof. $\mathbb{P}(X\in dx,\ \text{接受})=q(x)\cdot\frac{p(x)}{Mq(x)}\,dx=\frac{p(x)}{M}dx$。對 $x$ 積分得 $\mathbb{P}(\text{接受})=1/M$。條件機率 $\mathbb{P}(X\in dx\mid\text{接受})=\frac{p(x)/M}{1/M}dx=p(x)\,dx$。 ◻

$M$ 是包絡與目標之間的「浪費比例」。高維時即使 $q$ 每一維都只比 $p$ 胖 10%，$M\approx1.1^d$ 仍指數成長——$d=100$ 時要丟一萬個才中一個。這就是高維抽樣必須改用 MCMC（第 6 章）的原因。

準蒙地卡羅：低差異序列

若把 i.i.d. 均勻點換成刻意「均勻鋪滿」的確定性點列（Halton、Sobol、Faure 序列），Koksma–Hlawka 不等式給出誤差上界 \[\big|\hat I_N-I\big|\ \le\ V(f)\cdot D^*_N,\qquad D^*_N=O\!\Big(\frac{(\ln N)^d}{N}\Big),\] 其中 $V(f)$ 是 $f$ 的（Hardy–Krause）變差，$D^*_N$ 是點列的星差異。$N^{-1}$ 比 $N^{-1/2}$ 快得多——代價是失去「自帶誤差條」（點不是隨機的，沒有標準誤差可算），且 $d$ 大時 $(\ln N)^d$ 常數會咬人。金融業普遍用「隨機化準蒙地卡羅（RQMC）」兩全其美。本教材聚焦真正的隨機方法，準蒙地卡羅留作延伸閱讀。

蒙地卡羅積分：從 $\pi$ 到高維球

本章用四個經典例子把第 2 章的理論「跑出來」。所有數字來自配套程式 monte_carlo.cpp（第 10 章），固定種子、可重現。

估計 $\pi$：hit-or-miss

在 $[0,1]^2$ 均勻丟點，落入四分之一單位圓 $x^2+y^2\le1$ 的機率是 $\pi/4$。設命中比例 $\hat p$，則 $\hat\pi=4\hat p$，且 \[\mathrm{s.e.}(\hat\pi)=4\sqrt{\frac{\hat p(1-\hat p)}{N}}\approx\frac{1.64}{\sqrt N}.\]

=== Example 1: estimating pi by hit-or-miss, error ~ 1/sqrt(N) ===
N         estimate    |error|     std.err     95% CI
100       3.120000    0.021593    0.165699    [2.795231, 3.444769]
1000      3.204000    0.062407    0.050501    [3.105017, 3.302983]
10000     3.137200    0.004393    0.016452    [3.104954, 3.169446]
100000    3.132320    0.009273    0.005213    [3.122102, 3.142538]
1000000   3.142760    0.001167    0.001641    [3.139543, 3.145977]
10000000  3.141712    0.000119    0.000519    [3.140694, 3.142729]

藍線是理論預測的 $1/\sqrt N$（對數座標下斜率 $-1/2$ 的直線）。紅點是「真實誤差」，它是隨機的——有時比 s.e. 小（運氣好），有時大（$N=1000$ 時 1.2 個 s.e.），但全部落在灰色 95% 虛線之下。真實誤差不會單調下降（$N=10^5$ 比 $N=10^4$ 還差），這不是 bug，是隨機性的正常表現；單調下降的是誤差的分布寬度。

Buffon 投針

地板有間距 $d$ 的平行線，丟一根長 $l\le d$ 的針。設針心到最近線的距離 $D\sim U(0,d/2)$、針與線夾角 $\Theta\sim U(0,\pi/2)$，兩者獨立。針與線相交 $\iff D\le\frac{l}{2}\sin\Theta$，故 \[\mathbb{P}(\text{相交})=\frac{1}{(d/2)(\pi/2)}\int_0^{\pi/2}\!\!\int_0^{\frac{l}{2}\sin\theta}dD\,d\theta =\frac{4}{\pi d}\cdot\frac{l}{2}\int_0^{\pi/2}\sin\theta\,d\theta=\frac{2l}{\pi d}.\] 取 $l=1,d=2$ 得 $\mathbb{P}=1/\pi$，於是 $\hat\pi=1/\hat p$。

=== Example 2: Buffon's needle (l = 1, d = 2), P(cross) = 2l/(pi d) = 1/pi ===
N = 1000       P(cross) = 0.32100  ->  pi ~ 3.11526
N = 100000     P(cross) = 0.31664  ->  pi ~ 3.15816
N = 10000000   P(cross) = 0.31856  ->  pi ~ 3.13917

注意 $\hat\pi=1/\hat p$ 是 $\hat p$ 的非線性函數，因此有偏（Jensen 不等式：$\mathbb{E}[1/\hat p]>1/\mathbb{E}[\hat p]$），偏差 $O(1/N)$，遠小於標準誤差 $O(1/\sqrt N)$，實務可忽略；但這提醒我們——「先估機率再取倒數」與「直接估想要的量」不完全相同。

Hit-or-miss vs. 平均值法：$\int_0^1 e^{-x^2}dx$

這個積分沒有初等封閉解（$=\frac{\sqrt\pi}{2}\operatorname{erf}(1)\approx0.746824$），是典型的「非要數值方法」的例子。第 2.6 節證明了平均值法變異數必較小，實驗印證：

=== Example 3: hit-or-miss vs mean-value estimator for I = int_0^1 exp(-x^2) dx ===
exact        = 0.746824
hit-or-miss  = 0.746643  (per-sample var 0.189167, s.e. 0.000435)
mean-value   = 0.746708  (per-sample var 0.040382, s.e. 0.000201)
variance ratio (hit-or-miss / mean-value) = 4.68

理論值：$I(1-I)=0.1891$；$\int_0^1 f^2\,dx-I^2=\int_0^1 e^{-2x^2}dx-I^2=\frac{\sqrt{2\pi}}{4}\operatorname{erf}(\sqrt2)-I^2\approx0.5981-0.5577=0.0404$——與樣本變異數一致。變異數比 4.68 意味著 hit-or-miss 要用 4.68 倍的樣本才能達到同樣精度。

高維單位球的體積：維度災難的第一眼

$d$ 維單位球體積的精確公式 $V_d=\pi^{d/2}/\Gamma(d/2+1)$。用 hit-or-miss：在 $[-1,1]^d$（體積 $2^d$）均勻丟點，$\hat V_d=2^d\hat p$。

=== Example 4: volume of the unit d-ball, N = 1e6 points in [-1,1]^d ===
d   MC volume   exact       rel.err     hit rate
2   3.14221     3.14159     0.0002      7.86e-01
4   4.93438     4.93480     0.0001      3.08e-01
6   5.14413     5.16771     0.0046      8.04e-02
8   4.08269     4.05871     0.0059      1.59e-02
10  2.46579     2.55016     0.0331      2.41e-03
12  1.26157     1.33526     0.0552      3.08e-04

三個觀察：

$V_d$ 在 $d=5$ 附近達最大後趨近於 0——高維球「幾乎沒有體積」，這是高維幾何的反直覺事實之一。
命中率 $V_d/2^d$ 指數下降：$d=12$ 時一百萬個點只有 308 個命中。比例估計的相對標準誤差 $\sqrt{(1-p)/(Np)}\approx 1/\sqrt{Np}$——命中數才是有效樣本數。$d=20$ 時命中率約 $2.5\times10^{-8}$，一百萬個點很可能一個都沒中。
這不是蒙地卡羅的失敗，而是 hit-or-miss 這種抽樣策略的失敗：我們把絕大多數樣本浪費在球外。解法是「往重要的地方多抽」——第 5 章的重要性抽樣、或第 6 章直接在球內漫步。

(1) 永遠附上標準誤差或信賴區間，沒有誤差條的蒙地卡羅數字是沒有意義的；(2) 報告樣本數與隨機種子，讓別人能重現；(3) 若估計的是極小機率，報告命中數而非只報比例——35 個命中和 35000 個命中的可信度天差地遠。

變異數縮減技術

兩個學生都想估計全校學生的平均身高。甲隨機抓 1000 人量；乙先按「男/女、各年級」分層，每層按比例抽人，總數同樣 1000。乙的估計穩得多——因為他把「層間差異」這部分變異數消滅了，只剩層內的隨機性。變異數縮減就是這類用已知結構換精度的技巧集合。

蒙地卡羅的誤差是 $\sigma/\sqrt N$。多花 100 倍時間只能換一位小數；但若能把 $\sigma^2$ 砍掉 100 倍，效果相同、成本幾乎為零。本章五種技術的效率量度都是同一個： \[\text{效率}=\frac{1}{\operatorname{Var}(\hat I)\times\text{每樣本成本}} .\]

對偶變數（Antithetic Variates）

用 $U$ 與 $1-U$ 配對：兩者都服從 $U(0,1)$，但負相關。估計量 $\hat I=\frac{1}{2}\big(f(U)+f(1-U)\big)$， \[\operatorname{Var}\Big(\frac{f(U)+f(1-U)}{2}\Big)=\frac{\operatorname{Var}f(U)+\operatorname{Cov}\big(f(U),f(1-U)\big)}{2}.\] 若 $f$ 單調，$\operatorname{Cov}(f(U),f(1-U))\le0$（同一個 $U$ 讓一項偏大、另一項偏小），故變異數嚴格小於用兩個獨立樣本平均的 $\operatorname{Var}f(U)/2$。對接近線性的 $f$，抵消幾乎完美——$f(x)=e^x$ 的實驗中變異數降到 3.2%。

$f$ 非單調（例如對稱的 $f(x)=(x-\tfrac12)^2$）時 $\operatorname{Cov}$ 可能為正，對偶變數反而變差。用之前先想清楚 $f$ 的形狀。

控制變數（Control Variates）

若有另一個隨機變數 $X$ 與 $Y=f(U)$ 高度相關、且 $\mathbb{E}[X]$ 已知，令 \[Y_c=Y-c\,(X-\mathbb{E}[X]).\] 仍是無偏的（$\mathbb{E}[Y_c]=\mathbb{E}[Y]$），而 \[\operatorname{Var}(Y_c)=\sigma_Y^2-2c\operatorname{Cov}(X,Y)+c^2\sigma_X^2 \ \xrightarrow{\ c^*=\operatorname{Cov}(X,Y)/\sigma_X^2\ }\ \operatorname{Var}(Y_{c^*})=\sigma_Y^2\,(1-\rho_{XY}^2).\] 相關係數 $\rho=0.99$ 就砍掉 98% 變異數。$c^*$ 實務上以樣本估計（會引入 $O(1/N)$ 偏差，可用先導樣本估 $c$ 避免）。實驗中以 $X=U$（$\mathbb{E}U=1/2$）當 $e^U$ 的控制變數，$\hat c=1.689$，變異數降到 1.6%。

金融：用有封閉解的幾何平均亞式選擇權當沒有封閉解的算術平均亞式選擇權的控制變數，兩者相關性 $>0.99$，變異數縮減上百倍（Kemna & Vorst 1990）。

分層抽樣（Stratified Sampling）

把定義域切成 $K$ 層 $S_1,\dots,S_K$，機率 $p_k=\mathbb{P}(X\in S_k)$，層內平均 $\mu_k$、層內變異數 $\sigma_k^2$。在第 $k$ 層抽 $n_k$ 個點，估計量 $\hat I=\sum_k p_k\bar Y_k$， \[\operatorname{Var}(\hat I)=\sum_{k=1}^K\frac{p_k^2\sigma_k^2}{n_k}.\]

定理 5.1 (比例配置永不吃虧). 取 $n_k=p_kN$（比例配置），則 \[\operatorname{Var}(\hat I_{\text{strat}})=\frac{1}{N}\sum_k p_k\sigma_k^2 =\frac{1}{N}\Big(\sigma^2-\sum_k p_k(\mu_k-\mu)^2\Big)\ \le\ \frac{\sigma^2}{N}.\]

Proof. 全變異數分解（law of total variance）：$\sigma^2=\mathbb{E}[\operatorname{Var}(Y\mid\text{層})]+\operatorname{Var}(\mathbb{E}[Y\mid\text{層}])=\sum p_k\sigma_k^2+\sum p_k(\mu_k-\mu)^2$。 ◻

被消滅的是層間變異數。實驗中把 $[0,1]$ 等切 100 層估 $\int e^x$，變異數降到樸素法的 $1.1\times10^{-4}$——四個數量級。這是因為 $e^x$ 光滑：每層內幾乎是常數，層內變異數 $\sim(1/K)^2$。（最佳配置 $n_k\propto p_k\sigma_k$ 更佳，稱 Neyman 配置。）

重要性抽樣（Importance Sampling）

這是最重要、也最危險的技術。

對任意密度 $q$ 滿足「$f(x)p(x)\ne0\Rightarrow q(x)>0$」， \[I=\mathbb{E}_p[f(X)]=\int f(x)\,p(x)\,dx=\int f(x)\frac{p(x)}{q(x)}\,q(x)\,dx=\mathbb{E}_q\big[f(X)\,w(X)\big], \qquad w(x)=\frac{p(x)}{q(x)}.\] 估計量：$\hat I^{\mathrm{IS}}_N=\frac1N\sum_{i=1}^N f(X_i)w(X_i)$，$X_i\sim q$。

定理 5.2 (最佳提議分布). $\operatorname{Var}_q\big(f(X)w(X)\big)$ 在 $q^*(x)\propto|f(x)|\,p(x)$ 時最小。若 $f\ge0$，則 $q^*=f p/I$ 給出零變異數。

Proof. 先寫出變異數：$\operatorname{Var}_q(fw)=\mathbb{E}_q[f^2w^2]-I^2=\int\frac{f^2p^2}{q}\,dx-I^2$，故只需最小化第一項。由 Cauchy–Schwarz， \[\Big(\int|f|\,p\,dx\Big)^2=\Big(\int\frac{|f|\,p}{\sqrt q}\cdot\sqrt q\,dx\Big)^2 \ \le\ \int\frac{f^2p^2}{q}\,dx\cdot\underbrace{\int q\,dx}_{=1},\] 左邊與 $q$ 無關，因此右邊第一項的下界為 $\big(\int|f|p\big)^2$，等號成立當 $\sqrt q\propto|f|p/\sqrt q$，即 $q\propto|f|\,p$。 ◻

零變異數的 $q^*$ 需要知道 $I$——雞生蛋。但它指明方向：把樣本抽到 $|f|p$ 大的地方，也就是「重要」的地方。

範例 5.1 (極小尾機率 $\mathbb{P}(Z>4)$，$Z\sim\mathcal N(0,1)$). 真值 $3.167\times10^{-5}$。樸素法 $f=\mathbf{1}\{z>4\}$，一百萬個樣本期望只有 32 個命中，相對誤差 $\approx1/\sqrt{32}\approx18\%$。改從 $q=\mathcal N(4,1)$ 抽（一半的樣本落在 $z>4$），權重 \[w(z)=\frac{\varphi(z)}{\varphi(z-4)}=\exp\!\Big(-\frac{z^2}{2}+\frac{(z-4)^2}{2}\Big)=e^{-4z+8}.\]

=== Example 6: importance sampling for the tail probability P(Z > 4) ===
exact                 = 3.1671e-05
naive  (35 hits)     = 3.5000e-05  s.e. 5.9160e-06  rel.err 18.7%
importance sampling   = 3.1593e-05  s.e. 6.7166e-08  rel.err 0.212%
variance reduction factor = 7.76e+03

變異數縮減 7760 倍，等於樸素法多跑 7760 倍樣本。尾機率越小，效果越誇張——$\mathbb{P}(Z>6)\approx10^{-9}$ 用樸素法根本不可能，用 IS 幾秒鐘。

若 $q$ 的尾巴比 $p$ 輕（$w=p/q$ 在尾端爆炸），$\mathbb{E}_q[w^2]$ 可能無限大——估計量仍無偏，但變異數無限，CLT 失效，你會看到「大多數時候估得很準、偶爾蹦出一個天文數字」。診斷工具：有效樣本數 $\mathrm{ESS}=\dfrac{(\sum w_i)^2}{\sum w_i^2}$。$N=10^6$ 但 ESS $=50$ 表示只有 50 個樣本在做事。原則：$q$ 的尾巴永遠要比 $p$ 重（例如用 $t$ 分布代替常態）。

條件蒙地卡羅（Rao–Blackwell 化）

若 $Y$ 的一部分期望可以解析算出，就別抽它。設 $Y=f(X,Z)$，令 $g(X)=\mathbb{E}[f(X,Z)\mid X]$（對 $Z$ 解析積分），則 $\mathbb{E}[g(X)]=\mathbb{E}[Y]$ 且 \[\operatorname{Var}\big(g(X)\big)=\operatorname{Var}(Y)-\mathbb{E}\big[\operatorname{Var}(Y\mid X)\big]\le\operatorname{Var}(Y).\] 這是 Rao–Blackwell 定理：條件期望永遠不會增加變異數。Gibbs 抽樣裡「用條件分布的解析平均取代樣本平均」就是這招。

綜合比較

=== Example 5: variance reduction for I = int_0^1 e^x dx = e - 1 ===
exact                 = 1.718282
naive                 = 1.717677  s.e. 0.001101
antithetic            = 1.718582  s.e. 0.000198   (variance / naive = 0.0323)
control variate (c=1.689) = 1.718576  s.e. 0.000140   (variance / naive = 0.0163)
stratified (K=100)    = 1.718268  s.e. 0.000012   (variance / naive = 1.10e-04)

技術	需要知道什麼	適用場景	風險
對偶變數	無	$f$ 單調、接近線性	$f$ 非單調時反效果
控制變數	相關變數的 $\mathbb{E}[X]$	有「近似問題」有封閉解	估 $c$ 引入小偏差
分層抽樣	各層機率 $p_k$	低維、可切層	高維層數爆炸
重要性抽樣	可算 $p/q$	罕見事件、尾機率、高維	權重退化、無限變異數
條件 MC	部分解析積分	有可積的子結構	需要數學推導

馬可夫鏈蒙地卡羅（MCMC）

你想知道一座城市各區的人口密度，但沒有地圖也不能俯瞰。你的辦法：隨便選一條街站好，每一步隨機挑一個相鄰街區，如果那裡人比較多就一定過去，比較少就以「人數比」的機率過去。走得夠久之後，你出現在每個街區的頻率，會恰好正比於那裡的人口密度。你從來不需要知道全市總人口——只需要相鄰兩區的人口比。這就是 Metropolis 演算法。

動機：為什麼前面的方法不夠

貝氏統計的後驗分布 $\pi(\theta)\propto p(\text{data}\mid\theta)\,p(\theta)$、統計物理的 Boltzmann 分布 $\pi(x)\propto e^{-E(x)/T}$，有兩個共同特徵：

只知道到常數倍：正規化常數 $Z=\int\tilde\pi(x)\,dx$ 本身就是高維積分，算不出來。
高維：$\theta$ 可能有幾百維、Ising 模型有 $L^2$ 個自旋。逆變換不可能、拒絕法接受率 $\approx0$。

MCMC 的答案：建一條馬可夫鏈，讓它的平穩分布恰好是 $\pi$；讓鏈走很久，把軌跡當樣本。整個過程只用到 $\pi(y)/\pi(x)$——正規化常數約掉了。

馬可夫鏈基礎

定義 6.1. 狀態空間 $\mathcal X$ 上的馬可夫鏈由轉移核 $K(x,y)=\mathbb{P}(X_{t+1}=y\mid X_t=x)$ 決定。分布 $\pi$ 稱為平穩分布若 \[\sum_x\pi(x)K(x,y)=\pi(y)\quad\forall y\qquad\text{（連續情形把和換成積分）}.\]

定理 6.1 (遍歷定理). 若馬可夫鏈不可約（任兩狀態互通）、非週期、且有平穩分布 $\pi$，則對任何起點，$X_t$ 的分布收斂到 $\pi$，並且對任何可積 $f$， \[\frac1T\sum_{t=1}^T f(X_t)\ \xrightarrow[T\to\infty]{\text{a.s.}}\ \mathbb{E}_\pi[f].\]

這是「相依樣本版的大數法則」。代價：樣本相關，有效樣本數比 $T$ 少（第 6.5 節）。

定理 6.2 (細緻平衡 $\Rightarrow$ 平穩). 若 $\pi(x)K(x,y)=\pi(y)K(y,x)$ 對所有 $x,y$ 成立（細緻平衡 / detailed balance），則 $\pi$ 是 $K$ 的平穩分布。

Proof. $\sum_x\pi(x)K(x,y)=\sum_x\pi(y)K(y,x)=\pi(y)\sum_xK(y,x)=\pi(y)$。 ◻

細緻平衡是「$x\to y$ 的機率流等於 $y\to x$ 的機率流」，比平穩性強，但容易驗證——MCMC 演算法設計的核心就是構造滿足細緻平衡的轉移核。

Metropolis–Hastings 演算法

目標 $\tilde\pi(x)\propto\pi(x)$（不需正規化）、提議分布 $q(y\mid x)$、起點 $x_0$、步數 $T$ 抽提議 $y\sim q(\cdot\mid x_t)$ 計算接受率 $\displaystyle\alpha(x_t,y)=\min\Big\{1,\ \frac{\tilde\pi(y)\,q(x_t\mid y)}{\tilde\pi(x_t)\,q(y\mid x_t)}\Big\}$ 抽 $u\sim U(0,1)$ $x_{t+1}\gets y$ $x_{t+1}\gets x_t$ 軌跡 $x_1,\dots,x_T$（丟掉前段 burn-in）

若 $q$ 對稱（$q(y\mid x)=q(x\mid y)$，例如隨機漫步 $y=x+\sigma\varepsilon$），接受率簡化為 $\min\{1,\pi(y)/\pi(x)\}$——這是 1953 年的原始 Metropolis 演算法：往上一定走，往下以比例機率走。

定理 6.3. Metropolis–Hastings 的轉移核滿足對 $\pi$ 的細緻平衡，故 $\pi$ 為其平穩分布。

Proof. 對 $x\ne y$，$K(x,y)=q(y\mid x)\alpha(x,y)$。 \[\pi(x)K(x,y)=\pi(x)q(y\mid x)\min\Big\{1,\frac{\pi(y)q(x\mid y)}{\pi(x)q(y\mid x)}\Big\} =\min\big\{\pi(x)q(y\mid x),\ \pi(y)q(x\mid y)\big\},\] 右邊對 $x,y$ 對稱，故等於 $\pi(y)K(y,x)$。$x=y$ 時細緻平衡自動成立。由定理 6.2 得證。 ◻

注意證明中只用到 $\pi$ 的比值，因此 $\tilde\pi$ 與 $\pi$ 差常數倍不影響。

步長的藝術：實驗

用隨機漫步 Metropolis 抽 $\mathcal N(0,1)$，提議 $y=x+\sigma\varepsilon$，$\varepsilon\sim\mathcal N(0,1)$，起點故意設在 $x_0=5$：

=== Example 7: random-walk Metropolis for N(0,1): step size vs acceptance vs ESS ===
sigma   accept    mean      var       tau       ESS
0.1     0.968     -0.087    0.973     469.0     383
1.0     0.705     0.003     1.006     7.9       22797
2.4     0.442     -0.005    0.999     4.2       42398
10.0    0.126     0.007     1.014     10.8      16597
50.0    0.025     0.022     0.941     59.7      3016

$\sigma=0.1$：接受率 97%，看似很好，其實每步只挪一點點，自相關時間 $\tau=469$，18 萬個樣本只值 383 個獨立樣本，平均值還偏著（$-0.087$，尚未忘記起點）。
$\sigma=50$：幾乎都被拒絕（2.5%），鏈一直停在原地，同樣糟。
$\sigma=2.4$：接受率 44%，$\tau=4.2$，ESS 最高。理論最佳值：一維常態目標的最佳步長約 $2.38\,\sigma_\pi$、接受率 $\approx0.44$（Roberts & Rosenthal）；高維極限為 $2.38/\sqrt d$、接受率 $\approx0.234$。實驗與理論完全吻合。

接受率目標：一維 $\approx 0.44$，高維 $\approx 0.23$。太高 $\Rightarrow$ 步長加大；太低 $\Rightarrow$ 步長縮小。更好的辦法是自適應 MCMC：用已有樣本的共變異數矩陣當提議分布的形狀（Haario et al. 2001）。

收斂診斷

Burn-in：丟掉鏈「還記得起點」的前段。沒有公式，看軌跡圖（trace plot）或多條鏈是否會合。
積分自相關時間 $\tau=1+2\sum_{k\ge1}\rho_k$，$\rho_k$ 為 lag-$k$ 自相關。相依樣本的平均值變異數為 $\sigma^2\tau/N$，而非 $\sigma^2/N$。
有效樣本數 $\mathrm{ESS}=N/\tau$：真正「獨立等價」的樣本數，決定標準誤差 $\sigma/\sqrt{\mathrm{ESS}}$。
Gelman–Rubin $\hat R$：從分散起點跑多條鏈，比較「鏈間變異」與「鏈內變異」，$\hat R\approx1$ 才算收斂。這是實務界的標準檢查。

「鏈看起來已經穩定」$\ne$「鏈已經收斂到 $\pi$」。多峰分布中，鏈可能在一個峰待上一億步、完全沒發現另一個峰的存在——軌跡圖完美、$\hat R=1$（如果所有鏈都困在同一個峰）、答案全錯。對付多峰：平行回火（parallel tempering）、多起點、或改用 SMC。

Gibbs 抽樣

當聯合分布 $\pi(x_1,\dots,x_d)$ 難抽、但每個條件分布 $\pi(x_j\mid x_{-j})$ 好抽時：

抽 $x_j^{(t+1)}\sim\pi\big(x_j\ \big|\ x_1^{(t+1)},\dots,x_{j-1}^{(t+1)},x_{j+1}^{(t)},\dots,x_d^{(t)}\big)$

Gibbs 是 MH 的特例：以「完整條件分布」為提議時，接受率恆為 1（讀者可代入 MH 公式驗證）。

範例 6.1 (二維相關常態). $(X,Y)\sim\mathcal N(0,\Sigma)$，$\Sigma=\begin{pmatrix}1&\rho\\ \rho&1\end{pmatrix}$。條件分布 $X\mid Y=y\sim\mathcal N(\rho y,1-\rho^2)$，對稱地 $Y\mid X=x\sim\mathcal N(\rho x,1-\rho^2)$。

=== Example 9: Gibbs sampler for bivariate normal with rho = 0.8 ===
mean = (0.0133, 0.0127)   var = (1.0018, 1.0019)   corr = 0.8004  (target 0.8)

$\rho\to1$ 時 Gibbs 會變得極慢（每次只能沿座標軸走，而分布是一條斜的細長橢圓）——這時要重新參數化或改用 HMC。

應用：貝氏推論

範例 6.2 (硬幣的偏差). 擲 12 次得 9 次正面，先驗 $\theta\sim\mathrm{Beta}(2,2)$（「大概公平，但不確定」）。後驗 \[\pi(\theta\mid\text{data})\propto\theta^{9}(1-\theta)^{3}\cdot\theta^{1}(1-\theta)^{1}=\theta^{10}(1-\theta)^{4} \quad\Rightarrow\quad\theta\mid\text{data}\sim\mathrm{Beta}(11,5).\] 這例子刻意選有封閉解的共軛模型，才能驗證 MH 是否正確。程式只寫下對數後驗 $10\ln\theta+4\ln(1-\theta)$，完全不需要正規化常數 $B(11,5)$：

=== Example 8: Bayesian coin (9 heads / 12 tosses, Beta(2,2) prior) via Metropolis-Hastings ===
acceptance rate      = 0.5410
posterior mean  (MH) = 0.6879   exact Beta(11,5) mean = 0.6875
posterior var   (MH) = 0.0126   exact = 0.0126
95% credible interval (MH) = [0.4498, 0.8819]
P(theta > 0.5 | data)      = 0.9409

最後一行是貝氏推論的招牌能力：「硬幣偏向正面的機率是 94%」——這種敘述頻率學派無法給出，而在 MCMC 樣本上只是數一數 $\theta_i>0.5$ 的比例。

應用：統計物理的 Ising 模型

$L\times L$ 格子上每點自旋 $s_i\in\{-1,+1\}$，能量 $E(s)=-J\sum_{\langle i,j\rangle}s_is_j$（相鄰同向能量低），Boltzmann 分布 $\pi(s)\propto e^{-E(s)/T}$（取 $J=k_B=1$）。狀態數 $2^{L^2}$——$L=24$ 時 $2^{576}$，正規化常數 $Z$ 連神都算不出，但 Metropolis 只需要翻轉一個自旋的能量差 \[\Delta E=2\,s_i\sum_{j\in\mathrm{nb}(i)}s_j,\qquad \alpha=\min\{1,e^{-\Delta E/T}\}.\]

Onsager 於 1944 年解析證明二維 Ising 在 $T_c=2/\ln(1+\sqrt2)\approx2.269$ 有相變：低溫時自發磁化 $|m|>0$，高溫時 $m=0$。

=== Example 10: 2D Ising model (L=24), |magnetization| vs temperature, T_c ~ 2.269 ===
T       <|m|>     <E>/N
1.50    0.9865    -1.9512
2.00    0.9123    -1.7477
2.27    0.6990    -1.4489
2.60    0.2144    -1.0262
3.50    0.0816    -0.6599

模擬與 Onsager 解在低溫段吻合；$T>T_c$ 時 $\langle|m|\rangle$ 不為零是有限尺寸效應（$L=24$ 的系統漲落 $\sim1/L$），這正是計算物理裡「有限尺寸標度分析」的起點。另一個值得注意的現象：靠近 $T_c$ 時 Metropolis 的自相關時間暴增（臨界慢化），實務改用 Wolff / Swendsen–Wang 叢集演算法。

進階方向（速覽）

Hamiltonian Monte Carlo（HMC）：借用物理的哈密頓動力學，用梯度資訊提議「遠又高接受率」的點，是 Stan、PyMC 的引擎。
平行回火（Parallel Tempering）：同時跑多個溫度的鏈並交換狀態，高溫鏈負責跨峰、低溫鏈負責精細，治多峰。
序列蒙地卡羅 / 粒子濾波（SMC）：時間序列的線上貝氏推論（機器人定位、目標追蹤），本質是重要性抽樣 + 重抽樣。
模擬退火：把 MH 的溫度隨時間降到零，鏈就從「抽樣」變成「最佳化」——見第 8.1 節。

Monte Carlo 型隨機化演算法

本章切換到演算法理論的視角：用隨機性交換效率，容忍一個可控的錯誤機率。

錯誤類型與錯誤放大

定義 7.1. 決策問題的隨機化演算法 $A$：

單邊錯誤（one-sided）：答「否」一定對；答「是」可能錯，機率 $\le p<1$。（複雜度類 RP / co-RP）
雙邊錯誤（two-sided）：兩種答案都可能錯，錯誤機率 $\le p<1/2$。（複雜度類 BPP）
零錯誤（Las Vegas）：永不出錯，期望時間多項式。（複雜度類 ZPP）

定理 7.1 (錯誤放大). 單邊錯誤演算法獨立執行 $k$ 次、任一次答「否」即輸出「否」，錯誤機率降為 $p^k$。雙邊錯誤演算法執行 $k$ 次取多數決，錯誤機率 $\le e^{-2k(1/2-p)^2}$（Hoeffding）。

這就是「Monte Carlo 演算法」的哲學：單次可能錯，但錯誤機率隨重複次數指數下降。重複 50 次錯誤率 $4^{-50}\approx10^{-30}$，比硬體宇宙射線翻位元的機率還低。

Miller–Rabin 質數測試

為什麼 Fermat 測試不夠

Fermat 小定理：$n$ 為質數 $\Rightarrow a^{n-1}\equiv1\pmod n$ 對所有 $\gcd(a,n)=1$。逆命題不真：Carmichael 數（如 $561=3\cdot11\cdot17$）對所有互質底數都通過 Fermat 測試——Fermat 測試對它們的錯誤率不是 $\le1/2$，而是 100%。

強偽質數與演算法

定義 7.2 (強可能質數). 奇數 $n>2$，寫 $n-1=2^r d$（$d$ 奇）。稱 $a$ 為 $n$ 的強說謊者（或 $n$ 對底 $a$ 是強可能質數）若 \[a^d\equiv1\pmod n\quad\text{或}\quad a^{2^i d}\equiv-1\pmod n\ \text{對某個 }0\le i<r.\] 否則 $a$ 是 $n$ 為合成數的見證者（witness）。

質數一定讓所有 $a$ 成為強說謊者（因為 $\mathbb Z_p$ 中 $x^2\equiv1$ 只有 $\pm1$ 兩個根，反覆開根號只能停在 $1$ 或 $-1$）。關鍵是反方向：

定理 7.2 (Rabin 1980). 若 $n$ 是奇合成數，則 $\{1,\dots,n-1\}$ 中強說謊者至多 $\dfrac{n-1}{4}$ 個。

因此隨機挑 $a$，合成數被誤判為「可能是質數」的機率 $\le1/4$；$k$ 輪 $\le4^{-k}$。這是單邊錯誤：說「合成數」一定對（找到見證者是鐵證）。

=== Example 11: Miller-Rabin primality test (k = 20 random bases) ===
561                  -> composite        [Carmichael number (fools Fermat test for all coprime bases)]
1105                 -> composite        [Carmichael number]
1000000007           -> probably prime   [prime (1e9+7)]
2147483647           -> probably prime   [prime (2^31-1, Mersenne)]
3215031751           -> composite        [composite = 151*751*28351 (strong pseudoprime to bases 2,3,5,7)]
2305843009213693951  -> probably prime   [prime (2^61-1, Mersenne)]
1000000016000000063  -> composite        [composite = (1e9+7)*(1e9+9)]
single-base error rate on 561 = 0.0142  (theoretical bound 0.25; after k=20 rounds <= 4^-20 ~ 9.1e-13)

$3215031751$ 是能同時騙過底數 $2,3,5,7$ 的最小合成數——固定底數的測試會被精心構造的數欺騙，隨機底數則不會（對手不知道你會抽什麼）。這是隨機化的另一個好處：對抗性輸入免疫。實測 561 的單底錯誤率 1.4%，遠低於 $1/4$ 上界（上界在 $n=pq$ 且 $p,q$ 特殊時才逼近）。

64 位元以內用固定底數 $\{2,3,5,7,11,13,17,19,23,29,31,37\}$ 即為確定性演算法（已驗證無反例）。RSA 金鑰生成用的 2048 位元質數則靠隨機底數的 Miller–Rabin（OpenSSL 跑約 64 輪）。AKS（2002）是多項式時間確定性演算法，但實務上比 Miller–Rabin 慢幾個數量級。

Freivalds 演算法：驗證矩陣乘法

問題：給 $n\times n$ 矩陣 $A,B,C$，判斷 $AB=C$？直接乘要 $O(n^3)$（或 $O(n^{2.37})$）。

隨機抽 $r\in\{0,1\}^n$ 「$AB\ne C$」「$AB=C$（機率 $\ge1-2^{-k}$）」

每輪只做三次矩陣–向量乘法，$O(n^2)$。

引理 7.1. 若 $D=AB-C\ne0$，則 $\mathbb{P}_r(Dr=0)\le1/2$。

Proof. $D$ 有某非零元 $d_{ij^*}$。$(Dr)_i=\sum_j d_{ij}r_j=d_{ij^*}r_{j^*}+\sum_{j\ne j^*}d_{ij}r_j$。固定 $r_{j}$（$j\ne j^*$）之後，令 $(Dr)_i=0$ 的 $r_{j^*}$ 值最多一個（一元一次方程），而 $r_{j^*}\in\{0,1\}$ 兩種等機率，故 $\mathbb{P}((Dr)_i=0)\le1/2$，於是 $\mathbb{P}(Dr=0)\le1/2$。 ◻

=== Example 12: Freivalds' matrix product verification (n = 200) ===
correct C, k=10  -> accepted
corrupted C (one entry +1), k=10 -> rejected
single-round false-accept rate on corrupted C = 0.5155  (bound 0.5)

只改了 $C$ 的一個元素 $+1$，單輪偵測率約 50%（正是最壞情況——只有一列有非零元），10 輪偵測率 $1-2^{-10}=99.9\%$。若把 $r$ 改抽自 $\{0,\dots,q-1\}$，單輪錯誤率降到 $1/q$。

Karger 隨機收縮最小割

無向圖 $G$ 的最小割：把頂點分成兩非空集合，跨越的邊數最少。確定性演算法（Stoer–Wagner）$O(nm+n^2\log n)$。Karger（1993）的演算法簡單到不可思議：

均勻隨機挑一條邊 $(u,v)$，把 $u,v$ 合併成一個頂點（保留多重邊、刪除自環）剩下兩個超頂點之間的邊數

定理 7.3. 單次收縮回傳最小割的機率 $\ge\dfrac{2}{n(n-1)}=\dfrac{1}{\binom n2}$。

Proof. 固定一個最小割 $C$，$|C|=k$。收縮成功 $\iff$ 從頭到尾都沒挑到 $C$ 中的邊。每個頂點的度數 $\ge k$（否則該頂點單獨一邊就是更小的割），故當前有 $i$ 個頂點時 $|E_i|\ge ik/2$，挑到 $C$ 中邊的機率 $\le k/(ik/2)=2/i$。 \[\mathbb{P}(\text{成功})\ \ge\ \prod_{i=n}^{3}\Big(1-\frac2i\Big)=\prod_{i=3}^{n}\frac{i-2}{i} =\frac{1\cdot2\cdot3\cdots(n-2)}{3\cdot4\cdots n}=\frac{2}{n(n-1)}.\qedhere\] ◻

推論 7.1. 重複 $T=\binom n2\ln n$ 次取最小值，失敗機率 $\le\big(1-\tfrac{1}{\binom n2}\big)^T\le e^{-\ln n}=\tfrac1n$。

=== Example 13: Karger's randomized min-cut (two K5 cliques joined by 2 edges) ===
n = 10, |E| = 22, repetitions T = 104
best cut found = 2 (true min cut = 2)
empirical single-run success rate = 0.3077  (lower bound 2/(n(n-1)) = 0.0222)

實測單次成功率 31%，遠高於下界 2.2%——下界是對「最壞的圖」證的，這個圖的最小割邊只佔 22 條中的 2 條，很難被誤挑。Karger–Stein（1996）用「收縮到 $n/\sqrt2$ 就分岔遞迴」把總時間降到 $O(n^2\log^3n)$，仍是隨機化演算法教科書的明星案例。

對照：Las Vegas 演算法

隨機化快速排序（隨機挑 pivot）永遠排對，只是時間隨機——期望 $O(n\log n)$，且對任何輸入都成立（沒有「最壞輸入」，只有「最壞運氣」）。這與本章的 Monte Carlo 型演算法互補：前者用隨機性消滅對抗性輸入，後者用隨機性交換速度並接受可控錯誤。

蒙地卡羅最佳化：模擬退火與 MCTS

模擬退火（Simulated Annealing）

鐵匠鍛劍：把鐵燒紅（高溫，原子亂跳、什麼結構都能形成），然後慢慢冷卻（原子逐漸安定到能量最低的晶格）。急速冷卻會凍出雜亂的結構（局部最優）。Kirkpatrick 等人在 1983 年把這個物理過程變成了最佳化演算法。

把目標函數 $E(x)$ 當能量、目標分布 $\pi_T(x)\propto e^{-E(x)/T}$。$T$ 高時 $\pi_T$ 近乎均勻（到處亂走）；$T\to0$ 時 $\pi_T$ 集中在全域最小值。模擬退火 $=$ Metropolis 演算法 + 逐漸降溫：

初始解 $x$、初溫 $T_0$、降溫率 $\alpha<1$、終溫 $T_{\min}$ $T\gets T_0$ 產生鄰居 $y$（例如 TSP 的 2-opt 反轉） $\Delta\gets E(y)-E(x)$ $x\gets y$ $T\gets\alpha T$

定理 8.1 (Geman & Geman 1984，收斂性). 若降溫 $T_k\ge c/\ln(k+1)$（$c$ 足夠大，與能量地形的最深「谷深」有關），則 $x_k$ 以機率 1 收斂到全域最小值集合。

對數降溫慢到毫無實用價值（要達 $T=0.01$ 需 $e^{100c}$ 步），但定理說明了精神：降得太快就會被凍在局部最優。實務用幾何降溫 $T\gets\alpha T$，$\alpha\in[0.9,0.999]$。

=== Example 14: simulated annealing on a 15-city TSP vs Held-Karp optimum ===
cooling alpha = 0.900  steps = 10300     SA best = 324.57  OPT = 324.57  gap = +0.00%
cooling alpha = 0.990  steps = 107700    SA best = 324.57  OPT = 324.57  gap = +0.00%
cooling alpha = 0.999  steps = 1081500   SA best = 324.57  OPT = 324.57  gap = +0.00%

15 城對 SA 太簡單，三種降溫速度都命中 Held–Karp 最優解。真正的差別在 $n$ 上百時：$\alpha=0.9$ 常卡在離最優 5–10% 的局部解，$\alpha=0.999$ 則能逼到 1–2%。與遺傳演算法比較：SA 是單點搜尋、靠溫度逃離局部最優；GA 是族群搜尋、靠交配重組。兩者都是「有偏好的隨機探索」。

先隨機走幾百步、記錄 $\Delta$ 的典型大小 $\bar\Delta$，令 $T_0$ 使初始接受率約 80%：$T_0=-\bar\Delta/\ln0.8\approx4.5\bar\Delta$。終溫則設為「接受任何變差都幾乎不可能」：$T_{\min}\approx\Delta_{\min}/\ln 1000$。

蒙地卡羅樹搜尋（MCTS）

問題：博弈樹太大

井字棋的博弈樹約 $5.5\times10^5$ 節點，可以窮舉；圍棋約 $10^{170}$，不可能。傳統 Minimax + $\alpha$–$\beta$ 剪枝需要一個好的靜態評估函數，而圍棋很難寫出來。MCTS 的想法：評估一個局面，就從那裡隨機下到終局幾千次，看勝率。不需要任何領域知識。

四個階段

UCT：探索與利用的平衡

選擇階段要決定往哪個子節點走。只走勝率最高的（純利用）會錯過還沒試夠的好棋；均勻亂走（純探索）浪費模擬。Kocsis & Szepesvári（2006）把多臂吃角子老虎機的 UCB1 公式搬進樹裡：

\[\text{選擇 }\ j^*=\mathop{\mathrm{arg\,max}}_j\ \Bigg[\ \underbrace{\frac{w_j}{n_j}}_{\text{利用：勝率}}\ +\ C\underbrace{\sqrt{\frac{\ln N}{n_j}}}_{\text{探索：不確定性}}\ \Bigg],\] $w_j$、$n_j$ 為子節點 $j$ 的勝場數與訪問數，$N=\sum_j n_j$ 為父節點訪問數，$C=\sqrt2$ 為理論值（勝率在 $[0,1]$ 時）。

第二項來自 Hoeffding 不等式——它是勝率的信賴區間上界。試得少的節點區間寬、上界高、有機會被選到；試得多的節點區間收窄、只靠真實勝率競爭。UCB1 的累積遺憾（regret）為 $O(\ln N)$，是理論最優階。

=== Example 15: Monte Carlo Tree Search (UCT) playing tic-tac-toe ===
MCTS opening move from empty board (20000 iters): cell 4 (center)
MCTS(100)  vs Random     : A wins 186, B wins 3, draws 11  (200 games)
MCTS(500)  vs Random     : A wins 189, B wins 0, draws 11  (200 games)
MCTS(3000) vs Random     : A wins 187, B wins 0, draws 13  (200 games)
MCTS(3000) vs MCTS(3000): A wins 0, B wins 0, draws 50  (50 games)

沒有任何井字棋知識，MCTS 自己「發現」開局要佔中央。
100 次迭代偶爾輸給隨機（3/200），500 次以上零敗——迭代數就是棋力。
兩個 MCTS 對弈 50 局全和——井字棋雙方完美下法必和，MCTS 已達完美。

終局節點不可再擴展。棋局已分勝負的節點若仍列出「空格」當作可走步，MCTS 會在勝負已定後繼續「下棋」，讓兩方都可能連成線，汙染整棵樹的估值。本教材配套程式開發時就踩過這個坑：修正前 3000 次迭代仍會輸隨機玩家 5%，修正後零敗。

AlphaGo（2016）= MCTS + 兩個神經網路：策略網路取代均勻隨機挑選擴展的子節點（先驗機率進 UCT 公式，即 PUCT），價值網路取代隨機模擬到終局。AlphaZero 更進一步用自我對弈的 MCTS 結果當訓練標籤。骨架仍是本節的四個階段。

金融與隨機系統模擬

選擇權定價

模型

風險中性測度下股價服從幾何布朗運動（GBM）： \[dS_t=rS_t\,dt+\sigma S_t\,dW_t \quad\Rightarrow\quad S_T=S_0\exp\!\Big(\big(r-\tfrac{\sigma^2}{2}\big)T+\sigma\sqrt T\,Z\Big),\ Z\sim\mathcal N(0,1).\] 歐式買權（到期才能執行）的價格是折現後的期望報酬 \[C=e^{-rT}\,\mathbb{E}\big[\max(S_T-K,0)\big].\] 這正是一個「期望值」——蒙地卡羅的主場。此例有 Black–Scholes 封閉解可對照： \[C_{\mathrm{BS}}=S_0\Phi(d_1)-Ke^{-rT}\Phi(d_2),\quad d_1=\frac{\ln(S_0/K)+(r+\sigma^2/2)T}{\sigma\sqrt T},\ d_2=d_1-\sigma\sqrt T.\]

=== Example 16: European call by Monte Carlo vs Black-Scholes; Asian call by path simulation ===
Black-Scholes price = 10.4506
N = 10000     plain = 10.6140 (s.e. 0.1483)   antithetic = 10.2882 (s.e. 0.0989)
N = 100000    plain = 10.4413 (s.e. 0.0466)   antithetic = 10.4786 (s.e. 0.0331)
N = 1000000   plain = 10.4623 (s.e. 0.0147)   antithetic = 10.4514 (s.e. 0.0104)
Asian (arithmetic-average, 52 steps) call = 5.8682 (s.e. 0.0182)  [no closed form]

$S_0=K=100$、$r=5\%$、$\sigma=20\%$、$T=1$ 年，Black–Scholes 價格 10.4506。蒙地卡羅在 $N=10^6$ 時誤差 0.01、標準誤差 0.015，一致。對偶變數（$Z$ 與 $-Z$ 配對）把變異數再砍一半——報酬函數 $\max(S-K,0)$ 單調於 $Z$，符合對偶變數的適用條件。

為什麼金融業離不開蒙地卡羅

亞式選擇權的報酬取決於整條路徑的平均價 $\frac1n\sum S_{t_i}$，算術平均的對數常態和沒有封閉解；籃子選擇權依賴多個資產（高維積分）；路徑相依的障礙、回顧選擇權；隨機波動率模型（Heston）。這些都是「期望值」，都能用同一套路徑模擬 + 報酬平均解決，只是 $Y_i$ 的定義不同。變異數縮減（控制變數 $\to$ 幾何平均亞式、重要性抽樣 $\to$ 深度價外選擇權、準蒙地卡羅 $\to$ Sobol 序列）在此是產業標準配備。

隨機漫步與賭徒破產

賭徒本金 $i$，每局以機率 $p$ 贏 1 元、否則輸 1 元，達到 $N$ 元或破產即停。

定理 9.1 (破產機率). 令 $q=1-p$，$P_i=\mathbb{P}(\text{由 }i\text{ 出發破產})$。則 \[P_i=\begin{cases}1-\dfrac{i}{N}, & p=\tfrac12,\\[8pt] \dfrac{(q/p)^i-(q/p)^N}{1-(q/p)^N}, & p\ne\tfrac12.\end{cases} \qquad \mathbb{E}[\text{局數}]=i(N-i)\ \ (p=\tfrac12).\]

證明概要. 首步分析（first-step analysis）給出遞迴 $P_i=pP_{i+1}+qP_{i-1}$，邊界 $P_0=1,P_N=0$。這是二階線性差分方程，特徵根 $1$ 與 $q/p$（$p\ne\frac12$）或重根 $1$（$p=\frac12$），套邊界條件即得。 ◻

=== Example 17: gambler's ruin, start i = 10, target N = 50, 100000 simulated games ===
p = 0.50  P(ruin) MC = 0.8011  exact = 0.8000  mean game length = 401.2
p = 0.49  P(ruin) MC = 0.9233  exact = 0.9230  mean game length = 308.8
p = 0.45  P(ruin) MC = 0.9998  exact = 0.9997  mean game length = 100.1

$p$ 從 0.50 降到 0.49（賭場的 1% 優勢），破產率從 80% 跳到 92%；到 0.45（輪盤大約如此）幾乎必破產。模擬的價值在於：模型一旦複雜到沒有封閉解（變動賭注、多種賭局、停損策略），程式碼幾乎不用改。

其他重要應用領域

排隊與作業系統：呼叫中心人力、資料中心排程、交通號誌——以指數/一般分布抽事件間隔，統計等待時間分布。
可靠度工程：系統失效機率 $10^{-6}$ 級別的估計，重要性抽樣與子集模擬（subset simulation）是標準工具。
粒子傳輸：中子、光子在物質中的散射與吸收（MCNP、Geant4）——蒙地卡羅的誕生地，至今仍是核工與放射治療劑量計算的黃金標準。
電腦圖學：路徑追蹤（path tracing）用蒙地卡羅解渲染方程的高維積分，好萊塢電影與現代遊戲的光影都靠它；重要性抽樣 BRDF、多重重要性抽樣（MIS）是核心技術。
氣象與氣候：集合預報（ensemble forecasting）對初始條件加擾動、跑幾十次模型，「降雨機率 70%」就是這麼來的。
機器學習：隨機梯度下降（SGD）本質是梯度的蒙地卡羅估計；變分自編碼器的重參數化技巧、強化學習的策略梯度、Dropout 的貝氏解釋、擴散模型的抽樣——蒙地卡羅無所不在。
統計推論：Bootstrap（對資料重抽樣估計統計量的分布）、置換檢定——不需任何分布假設的推論工具。

完整 C++ 程式

以下為配套程式 monte_carlo.cpp 全文，包含本教材全部 17 個實例。單一檔案、僅用標準函式庫、固定隨機種子、每個實例皆與理論值或精確演算法對照。編譯與執行：

g++ -std=c++17 -O2 -Wall -Wextra -o monte_carlo monte_carlo.cpp
./monte_carlo          # 全部跑完約 3 秒

程式結構一覽：

實例	函數	主題	對照基準
1	`ex1_pi`	hit-or-miss 估 $\pi$、收斂率	$\pi$
2	`ex2_buffon`	Buffon 投針	$\pi$
3	`ex3_two_estimators`	hit-or-miss vs. 平均值法	$\frac{\sqrt\pi}{2}\operatorname{erf}(1)$
4	`ex4_ball_volume`	高維球體積、維度災難	$\pi^{d/2}/\Gamma(d/2+1)$
5	`ex5_variance_reduction`	對偶 / 控制 / 分層	$e-1$
6	`ex6_importance_sampling`	重要性抽樣、尾機率	$1-\Phi(4)$
7	`ex7_metropolis_normal`	隨機漫步 Metropolis、步長、ESS	$\mathcal N(0,1)$
8	`ex8_bayesian_coin`	MH 抽貝氏後驗	$\mathrm{Beta}(11,5)$
9	`ex9_gibbs`	Gibbs 抽樣	$\rho=0.8$
10	`ex10_ising`	Ising 模型、相變	Onsager $T_c$
11	`ex11_miller_rabin`	隨機化質數測試	已知質數/合成數
12	`ex12_freivalds`	矩陣乘法驗證	直接相乘
13	`ex13_karger`	隨機收縮最小割	構造圖的已知最小割
14	`ex14_simulated_annealing`	模擬退火 TSP	Held–Karp DP
15	`ex15_mcts`	MCTS / UCT 井字棋	完美下法必和
16	`ex16_option_pricing`	GBM 路徑、歐式/亞式選擇權	Black–Scholes
17	`ex17_gamblers_ruin`	隨機漫步、破產機率	差分方程封閉解

完整執行結果

以下為 ./monte_carlo 一次執行的完整輸出（Apple clang 21，-O2，約 3 秒），各段解讀已分散在前面各章。

所有隨機性來自單一 std::mt19937_64 引擎與固定種子；std::normal_distribution 的實作在不同標準函式庫（libc++ / libstdc++）間可能不同，因此換平台後數字可能略異但統計結論不變——這本身就是蒙地卡羅方法的重要一課：結論要建立在信賴區間上，而不是某一次的數字上。

練習題

（基礎）用 hit-or-miss 法估計 $\int_0^1\sqrt{1-x^2}\,dx=\pi/4$，並用 Hoeffding 不等式算出：以 99% 信心讓誤差小於 $10^{-3}$ 所需的 $N$。再以 CLT 估一次，比較兩者。
（估計量）證明 Buffon 投針的 $\hat\pi=1/\hat p$ 有正偏差（$\mathbb{E}[\hat\pi]>\pi$），並以 delta method（把 $g(p)=1/p$ 在 $p_0=1/\pi$ 二階展開）導出偏差的一階近似 $\mathbb{E}[\hat\pi]-\pi\approx\dfrac{g''(p_0)}{2}\cdot\dfrac{p_0(1-p_0)}{N}=\dfrac{1-p_0}{p_0^2N}=\dfrac{\pi(\pi-1)}{N}$，再用模擬（$N=10^2,10^3,10^4$，各重複一萬次）驗證偏差確為 $O(1/N)$。
（抽樣）用逆變換法實作 Cauchy 分布 $F^{-1}(u)=\tan(\pi(u-\frac12))$ 的抽樣；對其樣本平均畫「$N$ 對 $\bar X_N$」的圖，解釋為什麼它不收斂（大數法則的哪個假設被違反？）。
（拒絕法）以 $q=\mathcal N(0,1)$ 為包絡，用拒絕抽樣生成截斷常態 $\mathcal N(0,1)\mid X>2$。計算理論接受率，再解釋為什麼截斷點改成 $X>6$ 時這方法徹底失敗，並提出替代方案（提示：指數分布包絡或逆變換）。
（變異數縮減）對 $f(x)=(x-\frac12)^2$ 在 $[0,1]$ 上比較樸素法與對偶變數法的變異數。解釋結果為何與 $e^x$ 的例子相反。
（重要性抽樣）估計 $\mathbb{P}(Z>6)$。試提議分布 $\mathcal N(6,1)$ 與 $\mathcal N(6,0.5^2)$，比較 ESS 與標準誤差。哪一個有無限變異數的風險？為什麼？
（MCMC）用隨機漫步 Metropolis 抽樣雙峰分布 $\pi(x)\propto e^{-(x-5)^2/2}+e^{-(x+5)^2/2}$，步長 $\sigma=1$。畫軌跡圖，觀察鏈多久才會第一次跨峰；改 $\sigma=10$ 再觀察。這說明了收斂診斷的什麼困難？
（Gibbs）把 Gibbs 實例的 $\rho$ 改為 0.99，量測 $x$ 座標的積分自相關時間。與 $\rho=0.8$ 比較，解釋為何 Gibbs 在強相關時失效，並提出「旋轉座標」的補救。
（Ising）把 $L$ 改為 8、16、32，在 $T=2.27$ 量 $\langle|m|\rangle$。畫出 $\langle|m|\rangle$ 對 $L$ 的圖——這是有限尺寸標度分析的起點；同時量測自相關時間隨 $L$ 的成長（臨界慢化）。
（隨機化演算法）證明：若 Miller–Rabin 對 $n$ 用固定底數集合 $\{2,3,5,7\}$，則存在合成數使其失敗（$3215031751$ 即是），但用隨機底數則任何合成數的失敗率 $\le1/4$。討論「對抗性輸入」在兩種情況下的差別。
（Karger）修改 Freivalds 使 $r\in\{0,\dots,q-1\}^n$，證明單輪錯誤率 $\le1/q$。對 Karger，構造一張圖使單次成功率接近下界 $2/(n(n-1))$（提示：環狀圖 $C_n$）。
（MCTS）把 UCT 常數 $C$ 從 $\sqrt2$ 改為 $0.1$ 與 $10$，重跑 MCTS vs. Random。解釋為何兩個極端都變弱（純利用 vs. 純探索）。
（金融）用「幾何平均亞式選擇權」（有封閉解）當控制變數，重新估計算術平均亞式買權，報告變異數縮減倍數。
（挑戰）實作 Hamiltonian Monte Carlo 抽樣 100 維標準常態，與隨機漫步 Metropolis 比較 ESS/秒。這是為什麼 Stan 與 PyMC 都選 HMC 的原因。

小結與延伸閱讀

一張表看完全教材

你想做什麼	用什麼	關鍵公式 / 保證
估一個期望值 / 積分	樣本平均	誤差 $\sigma/\sqrt N$，與維度無關
知道需要多少樣本	Hoeffding	$N\ge\frac{(b-a)^2}{2\varepsilon^2}\ln\frac2\delta$
從已知 CDF 抽樣	逆變換	$X=F^{-1}(U)$
從逐點可算的密度抽樣（低維）	拒絕法	接受率 $1/M$
砍變異數（$f$ 單調）	對偶變數	$\operatorname{Cov}(f(U),f(1-U))<0$
砍變異數（有相關的已知量）	控制變數	$\operatorname{Var}\to\sigma^2(1-\rho^2)$
砍變異數（可切層）	分層抽樣	消滅層間變異數
罕見事件 / 尾機率	重要性抽樣	$q^*\propto\|f\|p$，小心權重退化
高維、未正規化的密度	MCMC（MH / Gibbs）	細緻平衡 $\Rightarrow$ 平穩；ESS $=N/\tau$
組合最佳化	模擬退火	Metropolis + 降溫
博弈決策	MCTS / UCT	$\frac{w}{n}+C\sqrt{\ln N/n}$
快速判定（容忍錯誤）	Monte Carlo 演算法	錯誤率 $p^k$ 指數下降

三個帶走的觀念

隨機不等於不精確。蒙地卡羅的誤差是可量化、可控制、可用重複次數指數壓低的；每個結果都該帶著它的標準誤差出門。
$1/\sqrt N$ 很慢，所以聰明抽樣比蠻力重要。變異數縮減與 MCMC 不是「進階選修」，是讓蒙地卡羅在真實問題上可行的必要條件。
蒙地卡羅是一種思考方式：遇到算不出來的量，先問「它是哪個隨機變數的期望值？我能抽樣嗎？」——從中子傳輸到 AlphaGo，這個問句七十年來從未失效。

維度 \(d\)	梯形法	Simpson	蒙地卡羅	勝者
1	\(N^{-2}\)	\(N^{-4}\)	\(N^{-1/2}\)	確定性方法（壓倒性）
2	\(N^{-1}\)	\(N^{-2}\)	\(N^{-1/2}\)	確定性方法
4	\(N^{-1/2}\)	\(N^{-1}\)	\(N^{-1/2}\)	梯形法打平；Simpson 略勝
8	\(N^{-1/4}\)	\(N^{-1/2}\)	\(N^{-1/2}\)	Simpson 打平
20	\(N^{-0.1}\)	\(N^{-0.2}\)	\(N^{-1/2}\)	蒙地卡羅（壓倒性）
100	\(N^{-0.02}\)	\(N^{-0.04}\)	\(N^{-1/2}\)	只有蒙地卡羅可行

技術	需要知道什麼	適用場景	風險
對偶變數	無	\(f\) 單調、接近線性	\(f\) 非單調時反效果
控制變數	相關變數的 \(\mathbb{E}[X]\)	有「近似問題」有封閉解	估 \(c\) 引入小偏差
分層抽樣	各層機率 \(p_k\)	低維、可切層	高維層數爆炸
重要性抽樣	可算 \(p/q\)	罕見事件、尾機率、高維	權重退化、無限變異數
條件 MC	部分解析積分	有可積的子結構	需要數學推導

實例	函數	主題	對照基準
1	`ex1_pi`	hit-or-miss 估 \(\pi\)、收斂率	\(\pi\)
2	`ex2_buffon`	Buffon 投針	\(\pi\)
3	`ex3_two_estimators`	hit-or-miss vs. 平均值法	\(\frac{\sqrt\pi}{2}\operatorname{erf}(1)\)
4	`ex4_ball_volume`	高維球體積、維度災難	\(\pi^{d/2}/\Gamma(d/2+1)\)
5	`ex5_variance_reduction`	對偶 / 控制 / 分層	\(e-1\)
6	`ex6_importance_sampling`	重要性抽樣、尾機率	\(1-\Phi(4)\)
7	`ex7_metropolis_normal`	隨機漫步 Metropolis、步長、ESS	\(\mathcal N(0,1)\)
8	`ex8_bayesian_coin`	MH 抽貝氏後驗	\(\mathrm{Beta}(11,5)\)
9	`ex9_gibbs`	Gibbs 抽樣	\(\rho=0.8\)
10	`ex10_ising`	Ising 模型、相變	Onsager \(T_c\)
11	`ex11_miller_rabin`	隨機化質數測試	已知質數/合成數
12	`ex12_freivalds`	矩陣乘法驗證	直接相乘
13	`ex13_karger`	隨機收縮最小割	構造圖的已知最小割
14	`ex14_simulated_annealing`	模擬退火 TSP	Held–Karp DP
15	`ex15_mcts`	MCTS / UCT 井字棋	完美下法必和
16	`ex16_option_pricing`	GBM 路徑、歐式/亞式選擇權	Black–Scholes
17	`ex17_gamblers_ruin`	隨機漫步、破產機率	差分方程封閉解

你想做什麼	用什麼	關鍵公式 / 保證
估一個期望值 / 積分	樣本平均	誤差 \(\sigma/\sqrt N\)，與維度無關
知道需要多少樣本	Hoeffding	\(N\ge\frac{(b-a)^2}{2\varepsilon^2}\ln\frac2\delta\)
從已知 CDF 抽樣	逆變換	\(X=F^{-1}(U)\)
從逐點可算的密度抽樣（低維）	拒絕法	接受率 \(1/M\)
砍變異數（\(f\) 單調）	對偶變數	\(\operatorname{Cov}(f(U),f(1-U))<0\)
砍變異數（有相關的已知量）	控制變數	\(\operatorname{Var}\to\sigma^2(1-\rho^2)\)
砍變異數（可切層）	分層抽樣	消滅層間變異數
罕見事件 / 尾機率	重要性抽樣	\(q^*\propto\|f\|p\)，小心權重退化
高維、未正規化的密度	MCMC（MH / Gibbs）	細緻平衡 \(\Rightarrow\) 平穩；ESS \(=N/\tau\)
組合最佳化	模擬退火	Metropolis + 降溫
博弈決策	MCTS / UCT	\(\frac{w}{n}+C\sqrt{\ln N/n}\)
快速判定（容忍錯誤）	Monte Carlo 演算法	錯誤率 \(p^k\) 指數下降