蒙特卡洛方法（Monte Carlo method）是一类通过随机采样来近似求解确定性问题的数值计算方法。估算圆周率 是蒙特卡洛方法最经典的教学案例之一，直观展示了"用随机性解决确定性问题的核心思想。
考虑一个单位正方形 ，其内切圆为半径为 1 的圆，圆心在原点。
因此，圆的面积与正方形面积之比为：
如果我们在正方形内随机均匀地投掷点，那么落在圆内的点的比例应当近似等于 。因此：
其中 是总投掷点数， 是落在圆内的点数。
对于随机点 ，其中 （均匀分布），点落在圆内的条件为：
假设我们投掷 20 个随机点。使用均匀分布生成 和 坐标：
| 点编号 | x | y | 在圆内？ | |
|---|---|---|---|---|
| 1 | 0.32 | -0.48 | 0.3328 | ✅ |
| 2 | -0.71 | 0.63 | 0.9010 | ✅ |
| 3 | 0.85 | 0.19 | 0.7586 | ✅ |
| 4 | -0.12 | -0.92 | 0.8608 | ✅ |
| 5 | 0.55 | -0.74 | 0.8501 | ✅ |
| 6 | -0.88 | -0.45 | 0.9769 | ✅ |
| 7 | 0.03 | 0.99 | 0.9810 | ✅ |
| 8 | 0.67 | 0.72 | 0.9673 | ✅ |
| 9 | -0.34 | 0.56 | 0.4292 | ✅ |
| 10 | 0.91 | 0.38 | 0.9725 | ✅ |
| 11 | -0.25 | -0.83 | 0.7514 | ✅ |
| 12 | 0.48 | 0.88 | 1.0048 | ❌（略大于1） |
| 13 | -0.65 | -0.23 | 0.4754 | ✅ |
| 14 | 0.73 | -0.67 | 0.9818 | ✅ |
| 15 | -0.55 | -0.55 | 0.6050 | ✅ |
| 16 | 0.89 | 0.45 | 0.9946 | ✅ |
| 17 | -0.41 | 0.29 | 0.2522 | ✅ |
| 18 | 0.12 | -0.96 | 0.9360 | ✅ |
| 19 | -0.78 | -0.62 | 0.9928 | ✅ |
| 20 | -0.05 | 0.15 | 0.0250 | ✅ |
结果： 19/20 点在圆内。
这个结果误差很大（真实值 ），因为样本太少。随着样本增加，结果会收敛。
| 投掷点数 | 圆内点数 | 估计值 | 误差百分比 |
|---|---|---|---|
| 100 | 79 | 3.16000 | 0.59% |
| 1,000 | 788 | 3.15200 | 0.33% |
| 10,000 | 7,848 | 3.13920 | -0.08% |
| 100,000 | 78,536 | 3.14144 | -0.005% |
| 1,000,000 | 785,354 | 3.14142 | -0.005% |
| 10,000,000 | 7,853,930 | 3.14157 | -0.0005% |
注意：由于随机性，每次运行结果会有差异。但 的收敛速度保证了误差随 增大而减小。
┌─────────────────────────────────┐
│ 开始 │
└─────────────────────────────────┘
│
▼
┌─────────────────────────────────┐
│ 设置总采样数 N │
│ 初始化计数器 count = 0 │
└─────────────────────────────────┘
│
▼
┌─────────────────────────────────┐
│ for i = 1 to N: │
│ x = random(-1, 1) │
│ y = random(-1, 1) │
└─────────────────────────────────┘
│
▼
┌─────────────────────────────────┐
│ 检查 x² + y² ≤ 1 │
│ 是 → count += 1 │
└─────────────────────────────────┘
│
▼
┌─────────────────────────────────┐
│ π ≈ 4 × count / N │
│ 输出结果 │
└─────────────────────────────────┘
蒙特卡洛方法的理论基础是大数定律（Law of Large Numbers）。设 为指示变量：
则 。
根据大数定律：
因此， 是 的一致估计量。
蒙特卡洛方法的收敛速度为 。这意味着：
误差分析：
估计值的标准差为：
代入 ：
| 样本量 N | 标准差 | 95% 置信区间 |
|---|---|---|
| 100 | 0.164 | |
| 1,000 | 0.052 | |
| 10,000 | 0.016 | |
| 100,000 | 0.005 | |
| 1,000,000 | 0.002 |
根据中心极限定理（Central Limit Theorem），当 足够大时，估计值 近似服从正态分布：
其中 。
import random
import math
import matplotlib.pyplot as plt
def estimate_pi(N):
"""使用蒙特卡洛方法估算圆周率。"""
count = 0
for _ in range(N):
x = random.uniform(-1, 1)
y = random.uniform(-1, 1)
if x*x + y*y <= 1:
count += 1
return 4.0 * count / N
# 测试不同精度
for N in [100, 1000, 10000, 100000, 1000000]:
pi_est = estimate_pi(N)
error = abs(pi_est - math.pi)
print(f"N={N:7d}: π≈{pi_est:.6f}, 误差={error:.6f}")
输出示例：
N= 100: π≈3.200000, 误差=0.058407
N= 1000: π≈3.140000, 误差=0.001593
N= 10000: π≈3.147200, 误差=0.005607
N= 100000: π≈3.142120, 误差=0.000527
N=1000000: π≈3.141704, 误差=0.000111
import numpy as np
def estimate_pi_vectorized(N):
"""向量化实现，比循环快 100 倍以上。"""
x = np.random.uniform(-1, 1, N)
y = np.random.uniform(-1, 1, N)
inside = (x*x + y*y) <= 1.0
return 4.0 * np.sum(inside) / N
# 性能对比
import time
N = 10000000 # 1000 万样本
t0 = time.time()
# pi_loop = estimate_pi(N) # 太慢，仅做基准
t1 = time.time()
t0 = time.time()
pi_vec = estimate_pi_vectorized(N)
t1 = time.time()
print(f"π ≈ {pi_vec:.8f}, 真实值 = {math.pi:.8f}")
print(f"误差 = {abs(pi_vec - math.pi):.2e}")
print(f"耗时: {t1-t0:.2f} 秒")
输出示例：
π ≈ 3.14154560, 真实值 = 3.14159265
误差 = 4.70e-05
耗时: 0.24 秒
from multiprocessing import Pool
import numpy as np
def monte_carlo_chunk(chunk_size):
"""计算单个 chunk 的圆内点数。"""
x = np.random.uniform(-1, 1, chunk_size)
y = np.random.uniform(-1, 1, chunk_size)
return np.sum(x*x + y*y <= 1.0)
def estimate_pi_parallel(total_N, num_workers=8):
"""多进程并行蒙特卡洛。"""
chunk_size = total_N // num_workers
with Pool(num_workers) as pool:
counts = pool.map(monte_carlo_chunk,
[chunk_size] * num_workers)
return 4.0 * sum(counts) / total_N
# 测试：1 亿样本
pi_parallel = estimate_pi_parallel(100_000_000, num_workers=8)
print(f"π ≈ {pi_parallel:.8f}")
print(f"误差 = {abs(pi_parallel - math.pi):.2e}")
将单位正方形划分为子区域，在每个子区域内均匀采样，保证采样空间覆盖更均匀。
def estimate_pi_stratified(N, grid=10):
"""分层采样：将正方形分为 grid×grid 个格子，每个格子内采样。"""
points_per_cell = N // (grid * grid)
count = 0
for i in range(grid):
for j in range(grid):
# 每个格子内的局部坐标
x_base = -1 + 2 * i / grid
y_base = -1 + 2 * j / grid
cell_size = 2.0 / grid
for _ in range(points_per_cell):
x = x_base + random.random() * cell_size
y = y_base + random.random() * cell_size
if x*x + y*y <= 1:
count += 1
return 4.0 * count / N
对比实验（N=10,000）：
| 方法 | 估计值 | 误差 | 标准差 |
|---|---|---|---|
| 朴素蒙特卡洛 | 3.13920 | -0.00239 | 0.0164 |
| 分层采样 (5×5) | 3.14231 | -0.00128 | 0.0102 |
| 分层采样 (10×10) | 3.14187 | -0.00028 | 0.0056 |
| 分层采样 (20×20) | 3.14156 | -0.00003 | 0.0028 |
分层采样可以将方差降低 50%-80%，大幅提升收敛速度。
如果知道哪部分区域"更重要"，可以增加在该区域的采样密度，然后通过权重修正偏差。
对于 的估算，由于均匀分布本身是 optimal 的（圆在正方形中均匀随机），重要采样的收益有限，此处仅做概念介绍。
使用低差异序列（Low-Discrepancy Sequence，如 Sobol 序列或 Halton 序列）代替纯随机数，可以获得更快的收敛速度 。
from scipy.stats.qmc import Halton
def estimate_pi_qmc(N):
"""使用 Halton 低差异序列。"""
sampler = Halton(d=2, scramble=False)
samples = sampler.random(N)
# 将 [0,1] 映射到 [-1,1]
x = 2 * samples[:, 0] - 1
y = 2 * samples[:, 1] - 1
inside = (x*x + y*y) <= 1.0
return 4.0 * np.sum(inside) / N
方法对比（N=1,000）：
| 方法 | 10 次运行平均误差 | 95% 置信区间宽度 |
|---|---|---|
| 朴素蒙特卡洛 | 0.0523 | 0.205 |
| 分层采样 (10×10) | 0.0315 | 0.098 |
| 拟蒙特卡洛 (Halton) | 0.0187 | 0.042 |
| 方法 | 收敛速度 | 实现复杂度 | 每步获得位数 | 适合场景 |
|---|---|---|---|---|
| 蒙特卡洛 | ⭐ 最低 | ~0.5位/10000步 | 教学、并行计算演示 | |
| 莱布尼茨级数 | ⭐⭐ 低 | ~0.3位/100步 | 数学推导 | |
| 高斯-勒让德算法 | ⭐⭐⭐⭐ 高 | 2倍位数/步 | 高精度计算 | |
| 丘德诺夫斯基算法 | ⭐⭐⭐⭐⭐ 很高 | ~14位/项 | 世界纪录 | |
| BBP 公式 | ⭐⭐⭐ 中 | ~1位/项 | 十六进制任意位 |
表示 位大数乘法的复杂度。
收敛速度可视化：
| 目标精度 | 蒙特卡洛所需样本 | 高斯-勒让德所需迭代 |
|---|---|---|
| 1 位小数 (误差 < 0.05) | 2 次 | |
| 3 位小数 (误差 < 0.0005) | 3 次 | |
| 6 位小数 (误差 < 5\times 10^{-7}) | 4 次 | |
| 10 位小数 (误差 < 5\times 10^{-11}) | 5 次 |
蒙特卡洛方法虽然收敛慢，但有其独特优势：实现简单、天然可并行、可以处理高维问题（如高维积分），而这些是确定性方法难以做到的。
估算 只是蒙特卡洛方法的一个教学案例。在实际工程和科研中，蒙特卡洛方法有广泛的应用：
对于 维积分 ：
这是蒙特卡洛方法在高维问题中无可替代的核心优势。
def monte_carlo_option_pricing(S0, K, r, sigma, T, N=100000):
"""使用蒙特卡洛模拟欧式看涨期权价格。"""
Z = np.random.standard_normal(N)
ST = S0 * np.exp((r - 0.5*sigma**2)*T + sigma*np.sqrt(T)*Z)
payoff = np.maximum(ST - K, 0)
return np.exp(-r*T) * np.mean(payoff)
# 示例：S0=100, K=105, r=0.05, σ=0.2, T=1
price = monte_carlo_option_pricing(100, 105, 0.05, 0.2, 1)
print(f"欧式看涨期权价格: ${price:.2f}")
在核反应堆设计中，蒙特卡洛方法用于模拟中子在介质中的随机行走：
这种"模拟物理过程"的方法被称为"直接模拟蒙特卡洛"（Direct Simulation Monte Carlo, DSMC）。
马尔可夫链蒙特卡洛（Markov Chain Monte Carlo, MCMC）是贝叶斯推断的核心工具：
import pymc3 as pm
# 贝叶斯线性回归的 MCMC 推断
with pm.Model():
alpha = pm.Normal('alpha', mu=0, sigma=10)
beta = pm.Normal('beta', mu=0, sigma=10)
sigma = pm.HalfNormal('sigma', 5)
mu = alpha + beta * x
y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=y)
trace = pm.sample(2000, tune=1000) # MCMC 采样
| 生成器 | 周期 | 速度（百万/秒） | 质量 |
|---|---|---|---|
Python random.random() (MT19937) |
~5 | 足够 | |
NumPy np.random.uniform (PCG64) |
~200 | 优秀 | |
secrets.randbits |
N/A | ~0.5 | 密码学安全 |
| Sobol 序列 | 确定性 | ~100 | 最佳收敛 |
推荐：使用 NumPy 的 PCG64 生成器，兼顾速度和统计质量。
| 技术 | 原理 | 方差降低倍数 |
|---|---|---|
| 分层采样 | 确保空间覆盖均匀 | 2-5× |
| 对偶变量（Antithetic Variables） | 使用负相关样本 | 2× |
| 控制变量（Control Variates） | 利用已知量修正 | 3-10× |
| 重要性采样 | 聚焦重要区域 | 10-100× |
rand() 的周期仅 ）总结： 蒙特卡洛方法估算 是最简单直观的并行计算演示之一。虽然对于计算 本身效率远不如解析算法，但它的通用性、可扩展性和对高维问题的适应性使其成为计算科学中不可或缺的工具。理解蒙特卡洛方法的原理、收敛性和方差控制技巧，是掌握现代数值计算的基础。