6.3 蒙地卡羅積分#

概念說明#

梯形與辛普森法在低維度非常優秀，但遇到高維積分就崩潰了：d 維空間中每邊切 n 段，就要 n^d 個格點——d = 10、n = 10 就是一百億個點。這叫維度詛咒。蒙地卡羅（Monte Carlo）方法用一個反直覺的策略突圍：不鋪格點，改用 N 個隨機抽樣點估計積分，誤差以 O(N^(−1/2)) 下降，而且這個速率與維度無關。d 越大，蒙地卡羅相對於格點法的優勢越懸殊，這使它成為金融定價、統計物理、貝葉斯推斷的標配工具。

最基本的「命中或 miss」法：想算 ∫∫_D f 的面積，就把區域丟進一個外接方框，隨機撒 N 個點，數落在 D 內的比例 × 方框面積。更聰明的是樣本平均法：

∫_Ω f(x) dx ≈ V(Ω) × (1/N) Σ f(Xᵢ)，其中 Xᵢ 均勻取自 Ω，V 是 Ω 的體積。

由大數法則，這個估計量期望值恰好等於真積分值；由中央極限定理，其標準差為 V·σ_f/√N，故要精度好一個位數，樣本得多 100 倍——慢，但穩定、可平行、對維度免疫。

第三層技巧是重要性抽樣（importance sampling）：若 f 有尖峰，均勻抽樣浪費了大多數點在「不重要」的地方。改從與 f 形狀相近的機率密度 g 抽樣，再除以 g 加權，可大幅降低變異數。這一招是馬可夫鏈蒙地卡羅（MCMC）的思想起點。

演算法步驟#

確定積分區域 Ω 與其體積（或外接框）。
產生 N 個均勻隨機點 X₁, …, Xₙ ∈ Ω。
計算樣本平均 f̄ = (1/N)Σf(Xᵢ)，估計值 = V(Ω)·f̄。
用樣本標準差 s 估計誤差：E ≈ V(Ω)·s/√N。
若誤差太大，增加 N（可邊抽邊累加，無需重來）。

Python 程式碼範例#

import numpy as np

rng = np.random.default_rng(0)

# --- 例 1：用蒙地卡羅估計 π（單位圓面積/外接方框） ---
def estimate_pi(N):
    x = rng.uniform(-1, 1, N)
    y = rng.uniform(-1, 1, N)
    inside = (x**2 + y**2 <= 1).sum()
    return 4 * inside / N

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

# --- 例 2：高維積分，維度詛咒的實證 ---
def integral_mc(f, a, b, d, N):
    """d 維區間 [a,b]^d 上的樣本平均積分"""
    pts = rng.uniform(a, b, size=(N, d))
    V = (b - a) ** d
    fx = f(pts)
    mean, std = fx.mean(), fx.std()
    return V * mean, V * std / np.sqrt(N)

f = lambda x: np.exp(-np.sum(x**2, axis=1))    # e^{-|x|^2}
for d in [1, 2, 5, 10]:
    est, err = integral_mc(f, -1, 1, d, 200_000)
    print(f"d={d:>2}: 積分 ≈ {est:.6f} ± {err:.1e}")

# --- 例 3：重要性抽樣 vs 均勻抽樣 ---
def uniform_sampling(N):
    x = rng.uniform(0, 10, N)
    I = 10 * np.mean(np.exp(-x**2) * np.cos(x)**2)
    return I

def importance_sampling(N):
    # 從 g(x)=λe^{-λx}（λ=0.5）抽樣，解析反函數可得 x = -2ln(u)
    u = rng.uniform(0, 1, N)
    x = -2 * np.log(1 - u)          # 反變換法
    x = x[x <= 10]                   # 截斷近似
    g = 0.5 * np.exp(-0.5 * x)
    I = np.mean(np.exp(-x**2) * np.cos(x)**2 / g)
    return I

true_val = 0.67156458                # 由 scipy.quad 求得的參考值
for N in [10**4, 10**5, 10**6]:
    e1 = abs(uniform_sampling(N) - true_val)
    e2 = abs(importance_sampling(N) - true_val)
    print(f"N={N:>8}: 均勻誤差 {e1:.2e} | 重要性誤差 {e2:.2e}")

誤差分析#

樣本平均估計量 Îₙ = V·f̄ 是無偏的，變異數 Var(Îₙ) = V²σ²_f/N，由 CLT 得置信區間 Îₙ ± 1.96·V·σ_f/√N。因此：誤差降到 1/10 需 100 倍樣本；誤差估計本身是「免費」的——樣本標準差當場就能算出來，這是蒙地卡羅對工程師特別友善的地方。與格點法比較：d 維張量積分要 O(N^d) 點達 O(N^(−2))（辛普森），蒙地卡羅則永遠 O(N^(−1/2))，交叉點大約在 d = 3~5。降低常數因子的三招：重要性抽樣（選 g 貼近 f）、分層抽樣（把區域切成塊分別估）、擬蒙地卡羅（QMC，用低差異序列如 Sobol，誤差可達 O((log N)^d/N)）。最後切記：random 與 np.random 是偽隨機數，統計性質夠用於積分，但密碼學場景不適用——這個主題在 9.1 節會完整展開。

練習#

畫出估計 π 的誤差對 N 的 log-log 圖，驗證斜率 ≈ −1/2；再跑 100 次不同 seed 取標準差，與理論值 V·σ/√N 比對。
用蒙地卡羅估計 5 維單位球體積（外接超方框 [−1,1]⁵，真值 = 8π²/15 ≈ 5.2638），並估算要到相對誤差 1% 需要多少樣本。
對 ∫₀¹ e^(1/(x+0.1)) dx 這類尖峰函數，比較均勻抽樣與重要性抽樣（建議 g 用平移縮放後的指數分布）同誤差下的樣本數差距。