6.3 蒙地卡羅與迭代取樣#

概念說明#

前面所有迭代都在「逼近確定答案」，蒙地卡羅法（Monte Carlo）則換了一套哲學：把答案表述為某隨機變數的期望值，用迭代取樣 + 統計平均來估計它。估 \(\pi\) 的經典例：在單位正方形內均勻撒 \(N\) 個點，落在四分之一圓內的比例趨近 \(\pi/4\) ，故 \(\pi \approx 4 \times (\text{命中數}/N)\) 。

收斂理論由大數法則與中央極限定理給出：估計量的標準誤差為

\[ \text{SE} = \frac{\sigma}{\sqrt{N}} \]

與維度無關——這是蒙地卡羅最深刻的性質。網格法求 \(d\) 維積分需 \(m^d\) 個點（維度災難），蒙地卡羅永遠 \(O(1/\sqrt{N})\) 。代價是收斂慢：精度要 10 倍，樣本要 100 倍；故它主宰高維問題（金融衍生品定價、統計物理、貝葉斯推論），低維問題讓位給確定性方法。

「迭代」體現在三個層次：

取樣迴圈： \(N\) 次獨立取樣、累積統計量，就是一個簡單迭代。
馬可夫鏈蒙地卡羅（MCMC）：樣本不獨立，而是由馬可夫鏈迭代產生——Metropolis-Hastings、Gibbs 取樣。鏈的平穩分布即目標分布，迭代越久越「混合」，樣本越有代表性。6.2 節 PageRank 正是馬可夫鏈的平穩分布。
變異數縮減：重要性取樣、對偶變數、控制變數——不改期望值、只縮小 \(\sigma\) ，等於「更聰明地迭代」。

實務注意：MCMC 需先丟棄燒入期（burn-in）樣本、檢查自相關（混合不良的鏈即使跑很久，有效樣本數仍很少）。

生活類比#

想知道整鍋湯的鹹度不必喝完整鍋：攪勻後舀幾匙嚐（取樣平均）。攪得越勻、舀得越多，估計越準——但每多舀一匙只多一點信心（ \(1/\sqrt{N}\) ）。選舉民調同理：抽樣數千人就能估數百萬人的意向，誤差界正是 \(\sigma/\sqrt{N}\) 。

演算法步驟#

蒙地卡羅積分（估 \(\pi\) 與一般積分）：

把目標量寫成期望值： \(\theta = E[g(X)]\) ， \(X\) 可取樣。
迭代 \(N\) 次：取樣 \(x_i\) ，累積 \(s \mathrel{+}= g(x_i)\) 與 \(s_2 \mathrel{+}= g(x_i)^2\) 。
估計 \(\hat\theta = s/N\) ，標準誤 \(\text{SE} = \sqrt{(s_2/N - \hat\theta^2)/N}\) 。
報告 \(\hat\theta \pm 1.96\,\text{SE}\) （95% 信賴區間）； \(N\) 由目標精度反推： \(N \approx (1.96\sigma/\delta)^2\) 。

Metropolis-Hastings 一輪：

由當前 \(x\) 提議 \(x' = q(x' \mid x)\) 。
計算接受率 \(\min\left(1, \dfrac{\pi(x')q(x \mid x')}{\pi(x)q(x' \mid x)}\right)\) 。
以該機率接受（否則留在原地），記錄樣本。

Python 範例#

import math
import random
import numpy as np

def mc_pi(N, seed=0):
    rng = np.random.default_rng(seed)
    x, y = rng.random(N), rng.random(N)
    hits = np.sum(x*x + y*y <= 1.0)
    return 4 * hits / N

for N in [10**3, 10**4, 10**6]:
    est = mc_pi(N)
    print(f"N={N:>8}: π ≈ {est:.6f}, 誤差 {abs(est - math.pi):.2e}")

# 蒙地卡羅積分：∫₀¹ e^(-x²) dx ≈ 0.7468
rng = np.random.default_rng(1)
N = 10**6
x = rng.random(N)
vals = np.exp(-x*x)
est, se = vals.mean(), vals.std() / math.sqrt(N)
print(f"積分 ≈ {est:.5f} ± {1.96*se:.5f}（95% CI）")

# Metropolis-Hastings：從標準常態取樣
def mcmc_normal(n_samples, burn=1000, step=1.0, seed=2):
    rng = np.random.default_rng(seed)
    x, samples = 0.0, []
    for i in range(n_samples + burn):
        prop = x + rng.normal(0, step)
        if rng.random() < math.exp(-(prop*prop - x*x) / 2):
            x = prop
        if i >= burn:
            samples.append(x)
    return np.array(samples)

s = mcmc_normal(50000)
print(f"MCMC 平均 {s.mean():.3f}（目標 0）、標準差 {s.std():.3f}（目標 1）")

複雜度與收斂性分析#

成本 \(O(N \cdot c)\) （ \(c\) 為單次取樣成本），精度 \(\delta\) 需 \(N = O(\sigma^2/\delta^2)\) 樣本，維度只影響 \(\sigma\) 不影響收斂率。MCMC 的有效樣本數 = \(N / (1 + 2\sum \rho_k)\) （自相關修正），混合慢的鏈可能有效樣本只有名目的百分之一。改善收斂的實用手法：重要性取樣把樣本集中在「貢獻大」的區域；對偶變數用正負相關的樣本對（例：估 \(\pi\) 時同時用 \(x\) 與 \(1-x\) ）壓低 \(\sigma\) 。判斷何時用蒙地卡羅：維度高、期望值可表述、允許統計誤差。

練習#

實驗驗證估 \(\pi\) 誤差 \(\propto 1/\sqrt{N}\) ：對 \(N = 10^3 \sim 10^7\) 各跑 20 次，畫誤差中位數對 \(N\) 的 log-log 圖，斜率應約 \(-0.5\) 。
用對偶變數估 \(\int_0^1 e^{x}\,dx\) ：樣本對 \((e^{U}, e^{1-U})\) 平均，比較變異數下降倍數。
用 MCMC 從 \(\pi(x) \propto e^{-x^4}\) 取樣，畫直方圖與數值積分結果比較，並測試不同 step 對混合速度的影響。