Monte Carlo Simulation_求圓周率π

- 8月 08, 2026

蒙特(地)卡羅模擬

是科學計算中最重要且最普遍的數值方法之一。利用電腦模擬一個理想環境，然後利用亂數產生大量隨機狀況，藉以探討許多難以進行實際試驗的理論。
最早應用於曼哈頓計畫（Manhattan Project），之後也被用於氫彈的研發。
後續時常用在經濟學、物理學(量子、熱力)、機器學習等領域
蒙地卡羅是摩納哥親王國最著名的一區，以豪華的賭場聞名於世。
20世紀40年代，在科學家馮·紐曼、斯塔尼斯拉夫·烏拉姆、尼古拉斯·梅特羅波利斯三位科學家在洛斯阿拉莫斯國家實驗室為核武器計劃工作時，發明了蒙地卡羅方法。
會取這個名稱，是因為烏拉姆的叔叔在摩納哥的蒙地卡羅賭場輸了很多錢。
1940年代美國研發核子武器時，採用電腦模擬進行爆炸威力研究，這個方法運用到機率及亂數來模擬中子碰撞所產生的能量，頗有賭博的意味，所以就用「蒙地卡羅」來命名。

優點

廣泛應用於量化金融（Quantitative Finance），近年隨著 GPU 的發展而更加普及。
當評價問題中包含大量隨機因素，多到無法直接使用一般數值方法進行評價時，蒙地卡羅模擬特別有用。
當標的變數的機率分布較為複雜，使得直接求解相當困難時，也適合使用蒙地卡羅模擬。
適合處理**報酬取決於標的資產價格路徑（Path-dependent）**的金融商品。

缺點

使用蒙地卡羅方法評價美式選擇權（American Options）並不容易，雖然仍有方法可以做到。
計算效率相對較低，通常需要進行大量模擬，才能使選擇權價格逐漸收斂到較穩定的結果。

求圓周率π

圓的周長和其直徑的比率，約等於3.14159265358979323846264....
它在18世紀中期之後一般用希臘字母 π 來表示
中國南宋朝數學家祖沖之，曾經用幾何方法將圓周率計算到小數點後7位數字。因此，數學界又將圓周率敬稱祖率。

測試程式碼(原子分布狀態模擬10000個點)

import numpy as np
import pylab as plt

batch = 10000
#隨機產生 10000 個二維座標點 (x, y)
#其中 x、y 都介於 0～1 之間，最後畫成散佈圖（Scatter Plot）。
#uniform(平均分布)->位於0~1之間機率都是一樣的
xs = np.random.uniform(0,1,batch)
ys = np.random.uniform(0,1,batch)
plt.scatter(xs, ys,1)
plt.show()

#normal(常態分布)->越接近中間0的位置機率越大，越遠則機率越小。
xs = np.random.normal(0,1,batch)
ys = np.random.normal(0,1,batch)
plt.scatter(xs, ys,1)
plt.show()

uniform(平均分布)->位於0~1之間機率都是一樣的

normal(常態分布)->越接近中間0的位置機率越大，越遠則機率越小。

圓的面積公式 = πr²

假設我們有一個單邊邊長為2r的矩形，矩形內有一個圓形，矩形四個邊與內圓相切。

以r作為圓的半徑，矩形面及為(2r)²=4r²，圓面積是 πr²

圓形與矩形的面積比就是π/4

如果我們能夠計算出這兩個面積的比例，那麼只要將這個比例乘以 4，就可以得到圓周率 π。

我們也可以在正方形範圍內隨機產生繪製一些點，再去計算落在圓內的點數和總點數的比。

蒙地卡羅求解圓周率測試程式

用1000萬測試200次

import numpy as np
import pylab as plt
import time

incircle=0
epochs = 200
batch = 100_000_000#一億個點
for e in range(epochs):
  t1 = time.time()
  #points = [x值位於0~1的隨機數,y值位於0~1的隨機數]
  points = [np.random.uniform(0,1,batch),np.random.uniform(0,1,batch)]
  dis = np.sqrt(np.square(points[0])+np.square(points[1]))
  #print("np.where(dis<=1):",np.where(dis<=1))
  #print("np.where(dis<=1)[0]:",np.where(dis<=1)[0])
  #print("np.where(dis<=1)[0].shape:",np.where(dis<=1)[0].shape)
  count=np.where(dis<=1)[0].shape[0]
  #print("in circle point count:",count)
  #print("==============================")
  incircle+=count
  pi=incircle/((e+1)*batch)*4
  t2 = time.time()
  print(f"Epoch:{e+1:03d}=>{t2-t1:.5f}秒,pi={pi}")