用 NumPy 實作蒙特卡羅方法:快速估算 π 值的實戰指南

作者: Calpa Liu
字數:1808
出版:April 21, 2025
分類: 後端開發 Python Numpy

想用最直觀的方法學會蒙特卡羅模擬?本篇教你如何結合 Python 的科學計算利器 NumPy,實作經典案例「估算 π 值」。從原理解析到圖形化展示,讓你快速掌握高效模擬技巧,適用於資料科學與數值分析入門。

蒙特卡羅模擬的基本概念

蒙特卡羅模擬是一種強大的計算技術,透過隨機抽樣來建模並分析具有不確定性的複雜系統。該方法得名於摩納哥著名的蒙特卡羅賭場,象徵著機會與隨機性,廣泛應用於金融、工程、物理、經濟學與風險分析等領域。

在實作上,分析者會先建立一個數學模型,模型中包含可能變動的輸入變數。透過大量的隨機抽樣與模擬,每一次迭代代表系統的一種可能狀況。透過累積模擬結果,可以估算系統的整體行為與結果分佈,特別適合處理具高不確定性或無法解析求解的問題。

蒙特卡羅模擬的工作流程

蒙特卡羅模擬通常遵循以下步驟:

  1. 定義問題與設定目標:確定要分析的專案方面,如成本、進度或績效,並為分析建立特定目標。
  2. 確定輸入變量及其分佈:這些變量可能包括任務持續時間、成本、資源可用性以及與專案風險相關的不確定性。還必須定義每個風險的概率和影響評分。
  3. 創建蒙特卡羅模擬模型:一旦確定了輸入變量及其分佈,可以使用專業軟體或工具創建蒙特卡羅模擬模型,這些工具可能包括電子表格、統計軟體包或專用的蒙特卡羅模擬軟體。
  4. 運行模擬並分析結果:模擬運行後,結果可以用來了解專案的可能結果範圍。
  5. 驗證和細化模型:最後,必須通過將模擬結果與實際專案數據或類似專案的歷史數據進行比較,以確保模型準確地表示專案的複雜性和不確定性。

如何用蒙特卡羅方法估算 π 值

使用蒙特卡羅方法估算 π 值是蒙特卡羅方法的一個經典應用。其基本原理是:

  1. 在一個邊長為 2 的正方形(範圍為[-1,1]×[-1,1])中畫一個半徑為 1 的圓
  2. 在正方形內隨機投擲大量點
  3. 計算落在圓內的點數佔總點數的比例
  4. 這個比例乘以 4 即為π的估計值

這是因為圓的面積為πr²,而正方形的面積為 4r²(其中 r=1),所以圓與正方形面積之比為π/4。通過蒙特卡羅方法,我們用落在圓內的點數與總點數之比來近似這個面積比。

用 NumPy 高效實作蒙特卡羅模擬

NumPy 是 Python 科學計算的基礎庫,為處理大型多維數組提供了高效的工具。在蒙特卡羅模擬中,NumPy 的優勢尤為明顯。

使用 NumPy 實現蒙特卡羅法估算π值的完整示例

以下是一個使用 NumPy 實現蒙特卡羅方法估算π值的完整示例,包括可視化。這個示例展示了如何利用 NumPy 的高效數組操作來生成大量隨機點,計算它們是否落在單位圓內,並最終估算出π值。同時,我們還使用 Matplotlib 來視覺化這個過程,直觀地展示了蒙特卡羅方法的工作原理。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib

# 使用默認英文字體並允許負號顯示
matplotlib.rcParams['font.family'] = 'sans-serif'
matplotlib.rcParams['axes.unicode_minus'] = False

# 模擬點的數量
n = 1000000

# 設置隨機種子以確保結果可重現
np.random.seed(42)

# 在 [-1, 1] x [-1, 1] 範圍內生成隨機點
points = np.random.uniform(-1, 1, (n, 2))

# 計算每個點到原點的距離
distances = np.sqrt(points[:, 0]**2 + points[:, 1]**2)

# 判斷點是否在單位圓內
inside_circle = distances <= 1

# 估算π值
pi_estimate = 4 * np.sum(inside_circle) / n

# 可視化(為了清晰起見,只繪製前 10000 個點)
plt.figure(figsize=(8, 8))
plt.scatter(points[inside_circle][:10000, 0], points[inside_circle][:10000, 1], 
            color='green', s=1, label='Inside Circle')
plt.scatter(points[~inside_circle][:10000, 0], points[~inside_circle][:10000, 1], 
            color='red', s=1, label='Outside Circle')

# 繪製單位圓邊界
theta = np.linspace(0, 2 * np.pi, 200)
plt.plot(np.cos(theta), np.sin(theta), color='green')

# 顯示標題和圖例
plt.title(f'Monte Carlo Estimation of π: {pi_estimate:.8f}, Actual π: {np.pi:.8f}')
plt.axis('equal')
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.legend()
plt.grid(True)
plt.show()

這段代碼不僅計算了π的估計值,還通過散點圖直觀地展示了蒙特卡羅方法的原理。

NumPy 與蒙特卡羅方法的優勢

計算效率

NumPy 的向量化操作極大地提升了蒙特卡羅模擬的效率。透過利用底層優化和並行處理能力,NumPy 能夠同時處理大量數據點,而無需使用耗時的 Python 循環。這種方法特別適合蒙特卡羅模擬,因為它們通常需要處理大量的隨機樣本。

相較於傳統的純 Python 實現,NumPy 的向量化操作可以將計算速度提高數十倍,甚至更多。這種顯著的性能提升使得更大規模和更複雜的蒙特卡羅模擬成為可能,從而提高了估算的精度和可靠性。對於需要進行大量迭代或處理海量數據的應用場景,NumPy 的高效計算能力尤為重要。

代碼簡潔性

使用 NumPy 進行蒙特卡羅模擬,可以用極少的代碼實現複雜的計算。例如,估算π值的核心代碼只需要幾行:

points = np.random.rand(n, 2)
inside_circle = np.sqrt(points[:, 0]**2 + points[:, 1]**2) < 1
pi_estimate = 4 * np.sum(inside_circle) / n

精度控制與提升

透過增加抽樣點數,蒙特卡羅方法的精度可以持續提高。根據大數定律,估計值會隨著樣本數的增加而更接近真實值。

使用 NumPy 還可以實現更複雜的隨機數生成策略,例如使用不同的隨機種子來減少隨機數序列之間的相關性:

# 使用兩個獨立的 RandomState 以減少 x 和 y 座標間的相關性
rng_x = np.random.RandomState(42)
rng_y = np.random.RandomState(123)

x = rng_x.uniform(-1, 1, n)
y = rng_y.uniform(-1, 1, n)

inside_circle = (x**2 + y**2) <= 1
pi_estimate = 4 * np.sum(inside_circle) / n

這種方法可以提高估計的準確性,特別是對於較小樣本量的情況。

結論

NumPy 與蒙特卡羅方法的結合為解決複雜數值問題提供了一個強大而靈活的工具。在估算π值的案例中,我們可以清晰地看到這種方法的原理、實現和優勢。NumPy 的向量化操作大大提高了計算效率,使得即使在普通的個人電腦上也能進行大規模的模擬計算。

蒙特卡羅方法的優勢在於其簡單性和通用性,適用於許多難以用解析方法解決的問題。而 NumPy 則通過其高效的數組操作和豐富的數學函數庫,為實現這些方法提供了堅實的基礎。

隨著計算能力的不斷提升和算法的改進,NumPy 與蒙特卡羅方法的結合將在科學研究、工程應用、金融分析等眾多領域發揮越來越重要的作用。對於程式開發者和研究人員來說,掌握這一工具組合無疑是一項寶貴的技能。

感謝您閱讀我的文章。歡迎隨時分享你的想法。
關於 Calpa

Calpa 擅長使用 TypeScriptReact.jsVue.js 建立 Responsive Website。

他積極參與開源社區,曾在 2019 年的香港開源大會上擔任講者,提供工作經驗和見解。此外,他也在 GitHub 上公開分享個人博客程式碼,已獲得超過 300 顆星星和 60 個分支的支持。

他熱愛學習新技術,並樂意分享經驗。他相信,唯有不斷學習才能跟上快速演變的技術環境。