用蒙地卡羅方法來計算圓周率

一、介紹蒙地卡羅方法

當一個問題太複雜、沒辦法用公式直接算出答案時，該怎麼辦？蒙地卡羅方法（Monte Carlo Method）的答案很簡單：算不出來，就「試出來」。用隨機抽樣反覆模擬成千上萬次，再從結果的統計分布去逼近答案；模擬次數越多，結果通常就越準。它給的不是一個絕對精確的數字，而是「最可能的答案」以及「可能的範圍」——這是一種用隨機性去馴服複雜性的方法。

由來：1946 年，數學家斯坦尼斯瓦夫．烏拉姆（Stanisław Ulam）在病中玩單人紙牌遊戲，想算出贏牌的機率，他發現與其埋頭推 導，不如實際玩很多次再統計。之後他與約翰．馮紐曼（John von Neumann）在洛斯阿拉莫斯國家實驗室，把這個想法用在中子擴散等核物理計算上——這些問題正是源自曼哈頓計畫的研究。同事尼古拉斯．梅特羅波利斯（Nicholas Metropolis）為它取了「蒙地卡羅」這個代號，名稱 來自摩納哥以賭場聞名的蒙地卡羅，因為方法的核心就是「機率」與「隨機」。而 ENIAC 等早期電子計算機的出現，讓大量模擬真正變得可行 ；1949 年 Metropolis 與 Ulam 發表了論文〈The Monte Carlo Method〉。

基本步驟：

1. 定義問題中的不確定因素與它們的機率分布。
2. 依照這些分布隨機抽樣，產生一組情境。
3. 對每個情境計算結果。
4. 重複數千到數百萬次。
5. 統計所有結果，得出平均值、分布與機率。

它的優點是能處理高維度、非線性、沒有公式解的問題，而且概念直覺；限制則是需要大量運算，而且結果好壞取決於假設的機率分布是否 正確 —— 假設錯了，模擬再多次也是錯的（垃圾進，垃圾出）。

二、應用領域

• 物理：中子擴散、輻射傳遞（最早的用途）、粒子物理實驗的偵測器模擬。
• 醫學：放射治療的劑量計算。
• 工程可靠度：將零件公差、材料強度、負載等不確定性隨機抽樣，估算系統失效機率與誤差範圍。
• 電腦圖學：路徑追蹤（Path Tracing）以隨機取樣光線路徑來計算光照積分，產生逼真的影像，電影動畫與遊戲渲染都大量使用。
• 人工智慧：蒙地卡羅樹搜尋（MCTS）以大量隨機對局評估棋步優劣，是 AlphaGo 的核心技術之一。
• 統計與貝氏推論：馬可夫鏈蒙地卡羅（MCMC，例如 Metropolis–Hastings 演算法、Gibbs 抽樣）用來從複雜的後驗分布中抽樣。
• 氣象與氣候：系集預報，以多組略有差異的初始條件反覆模擬。
• 複雜系統：模擬大量個體互動的行為。

三、金融領域應用簡介

• 選擇權定價：假設股價服從幾何布朗運動（GBM），模擬成千上萬條可能的價格路徑，計算每條路徑到期時的報酬（payoff），以無風險利率折現後取平均，即得到選擇權價格的估計。對於路徑相依的選擇權（例如以期間平均價格結算的亞式選擇權 Asian option>），往往沒有簡單的封閉解，蒙地卡羅方法特別好用。
• 風險值（VaR）與風險管理：模擬投資組合在未來一段時間的損益分布，取其分位數（例如 99%）來估算「在一定信心水準下最多可能損失多少」，也可延伸計算預期損失（CVaR）與壓力測試。
• 投資組合與退休規劃：以隨機報酬模擬未來數十年的資產變化，估算「資金撐到退休後 30 年」的成功機率，協助決定提領率與資產配置。

⚠️ 提醒：金融模擬常假設報酬呈 常態分布，但真實市場存在肥尾與極端事件（黑天鵝），模型假設是否合理，與模擬次數同樣重要。

四、計算圓周率的原理

想像一個邊長為 1 的正方形，範圍是 [0,1]×[0,1]。以原點為圓心、半徑 1 畫圓，落在這個正方形裡的正好是四分之一圓。

四分之一圓面積 ÷ 正方形面積 = (π·1²/4) ÷ 1² = π/4

為什麼用四分之一圓？完整的做法是在 2×2 的正方形 [−1,1]×[−1,1] 裡放一個半徑 1 的整圓，面積比為 π/4。但圓與正方形都對 稱於兩軸，四個象限的情況完全相同，所以只取第一象限（x, y ∈ [0,1]）比例不變，而且亂數直接取 0～1 之間即可，更簡單。

一個隨機點 (x, y) 落在四分之一圓內的條件是：

x² + y² ≤ 1

隨機丟 N 個點，數出落在圓內的點數，則「圓內點數 / N」會逼近 π/4，因此：

π ≈ 4 × 圓內點數 ÷ N

五、標準差與誤差

每一個點只有兩種結果：落在圓內（機率 p = π/4 ≈ 0.7854）或圓外，這是一次伯努利試驗。N 個點中圓內比例的變異數為 p(1−p)/N，而估計值是比例乘以 4，所以估計值的標準差為：

σ = 4 × √( p(1−p) / N ) ≈ 1.642 / √N

誤差與 1/√N 成正比：點數增加 100 倍，精度才提高 10 倍（多一位有效數字）。這正是蒙地卡羅方法的特性——簡單通用，但收斂 較慢。

依常態分布近似，真實值落在估計值附近的機率：

範圍	機率
±1σ	68.3%
±2σ	95.4%
±3σ	99.7%

程式中顯示的 ssss ≤ xxxx ≤ llll，機率 99.7%，就 是「估計值 ± 3σ」的區間：中間是估計值，左右是下限與上限。

• 若希望整個 ±3σ 區間寬度（6σ）小於 0.01，需要 N ≈ (6×1.642/0.01)² ≈ 97 萬個點。
• 若希望上下限都以 3.14 開頭（下限 ≥ 3.14，即 3σ ≤ π−3.14 ≈ 0.00159），需要約 960 萬個點。

💡 小插曲：品質管理中的「六標準差（Six Sigma）」指規格界限距平均值 ±6σ；考慮製程長期約 1.5σ 的偏移後，不良率約為每百萬個 3.4 個。

六、程式與線上執行

我用 Kotlin Multiplatform 寫了一個小程式，同時支援桌面版與網頁版：

• 原始碼：https://github.com/neojou/MonteCarlo_Pi
• 線上執行：https://neojou.github.io/MonteCarlo_Pi/

功能包括：以滑桿調整點數與速度、Start／Pause／Reset 控制、可指定亂數種子（seed）以重現結果、即時顯示估計值、σ、誤差等 KPI>，以及收斂圖表。

七、執行結果

圖一：N = 10,000,000、seed = 0。估計值 3.141047，區間 3.139489 ≤ 3.141047 ≤ 3.142605，機率 99.7%；圓內點數 7,852,617，σ = 0.000519，誤差 0.000546，真實 π 落在區間內 ✓

圖二：收斂圖（捲動檢視）。橫軸為點數 N（對數座標 10¹～10⁷），>縱軸為估計值；±3σ 帶隨 N 增加而變窄，估計值逐漸收斂到 π。

解讀：4 × 7,852,617 ÷ 10,000,000 = 3.1410468；σ = 1.642/√10⁷ ≈ 0.000519，±3σ 即 ±0.001558。實際誤差 |3.1410468 − π| ≈ 0.000546，約為 1.05σ，完全在預期的統計波動範圍內。收斂圖也清楚呈現 1/√N 的特性：每往右一個數量級，帶寬約 縮小為 1/√10 ≈ 0.32 倍。

八、結語

從烏拉姆的一副紙牌，到核物理、電腦圖學、AlphaGo 與金融風險管理，蒙地卡羅方法用「大量隨機」換來「可靠的近似」。用它計算圓周 率雖然不是最有效率的方法，卻是理解這套思想最直覺的入口：丟一千萬個點，我們不只得到 3.141047，還知道這個答案有多可信。歡迎到線上版親自試試，看看 π 如何在隨機中慢慢浮現。