1946：ENIAC 與 Monte Carlo 誕生 — 病房裡的一局紙牌#

副標題：當中子輸運遇上接龍牌，暴力隨機成為一門科學。

案發現場#

1946 年的洛斯阿拉莫斯，空氣裡還留著戰爭的硝煙味。原子彈已經爆炸，卻沒有人真正「算」過中子在核材料裡的旅程。一顆中子誕生、散射、吸收、分裂，每一步都是機率。玻爾茲曼輸運方程寫得出來，卻解不出來，七維相空間讓解析解舉白旗。同一時間，費城的 ENIAC 剛問世，重三十噸，有一萬八千根真空管。它是為彈道表而生的巨獸，只會加減乘除，不懂什麼叫機率。案件成立：如何讓決定性的機器，算出本質隨機的中子雪崩？被害人是氫彈設計的關鍵數據，中子增殖率算不準，試爆就是賭命。 1946 年秋天，烏拉姆因腦部手術住院，只能玩接龍紙牌打發時間。他忽然冒出偵探式疑問：某種牌局能成功的機率是多少？組合數學會爆炸，但「玩一百局數成功幾次」不就行了嗎？第一條線索就此出現：用抽樣代替窮舉。

偵查過程#

線索一：從紙牌到中子#

烏拉姆把直覺寫成信寄給馮紐曼：用隨機抽樣估計解不出的積分。馮紐曼正為氫彈的中子擴散焦頭爛額，兩人一拍即合。平均值可以繞過分佈的全部細節，這是關鍵轉折。想求體積，往包圍盒丟隨機點，數落在體內的比例即可。想求中子行為，讓虛擬中子在電腦裡隨機行走即可。 | 紙牌世界 | 中子世界 | |----------|----------| | 洗牌發牌 | 抽取中子初速與方向 | | 一局牌的輸贏 | 一顆中子的吸收或洩漏史 | | 百局勝率 | 十萬顆中子的增殖率 |

線索二：ENIAC 的隨機數詭計#

ENIAC 沒有隨機數產生器，這是第二道鎖。馮紐曼發明平方取中法：平方後截取中間幾位當下一個數。方法不完美但夠快，後來被同餘產生器取代。設種子為 \(x_0\) ，乘子為 \(a\) ，模數為 \(m\) ，迭代即乘加取餘。配合均勻隨機數，反函數法可抽取任意分佈。指數分佈的自由程、散射角餘弦，都化為查表加內插。 ENIAC 為此重插接線板，操作員像電話接線生一樣工作。

線索三： \(1/\sqrt N\) 的收斂契約#

真正的突破是證明這種賭博會收斂。設待估均值為 \(I\) ，每次試驗為 \(X_i\) ，樣本平均即估計值。大數法則保證樣本平均依機率收斂到 \(I\) 。中央極限定理給出誤差帶，標準差按樣本數開方反比縮小。數量級可寫成正比於 \(1/\sqrt N\) 的形式，其中 \(N\) 為樣本數。精度提高十倍，樣本數須增百倍，代價明確。但收斂速度與維度無關，傳統網格法的高維爆炸被繞過。 1947 年首次中子試算與可解簡化情形吻合。 1949 年 Metropolis 與烏拉姆發表論文，Monte Carlo 正式命名。名字本身就是玩笑：摩納哥賭場正是隨機的聖地。

結案報告#

蒙地卡羅在 1946 至 1949 年間誕生，把積分變成頻率。中子輸運首次被電腦攻克，氫彈設計拿到關鍵參數。哲學遺產更深：決定論電腦可用來模擬機率世界。此後分岔出兩條大路，一條是 1953 年重要性抽樣。另一條是量子蒙地卡羅，乃至今日的擴散生成模型。 ENIAC 退役了，它擲出的第一把骰子仍在每台電腦裡滾動。

證據與工具#

核心估計式可用獨立公式表示如下。

\[ \hat I_N = \frac{1}{N} \sum_{i=1}^{N} f(X_i) \]

其中 \(X_i\) 為獨立同分佈樣本， \(f\) 為被積函數， \(N\) 為總數。誤差縮放關係可記為以下正比式。

\[ \sigma_{\hat I} \propto \frac{\sigma_f}{\sqrt N} \]

其中 \(\sigma_f\) 為單次標準差， \(N\) 為抽樣總數。 | 符號 | 意義 | 備註 | |------|------|------| | \(N\) | 樣本總數 | 越大越準，代價線性成長 | | \(I\) | 待求真值 | 高維積分或增殖率 | | \(X_i\) | 單次試驗 | 中子歷史或隨機點 | 最小虛擬碼：正方形內抽 \(x\) 與 \(y\) ，計 \(x^2 + y^2 \le 1\) 的點數 \(M\) 。圓周率估計為 \(4M/N\) ，誤差約 \(1/\sqrt N\) 量級。中子版把圓換成截面，把計數換成分裂中子數。關鍵文獻：Metropolis 與 Ulam 1949 年論文，ENIAC 中子計算報告。

補充：程式實作#

對應程式：1946-monte_carlo_pi.py

理論呼應#

本文斷言高維積分可用頻率代替窮舉，樣本平均即為估計值。

程式在單位正方形內投下均勻隨機點，以落入四分之一圓的比例推估圓周率。

其估計式正是正文核心公式 \(4M / N\) ，其中 \(M\) 為命中數而 \(N\) 為總數。

誤差依 \(1 / \sqrt N\) 收縮且與維度無關，這正是蒙地卡羅繞過網格爆炸的關鍵。

嫌犯藏在高維相空間裡，唯有隨機撒網能逼它現出原形，這一局賭得有理。

如何執行#

在本文所屬目錄下執行以下指令即可重現全部數字。

python3 _code/1946-monte_carlo_pi.py

程式只依賴 numpy ，已固定亂數種子，任何機器重跑都應一致。

預設樣本數為 N = 200000 ，統計漲落屬於正常現象而非程式錯誤。

實測結果抄錄#

以下數字為本次實際執行所得，未經修飾，直接抄錄關鍵行。

N=200000 inside=157168 pi_est=3.143360 true=3.141593 err=1.767346e-03

理論標度 1/sqrt(N)=2.236068e-03

理論標準差約 3.672034e-03

err/(1/sqrt(N)) = 0.7904

驗證行顯示誤差小於 0.01 且小於五倍理論標度，程式最後印出 PASS 。

由此可見誤差確實落在中央極限定理預言的量級之內。

動手實驗#

讀者可把 N 改為 20000 或 2000000 ，觀察誤差是否大致按十倍樣本縮小約三倍。

完整程式如下：

# -*- coding: utf-8 -*-
"""1946 Monte Carlo 估 pi（對應 wiki：計算模擬學 / Monte Carlo 方法）
在 [0,1]² 均勻投 N=2e5 點，落入 x²+y²≤1（1/4 圓）比例 p→pi/4，
pi_est=4·count/N。統計誤差標度 ~1/sqrt(N)，
理論標準差 4·sqrt(p(1-p)/N)≈1.64/sqrt(N)。
驗證 |est-pi|<0.01。
只用 numpy，固定種子。
"""
import numpy as np

np.random.seed(0)

N = 200000
xy = np.random.rand(N, 2)
inside = np.sum(xy[:, 0] ** 2 + xy[:, 1] ** 2 <= 1.0)
pi_est = 4.0 * inside / N
err = abs(pi_est - np.pi)
scale = 1.0 / np.sqrt(N)
p = np.pi / 4.0
theory_std = 4.0 * np.sqrt(p * (1.0 - p) / N)

print(f"N={N} inside={inside} pi_est={pi_est:.6f} true={np.pi:.6f} err={err:.6e}")
print(f"理論標度 1/sqrt(N)={scale:.6e}, 理論標準差≈{theory_std:.6e} (≈1.64/sqrt(N))")
print(f"err/(1/sqrt(N)) = {err / scale:.4f} (應為 O(1) 量級)")

# 驗證數字：理論值 vs 實測
ok = err < 0.01
ok_scale = err < 5.0 * scale
print(f"VERIFY err={err:.6e} (<0.01? {ok})")
print(f"VERIFY err={err:.6e} (<5/sqrt(N)={5*scale:.6e}? {ok_scale})")
assert ok, f"Monte Carlo 誤差太大: {err}"
assert ok_scale, f"誤差偏離 1/sqrt(N) 標度: {err / scale:.2f}σ單位"
print("PASS")