跳轉到

5.3 大數法則與蒙地卡羅方法#

本節目標#

  • 分辨大數法則（LLN）與中央極限定理（CLT）：平均值收斂 vs 收斂到什麼形狀
  • 用 LLN 估 \(\pi\) ：Buffon 投針（1777）與四分之一圓法
  • 學會蒙地卡羅積分的基本套路
  • 認識「蒙地卡羅」名字的由來（Ulam, 1946）
  • 用模擬驗證蒙地卡羅誤差 \(\sim 1/\sqrt{N}\)

問題：可以用「丟亂數」算數學嗎？#

5.2 節說過：LLN 保證 \(\bar{X} \xrightarrow{P} \mu\) ——樣本平均會收斂到母體平均。這句話反過來讀就是一套計算方法：

想算一個數 \(\theta\) ？若能把 \(\theta\) 寫成某個隨機變數的期望值 \(\theta = E[Y]\) ，就大量抽 \(Y_1, \ldots, Y_N\) ，用 \(\bar{Y}\) 當答案。 \(N\) 越大越準——這就是蒙地卡羅方法。

LLN 回答「會不會準」（會，機率收斂），CLT 回答「誤差多大」—— \(\bar{Y}\) 的波動約 \(\sigma_Y/\sqrt{N}\) ，所以每多花 100 倍算力，只換到 10 倍精度。這條 \(1/\sqrt{N}\) 定律是蒙地卡羅方法的宿命，本節最後會實測它。

歷史：從蒲豐的針到蒙地卡羅#

  • 1777 — Buffon 投針問題。法國博物學家蒲豐（Georges-Louis Leclerc, Comte de Buffon）提出：在間距 \(D\) 的平行線上隨機投一根長 \(L \le D\) 的針，針與線相交的機率是 $$ p = \frac{2L}{\pi D}

$$ 這是史上第一個「用隨機實驗逼近數學常數」的範式：反解出 \(\pi = \frac{2LN}{D \cdot \text{hits}}\) 。比電腦早了近兩百年。 - 1946 — 「蒙地卡羅」誕生。波蘭裔數學家 Ulam（Stanisław Ulam）在 Los Alamos 養病時玩接龍牌局，想算中獎機率，發現窮舉不可行、改用「重複抽樣統計次數」就有希望。他與 von Neumann 把這方法用在中子輸運計算，同事 Nicholas Metropolis 以摩納哥的賭城蒙地卡羅（Ulam 叔父常去借钱賭博的地方）命名——「用偶然性計算確定性的量」從此有了名字。

理論：兩個估 \(\pi\) 的配方#

四分之一圓法：在 \([0,1] \times [0,1]\) 均勻撒點 \((X, Y)\) ，落在四分之一圓內的機率是圓面積：

$$ P(X^2 + Y^2 \le 1) = \frac{\pi}{4} \quad\Rightarrow\quad \hat{\pi} = 4 \cdot \frac{\text{圓內點數}}{\text{總點數}}

$$ 蒙地卡羅積分：想算 \(I = \int_a^b f(x)\,dx\) ，把 \(f\) 當成某均勻隨機變數的期望值：

$$ I = (b-a) \cdot E[f(U)], \quad U \sim U(a, b) \quad\Rightarrow\quad \hat{I} = (b-a)\cdot\frac{1}{N}\sum_{k=1}^N f(U_k)

$$

動手（一）：LLN 現場——骰子平均的收斂#

import math
import rstatistics as r

# ---- (1) LLN：擲公平骰子的樣本平均收斂 ----
g = r.set_seed(3)
rolls = [float(v) for v in g.sample(6, 10000, replace=True)]  # 1..6
print("擲公平骰子，樣本平均隨 n 變化（母體平均 = 3.5）：")
running = 0.0
for n, v in enumerate(rolls, 1):
    running += v
    if n in (10, 100, 1000, 10000):
        print(f"  n = {n:5d}   xbar = {running/n:.4f}")

# ---- (2) Buffon 投針 (1777) 估 pi ----
print()
L = 1.0   # 針長
D = 2.0   # 平行線間距（>= L）
N = 200000
g = r.set_seed(11)
hits = 0
for _ in range(N):
    x = g.runif(1, 0.0, D / 2)[0]        # 針中點到最近線的距離
    theta = g.runif(1, 0.0, math.pi / 2)[0]
    if x <= (L / 2) * math.sin(theta):
        hits += 1
pi_hat = (2 * L * N) / (D * hits)
print(f"Buffon 投針 N={N}: 相交 {hits} 次, 估計 pi = {pi_hat:.4f}")

輸出（實際執行）：

擲公平骰子，樣本平均隨 n 變化（母體平均 = 3.5）：
  n =    10   xbar = 3.3000
  n =   100   xbar = 3.3900
  n =  1000   xbar = 3.4750
  n = 10000   xbar = 3.4827

Buffon 投針 N=200000: 相交 63666 次, 估計 pi = 3.1414

樣本平均一步步滑向 3.5（LLN），20 萬根針把 \(\pi\) 鎖到小數第二位（Buffon）。注意骰子在 \(n=100\) 時還差 0.11、 \(n=10000\) 時差 0.017——縮小的比例正是 \(\sqrt{10}\) 倍左右，CLT 的 \(1/\sqrt{N}\) 在背後運作。

動手（二）：四分之一圓法、蒙地卡羅積分與 \(1/\sqrt{N}\) 定律#

import math
import rstatistics as r

def quarter_circle(N, seed):
    g = r.set_seed(seed)
    hits = 0
    for _ in range(N):
        x = g.runif(1, 0.0, 1.0)[0]
        y = g.runif(1, 0.0, 1.0)[0]
        if x * x + y * y <= 1.0:
            hits += 1
    return 4.0 * hits / N

print("=== 四分之一圓法估 pi ===")
for N in (1000, 10000, 100000, 1000000):
    print(f"N = {N:8d}   pi_hat = {quarter_circle(N, 5):.4f}   誤差 = {abs(quarter_circle(N, 5) - math.pi):.4f}")

# ---- 蒙地卡羅估積分 ∫0^1 e^x dx = e - 1 ----
print()
print("=== 蒙地卡羅估積分 I = ∫_0^1 e^x dx = e - 1 ≈ 1.718282 ===")
def mc_integral(N, seed):
    g = r.set_seed(seed)
    tot = 0.0
    for _ in range(N):
        tot += math.exp(g.runif(1, 0.0, 1.0)[0])
    return tot / N

for N in (1000, 10000, 100000, 1000000):
    print(f"N = {N:8d}   I_hat = {mc_integral(N, 8):.6f}   誤差 = {abs(mc_integral(N, 8) - (math.e - 1)):.6f}")

# ---- 誤差 ~ 1/sqrt(N) 的驗證：重複 500 次估同一積分，看各 N 的標準差 ----
print()
print("=== 誤差標準差 vs 1/sqrt(N)（每個 N 重複 500 次，seed 1..500）===")
for N in (1000, 10000, 100000):
    errs = []
    for seed in range(1, 501):
        errs.append(mc_integral(N, seed) - (math.e - 1))
    sd = r.sd(errs)
    print(f"N = {N:7d}: 誤差的 sd = {sd:.6f}   0.5/sqrt(N) = {0.5/N**0.5:.6f}   比值 = {sd*N**0.5:.4f}")

輸出（實際執行）：

=== 四分之一圓法估 pi ===
N =     1000   pi_hat = 3.1240   誤差 = 0.0176
N =    10000   pi_hat = 3.1492   誤差 = 0.0076
N =   100000   pi_hat = 3.1506   誤差 = 0.0090
N =  1000000   pi_hat = 3.1431   誤差 = 0.0015

=== 蒙地卡羅估積分 I = ∫_0^1 e^x dx = e - 1 ≈ 1.718282 ===
N =     1000   I_hat = 1.714657   誤差 = 0.003625
N =    10000   I_hat = 1.719690   誤差 = 0.001408
N =   100000   I_hat = 1.718481   誤差 = 0.000199
N =  1000000   I_hat = 1.718043   誤差 = 0.000239

=== 誤差標準差 vs 1/sqrt(N)（每個 N 重複 500 次，seed 1..500）===
N =    1000: 誤差的 sd = 0.015670   0.5/sqrt(N) = 0.015811   比值 = 0.4955
N =   10000: 誤差的 sd = 0.004944   0.5/sqrt(N) = 0.005000   比值 = 0.4944
N =  100000: 誤差的 sd = 0.001529   0.5/sqrt(N) = 0.001581   比值 = 0.4836

讀結果：

  • 單次估計的誤差有起有落（四分之一圓法 \(N=10^5\) 那次的誤差反而比 \(10^4\) 大）——蒙地卡羅的誤差是隨機的，一次看不出定律。
  • 「誤差的 sd」那張表才是定律本尊： \(N \times 100\) 時 sd 恰好除以 10，且 \(\text{sd}\cdot\sqrt{N} \approx 0.49\) 幾乎不變。理論上這個常數是 \(\sqrt{\mathrm{Var}(f(U))} = \sqrt{\mathrm{Var}(e^U)} \approx 0.4958\) ——模擬與 CLT 完全吻合。

為什麼還要用蒙地卡羅算積分？ \(\int_0^1 e^x dx\) 用數值積分（梯形法）遠比蒙地卡羅準。但一旦維度變高——積分變成 \(\int_{[0,1]^d}\) ——數值積分的格點數是指數成長（維度災詛咒），蒙地卡羅的 \(1/\sqrt{N}\) 卻與維度無關。Ulam 當年算的中子輸運正是這種高維問題；這也是它至今是財務、物理、貝葉斯統計主力工具的原因。

與 R 對照#

set.seed(5)
N <- 100000
pts <- matrix(runif(2 * N), ncol = 2)
4 * mean(pts[, 1]^2 + pts[, 2]^2 <= 1)

set.seed(8)
mean(exp(runif(N)))          # MC 估 ∫0^1 e^x dx
exp(1) - 1

輸出（實際執行）：

[1] 3.14112
[1] 1.718481
[1] 1.718282

小結#

  1. LLN： \(\bar{Y} \xrightarrow{P} E[Y]\) ——蒙地卡羅方法「為什麼會準」的全部根據。
  2. CLT：誤差尺度 \(\sigma_Y/\sqrt{N}\) ——「有多準、值多少算力」的根據；精度一個數位要 100 倍樣本。
  3. Buffon（1777）示範了「隨機實驗算常數」的範式；Ulam（1946）把它變成系統性方法並得名蒙地卡羅。
  4. 蒙地卡羅積分在高維度打敗數值積分，因為 \(1/\sqrt{N}\) 不怕維度災詛咒。
  5. 本章所有「模擬驗證」其實都是蒙地卡羅——你已經用了一整章了。

練習#

  1. 用四分之一圓法估 \(\pi\) 到誤差 0.001 以內，約需要多少 \(N\) ？先用 CLT 預測（ \(Y = 4\cdot\mathbf{1}\{X^2+Y^2\le1\}\) 的 sd 是 \(\pi\sqrt{4-\pi}/2 \approx 1.642\) ），再實測。
  2. 把 Buffon 的 \(N\) 改成 20000 與 2000000，估計值與 \(1/\sqrt{N}\) 的預測一致嗎？ （提示：注意 Buffon 估的是「相交次數」的比例，誤差常數與 \(p\) 有關。）
  3. 用蒙地卡羅估 \(\int_0^1 \cos(x^2)\,dx\) （此積分沒有初等封閉形式）， \(N = 10^5\) ，並由「500 次重複的 sd」給出誤差範圍。 （提示：改被積函數即可；誤差範圍寫成 \(\hat{I} \pm 2\,\text{sd}\) 。）
  4. 比較「一次 \(N=10^6\) 」與「100 次各 \(N=10^4\) 再平均」的估計精度，哪個好？誤差理論怎麼解釋？ （提示：兩者樣本總量一樣；由 CLT 兩者 sd 都 \(\approx 0.5/\sqrt{10^6}\) ——但後者還能給出誤差估計，這是分批跑的隱藏好處。）