1953：Metropolis Monte Carlo — MANIAC 上的狀態方程#

副標題：按玻爾茲曼的指引，去重要的地方抽樣。

案發現場#

1953 年，洛斯阿拉莫斯 MANIAC 電腦前站著五個人。 Metropolis、Rosenbluth 夫婦、Teller 夫婦，盯著二維流體狀態方程。問題具體：給定溫度密度，壓強與內能是多少？均勻抽樣在此撞牆：高能位形量多但權重小。直接抽樣抽不到有貢獻的低能區，估計抖動巨大。癥結是維度災難加玻爾茲曼因子的尖銳峰。配分函數是 \(3N\) 維積分， \(N\) 為粒子數，網格毫無機會。均勻蒙地卡羅方差太大， \(1/\sqrt N\) 的常數因子太大了。偵探需要新抽樣：自動走向重要的地方。

偵查過程#

線索一：重要性採樣的直覺#

平衡分佈正比於玻爾茲曼因子，能低權重大。與其均勻撒網，不如讓隨機行走服從目標分佈。期望值變成沿馬可夫鏈的算術平均。關鍵是構造鏈，使平穩分佈恰為玻爾茲曼分佈。細緻平衡是充分條件，也是破案鑰匙。設躍遷機率為 \(W\) ，目標機率為 \(p\) ，流量須兩兩相抵。满足細緻平衡，鏈跑久必收斂到 \(p\) 。

線索二：接受率的妙手#

方案分兩步：提議加接受，隨機位移一粒子。算能量差 \(\Delta E\) ，降低直接接受，升高按機率接受。接受率行內寫成 \(\min(1, e^{-\beta\Delta E})\) 的形式。其中 \(\beta\) 為逆溫度， \(\Delta E\) 為試探能量差。設計恰满足細緻平衡，兩行代數即可驗證。設提議對稱，接受率之比等於目標機率之比。溫度藏在 \(\beta\) 裡：高溫敢爬坡，低溫只走下坡。

線索三：MANIAC 的實證#

模擬二維硬盤與方阱流體，粒子數約數百。步長須調校：太大一直被拒，太小走得太慢。接受率三到五成是經驗甜蜜點。數萬步後能量壓強趨穩，壓強出現疑似相變拐折。低密度與維里展開吻合，高密度給出全新預言。論文寥寥數頁，卻是史上被引用最多的卷宗之一。平衡證明：正逆向接受率代入，細緻平衡逐對成立。

結案報告#

1953 年算法破解平衡態抽樣之謎，方差常數壓小數個量級。狀態方程第一次由電腦直接「量」出，而非解析近似。 Hastings 1970 年推廣為非對稱提議的 MH 算法。 Gibbs 抽樣、叢集算法、混合蒙地卡羅接連誕生。貝葉斯統計、Ising 模型、蛋白摺疊處處都是它的指紋。結語：聰明的隨機勝過勤奮的窮舉。

證據與工具#

接受率獨立公式如下，為核心證物。

\[ A(o \to n) = \min(1, e^{-\beta\Delta E}) \]

其中 \(\Delta E\) 為能量差， \(\beta\) 為逆溫度， \(A\) 為接受機率。目標分佈獨立公式如下，為平衡態通緝令。

\[ p(x) = \frac{e^{-\beta E(x)}}{Z} \]

其中 \(E\) 為位形能量， \(Z\) 為配分函數， \(x\) 為全部座標。 | 符號 | 意義 | 調校建議 | |------|------|----------| | \(\beta\) | 逆溫度 | 等於 \(1/k_BT\) ，低溫更難 | | \(\Delta E\) | 能量差 | 負值直接接受 | | \(A\) | 接受率 | 維持約 0.3 至 0.5 | 虛擬碼：初始化 \(x\) 算 \(E\) ，位移一粒子得 \(\Delta E\) 。若 \(\Delta E \le 0\) 接受，否則以 \(e^{-\beta\Delta E}\) 機率接受。捨棄前段熱化，後段取平均即得狀態方程。因果鏈：1946 年均勻抽樣開案，本案升級為重要性抽樣。 Hastings 1970 年推廣非對稱提議，適用範圍大增。 Gibbs 抽樣逐座標更新，高維貝葉斯推斷全靠它。叢集算法整塊翻轉自旋，臨界慢化得以緩解。 1987 年直方圖重加權，一個溫度跑出全相圖。本案檔案編號 MC-1953 ，狀態已結案，抽樣仍在繼續。調校格言：接受率太高走太慢，太低原地踏步。五人署名論文僅數頁，引用數卻居蒙地卡羅之冠。 MANIAC 之名正是數學分析儀，名副其實的第一案發現場。

補充：程式實作#

對應程式：1953-metropolis_ising.py

理論呼應#

本文斷言聰明抽樣勝過均勻撒網，關鍵是讓鏈的平穩分佈等於目標分佈。

程式以二維伊辛模型示範重要性抽樣，提議翻轉加接受率決定去留。

核心接受率為 \(\min(1, e^{-\beta \Delta E})\) ，能量下降直行而上升按機率放行。

逆溫度 \(\beta\) 藏著高低溫的性格，低溫有序而高溫無序正好呼應相變直覺。

證物只有一串磁化強度起伏，而細緻平衡正是讓它開口招供的審訊技巧。

如何執行#

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

python3 _code/1953-metropolis_ising.py

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

預設為 20x20 格點與 2000 次掃掠，低溫與高溫各跑一條鏈。

實測結果抄錄#

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

T=1.0 |M|=1.0000

T=4.0 |M|=0.2250

VERIFY: M_low=1.000000, M_high=0.225000

驗證行顯示低溫磁化大於零點九而高溫小於零點四。

程式最後印出 ALL CHECKS PASS ，表示有序與無序兩相皆被正確捕捉。

由此可見接受率設計確實把隨機行走牽引到玻爾茲曼權重最大的區域。

動手實驗#

讀者可把溫度 1.0 改為 2.0 或 2.3 附近，觀察磁化如何從有序滑向臨界起伏。

完整程式如下：

# -*- coding: utf-8 -*-
"""1953 Metropolis + Ising 2D 模型 (對應 wiki: 計算模擬學 / 蒙地卡羅統計物理)
只用 numpy,固定種子。20x20,T=1.0(有序)與 T=4.0(無序),Metropolis 各 2000 sweeps。
驗證:低溫 |M|>0.9、高溫 |M|<0.4。
"""
import numpy as np

np.random.seed(0)
L = 20
SWEEPS = 2000

def run_ising(T):
    spins = np.random.choice(np.array([-1, 1]), size=(L, L))
    # dE 只可能為 -8,-4,0,4,8,預先算接受率
    p4 = np.exp(-4.0 / T)
    p8 = np.exp(-8.0 / T)
    for _ in range(SWEEPS):
        ii = np.random.randint(0, L, size=L * L)
        jj = np.random.randint(0, L, size=L * L)
        rr = np.random.rand(L * L)
        for k in range(L * L):
            i = int(ii[k]); j = int(jj[k])
            s = spins[i, j]
            nb = spins[(i + 1) % L, j] + spins[(i - 1) % L, j] + spins[i, (j + 1) % L] + spins[i, (j - 1) % L]
            dE = 2 * s * nb
            if dE <= 0:
                spins[i, j] = -s
            elif dE == 4:
                if rr[k] < p4:
                    spins[i, j] = -s
            else:  # dE == 8
                if rr[k] < p8:
                    spins[i, j] = -s
    mag = abs(float(np.sum(spins))) / (L * L)
    return mag

m_low = run_ising(1.0)
m_high = run_ising(4.0)
ok_low = m_low > 0.9
ok_high = m_high < 0.4
print(f"T=1.0 |M|={m_low:.4f} (驗證 >0.9: {'PASS' if ok_low else 'FAIL'})")
print(f"T=4.0 |M|={m_high:.4f} (驗證 <0.4: {'PASS' if ok_high else 'FAIL'})")
print(f"VERIFY: M_low={m_low:.6f}, M_high={m_high:.6f}")
assert ok_low, "低溫未有序"
assert ok_high, "高溫未無序"
print("ALL CHECKS PASS")