概率方法与计算机模拟

邱飞

2026-01-01

引言：从经典问题到现代方法

为什么要学习概率方法?

概率方法是现代科学计算的核心工具之一:

理论价值: 将确定性问题转化为概率模型
计算优势: 通过随机模拟逼近复杂的解析解
广泛应用: 从物理模拟到金融工程,从机器学习到运筹优化

本讲座目标: 通过经典的Buffon投针实验,展示如何用概率方法和计算机模拟解决实际问题。

案例研究：估算圆周率

Buffon投针实验简介

18世纪法国博物学家Buffon提出的经典概率问题:

核心思想: 通过几何概率模型,将物理实验与数学常数\(\pi\)联系起来。

现代意义: 这是蒙特卡洛方法的早期原型——用随机抽样估计确定性量。

几何概型定义

当随机试验的样本空间是某个可度量的区域（长度、面积或体积），且样本点落入区域内各处的可能性相等时：

\[ P(A) = \frac{S_A}{S} = \frac{\text{构成事件 A 的子区域度量}}{\text{样本空间的总度量}} \tag{1}\]

Note

Understanding Geometric Probability

Geometric probability differs from classical probability in that the sample space is continuous rather than discrete. The key assumption is uniform distribution — every point in the space has equal likelihood. This is analogous to throwing a dart at a board: the probability of hitting region A depends solely on its area relative to the total board area.

本讲座将复现经典的 Buffon (蒲丰) 投针问题,展示如何通过几何概率模型和计算机模拟来估算圆周率 \(\pi\)。

物理模型构建

实验设定：在一个平面上画有间距为 \(a\) 的平行线，随机投掷长度为 \(l\) (\(l < a\)) 的针。

关键变量定义：

\(x\): 针的中心点 \(M\) 到最近一条平行线的垂直距离,取值范围 \([0, a/2]\)
\(\theta\): 针与平行线的夹角,取值范围 \([0, \pi]\)

Tip

Why These Variables?

We only need to consider \(x \in [0, a/2]\) due to symmetry — the probability is the same whether the needle is closer to the line above or below. Similarly, \(\theta \in [0, \pi]\) covers all possible orientations since angles beyond \(\pi\) are equivalent due to the needle’s symmetry.

Figure 1: Buffon投针实验的物理设置图示

数学建模

为了利用 Equation 1 中的几何概型公式,我们需要将物理问题转化为平面区域问题。

样本空间构建

由于 \(\theta\) 和 \(x\) 都是连续型随机变量且相互独立,样本空间 \(S\) 是一个矩形区域:

\[ S = \{(\theta, x) : 0 \le \theta \le \pi, \; 0 \le x \le a/2\} \]

样本空间的总度量(面积)为: \(|S| = \pi \cdot \frac{a}{2} = \frac{\pi a}{2}\)

相交事件的几何条件

事件 \(A\): 针与某条平行线相交

几何分析: 如 Figure 1 所示,针在垂直方向的投影长度为 \(\frac{l}{2}\sin\theta\)。当针中心到最近直线的距离 \(x\) 小于这个投影时,针必然与直线相交。

相交条件: \[ x \le \frac{l}{2} \sin \theta \tag{2}\]

Note

Geometric Intuition

Think of the needle as a rotating line segment. At angle \(\theta\), its ‘reach’ toward the nearest line is \(\frac{l}{2}\sin\theta\) (half the needle’s length times the sine component). If the center is within this reach distance from the line, intersection occurs.

在 \((\theta, x)\) 坐标系中, Equation 2 对应于正弦曲线下方的区域。

理论可视化：相交区域

下图展示了样本空间 \(S\) 和事件 \(A\) 的几何关系。阴影区域的面积即为相交概率的几何度量。

Code

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

# 实验参数设置
l = 0.8  # 针长
a = 1.0  # 平行线间距

# 定义样本空间: theta in [0, pi], x in [0, a/2]
theta_range = np.linspace(0, np.pi, 500)
x_boundary = (l / 2) * np.sin(theta_range)  # 相交条件的临界曲线

# 创建图形
fig, ax = plt.subplots(figsize=(11, 6))

# 绘制样本空间的矩形边界
rect_x = [0, np.pi, np.pi, 0, 0]
rect_y = [0, 0, a/2, a/2, 0]
ax.plot(rect_x, rect_y, 'k-', linewidth=2.5, label='样本空间边界')

# 绘制相交条件的临界曲线
ax.plot(theta_range, x_boundary, 'r-', linewidth=3, 
        label=r'相交边界: $x = \frac{l}{2}\sin\theta$')

# 填充相交区域 (事件A)
ax.fill_between(theta_range, 0, x_boundary, 
                color='#E3120B', alpha=0.25, 
                label=r'事件$A$: 针与线相交')

# 坐标轴设置
ax.set_xlim(-0.1, np.pi + 0.1)
ax.set_ylim(-0.05, a/2 + 0.1)
ax.set_xticks([0, np.pi/4, np.pi/2, 3*np.pi/4, np.pi])
ax.set_xticklabels(['0', r'$\pi/4$', r'$\pi/2$', r'$3\pi/4$', r'$\pi$'], fontsize=14)
ax.set_yticks([0, a/4, a/2])
ax.set_yticklabels(['0', r'$a/4$', r'$a/2$'], fontsize=14)
ax.set_xlabel(r'针的夹角 $\theta$ (弧度)', fontsize=16, fontweight='bold')
ax.set_ylabel(r'针心到直线的距离 $x$', fontsize=16, fontweight='bold')
ax.set_title('Buffon投针问题的几何概率模型', fontsize=18, fontweight='bold', pad=20)

# 网格和图例
ax.grid(True, alpha=0.3, linestyle='--')
ax.legend(loc='upper right', fontsize=13, framealpha=0.95)

# 添加关键点标注
ax.annotate('最大投影点', xy=(np.pi/2, l/2), xytext=(np.pi/2 + 0.5, l/2 + 0.05),
            arrowprops=dict(arrowstyle='->', color='#2C3E50', lw=2),
            fontsize=12, color='#2C3E50', fontweight='bold')

plt.tight_layout()
plt.show()

关键观察: 当 \(\theta = \pi/2\) (针垂直于平行线)时,投影达到最大值 \(l/2\),相交概率最高。

π估计公式的完整推导

步骤1: 计算相交区域面积

根据 Equation 2,事件 \(A\) 对应的区域为正弦曲线下方的面积:

\[ S_A = \int_0^{\pi} \frac{l}{2}\sin\theta \, d\theta \]

Tip

Integration Details

Recall that \(\int \sin\theta \, d\theta = -\cos\theta + C\). Applying the fundamental theorem of calculus:

\[ S_A = \frac{l}{2} \left[-\cos\theta\right]_0^{\pi} = \frac{l}{2}\left[(-\cos\pi) - (-\cos 0)\right] = \frac{l}{2}[1 - (-1)] = l \]

因此,相交区域的面积为: \(S_A = l\)

步骤2: 推导理论概率

样本空间的总面积为 \(|S| = \frac{\pi a}{2}\),根据 Equation 1:

\[ P(A) = \frac{S_A}{|S|} = \frac{l}{\frac{\pi a}{2}} = \frac{2l}{\pi a} \tag{3}\]

步骤3: 从频率到概率的桥梁

根据大数定律,当实验次数 \(n\) 足够大时,事件 \(A\) 的频率 \(\frac{k}{n}\) 收敛于其理论概率:

\[ \frac{k}{n} \xrightarrow{n \to \infty} P(A) = \frac{2l}{\pi a} \tag{4}\]

步骤4: 反解圆周率

从 Equation 4 中解出 \(\pi\):

\[ \pi \approx \frac{2ln}{ak} \tag{5}\]

Caution

Estimation Accuracy

The approximation quality in Equation 5 depends critically on \(n\) (sample size) and \(k\) (intersection count). The standard error decreases as \(O(1/\sqrt{n})\). For accurate \(\pi\) estimation, we typically need \(n \geq 10^4\) trials.

结论: 通过随机投针实验,我们建立了物理过程(投针)与数学常数(\(\pi\))之间的桥梁。

蒙特卡洛模拟实现

算法原理

蒙特卡洛方法的核心思想:通过随机抽样估计确定性量。

具体步骤:

在样本空间 \((\theta, x)\) 矩形区域内随机生成 \(n\) 个点
检验每个点是否满足 Equation 2 中的相交条件
统计相交次数 \(k\),利用 Equation 5 估算 \(\pi\)

Python代码实现

Listing 1

Code

import numpy as np

def buffon_simulation(n, l, a):
    """
    执行Buffon投针实验的蒙特卡洛模拟。
    
    Parameters
    ----------
    n : int
        模拟次数(投针次数)
    l : float
        针的长度
    a : float
        平行线之间的间距(必须满足 l < a)
    
    Returns
    -------
    pi_estimate : float
        根据公式估算的圆周率值
    theta_samples : ndarray
        所有针的角度样本 (shape: (n,))
    x_samples : ndarray
        所有针心到直线距离的样本 (shape: (n,))
    intersection_mask : ndarray
        布尔数组,标记每根针是否与线相交 (shape: (n,))
    """
    # 步骤1: 生成服从均匀分布的随机样本
    # theta ~ U(0, π): 针与平行线的夹角
    theta_samples = np.random.uniform(0, np.pi, n)
    
    # x ~ U(0, a/2): 针心到最近直线的距离
    x_samples = np.random.uniform(0, a/2, n)
    
    # 步骤2: 计算相交条件的临界值
    # 针在垂直方向的投影长度
    critical_distance = (l / 2) * np.sin(theta_samples)
    
    # 步骤3: 判断是否相交
    # 当实际距离x小于等于临界投影时,针与线相交
    intersection_mask = x_samples <= critical_distance
    
    # 步骤4: 统计相交次数
    k = np.sum(intersection_mask)
    
    # 步骤5: 根据公式估算π
    # 防止除零错误:当k=0时返回NaN
    pi_estimate = (2 * l * n) / (a * k) if k > 0 else np.nan
    
    return pi_estimate, theta_samples, x_samples, intersection_mask

# 执行模拟实验
np.random.seed(42)  # 设置随机种子以保证可重复性
n_trials = 20000    # 模拟次数
needle_length = 0.8 # 针长
line_spacing = 1.0  # 平行线间距

# 调用模拟函数
pi_estimated, theta_vals, x_vals, hit_mask = buffon_simulation(
    n_trials, needle_length, line_spacing
)

# 输出结果
print(f'模拟次数: {n_trials:,}')
print(f'相交次数: {np.sum(hit_mask):,}')
print(f'相交频率: {np.sum(hit_mask)/n_trials:.4f}')
print(f'估算的π值: {pi_estimated:.6f}')
print(f'真实π值: {np.pi:.6f}')
print(f'相对误差: {abs(pi_estimated - np.pi)/np.pi * 100:.2f}%')

模拟次数: 20,000
相交次数: 10,263
相交频率: 0.5131
估算的π值: 3.117997
真实π值: 3.141593
相对误差: 0.75%

模拟结果可视化

下图展示了 Listing 1 中模拟生成的20,000个样本点在 \((\theta, x)\) 空间中的分布。

颜色编码:

红色点: 针与平行线相交(位于正弦曲线下方)
青色点: 针未与平行线相交(位于正弦曲线上方)

Code

# 创建高质量可视化图形
fig, ax = plt.subplots(figsize=(12, 6))

# 绘制理论相交边界(正弦曲线)
theta_theory = np.linspace(0, np.pi, 300)
boundary_theory = (needle_length / 2) * np.sin(theta_theory)
ax.plot(theta_theory, boundary_theory, 'k-', linewidth=3, 
        label=r'理论边界: $x = \frac{l}{2}\sin\theta$', zorder=5)

# 绘制未相交的样本点(背景层)
ax.scatter(theta_vals[~hit_mask], x_vals[~hit_mask], 
           color='#008080', s=8, alpha=0.4, 
           label=f'未相交: {np.sum(~hit_mask):,} 次', zorder=2)

# 绘制相交的样本点(前景层)
ax.scatter(theta_vals[hit_mask], x_vals[hit_mask], 
           color='#E3120B', s=10, alpha=0.6, 
           label=f'相交: {np.sum(hit_mask):,} 次', zorder=3)

# 坐标轴设置
ax.set_xlim(-0.1, np.pi + 0.1)
ax.set_ylim(-0.02, line_spacing/2 + 0.05)
ax.set_xticks([0, np.pi/4, np.pi/2, 3*np.pi/4, np.pi])
ax.set_xticklabels(['0', r'$\pi/4$', r'$\pi/2$', r'$3\pi/4$', r'$\pi$'], fontsize=13)
ax.set_yticks([0, line_spacing/4, line_spacing/2])
ax.set_yticklabels(['0', r'$a/4$', r'$a/2$'], fontsize=13)
ax.set_xlabel(r'针的夹角 $\theta$ (弧度)', fontsize=15, fontweight='bold')
ax.set_ylabel(r'针心到直线的距离 $x$', fontsize=15, fontweight='bold')

# 标题包含估算结果
ax.set_title(
    f'蒙特卡洛模拟结果 (n={n_trials:,})\n'
    f'估算值: π ≈ {pi_estimated:.6f} | 真实值: π = {np.pi:.6f} | '
    f'误差: {abs(pi_estimated - np.pi):.6f} ({abs(pi_estimated - np.pi)/np.pi * 100:.2f}%)',
    fontsize=14, fontweight='bold', pad=15
)

# 网格和图例
ax.grid(True, alpha=0.3, linestyle='--', linewidth=0.5)
ax.legend(loc='upper right', fontsize=12, framealpha=0.95, edgecolor='black')

# 边框设置
for spine in ax.spines.values():
    spine.set_linewidth(1.5)
    spine.set_color('#2C3E50')

plt.tight_layout()
plt.show()

观察结果: 红色点密集分布在正弦曲线下方,验证了理论模型的正确性。相交频率 \(k/n\) 与理论概率 Equation 3 高度吻合。

收敛性分析

让我们通过增加模拟次数来观察估算精度的提升:

Code

# 进行多次不同规模的模拟
sample_sizes = np.logspace(2, 5, 50, dtype=int)  # 从100到100,000
pi_estimates = []

np.random.seed(42)
for n in sample_sizes:
    pi_est, _, _, _ = buffon_simulation(n, needle_length, line_spacing)
    pi_estimates.append(pi_est)

pi_estimates = np.array(pi_estimates)

# 绘制收敛图
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# 左图: π估算值随样本量的变化
ax1.semilogx(sample_sizes, pi_estimates, 'o-', color='#E3120B', 
             markersize=4, linewidth=1.5, alpha=0.7, label='蒙特卡洛估算值')
ax1.axhline(y=np.pi, color='#2C3E50', linestyle='--', linewidth=2.5, label=r'真实值: $\pi$')
ax1.fill_between(sample_sizes, np.pi - 0.1, np.pi + 0.1, 
                  color='#2C3E50', alpha=0.1, label='±0.1 容差带')
ax1.set_xlabel('模拟次数 (对数尺度)', fontsize=13, fontweight='bold')
ax1.set_ylabel(r'估算的 $\pi$ 值', fontsize=13, fontweight='bold')
ax1.set_title('π估算值的收敛性', fontsize=14, fontweight='bold')
ax1.grid(True, alpha=0.3, which='both', linestyle=':')
ax1.legend(fontsize=11)

# 右图: 绝对误差随样本量的变化
errors = np.abs(pi_estimates - np.pi)
ax2.loglog(sample_sizes, errors, 's-', color='#008080', 
           markersize=4, linewidth=1.5, alpha=0.7, label='绝对误差')
# 添加理论O(1/√n)参考线
theoretical_error = 2 / np.sqrt(sample_sizes)
ax2.loglog(sample_sizes, theoretical_error, '--', color='#F0A700', 
           linewidth=2, label=r'理论: $O(1/\sqrt{n})$')
ax2.set_xlabel('模拟次数 (对数尺度)', fontsize=13, fontweight='bold')
ax2.set_ylabel('绝对误差 (对数尺度)', fontsize=13, fontweight='bold')
ax2.set_title('误差收敛速率', fontsize=14, fontweight='bold')
ax2.grid(True, alpha=0.3, which='both', linestyle=':')
ax2.legend(fontsize=11)

plt.tight_layout()
plt.show()

Note

Convergence Rate Analysis

The Monte Carlo method has a convergence rate of \(O(1/\sqrt{n})\), which is independent of the problem’s dimensionality. This is both a strength (works well in high dimensions) and a limitation (requires 100× more samples to improve accuracy by 10×).

总结

核心成果

通过本案例,我们完成了从理论建模到计算实现的完整流程:

几何建模: 将物理投针问题抽象为二维参数空间 \((\theta, x)\) (见 Figure 1)
概率分析: 推导相交条件 Equation 2 和理论概率 Equation 3
算法设计: 实现蒙特卡洛模拟 (Listing 1)
数值求解: 通过 Equation 5 反解圆周率,误差控制在1%以内

方法论意义

蒙特卡洛方法的核心优势:

通用性: 适用于难以解析求解的复杂问题
可扩展性: 易于并行化,适合大规模计算
维度鲁棒性: 收敛速度不随问题维度显著恶化

实际应用领域

金融工程: 期权定价、风险估计(VaR)、投资组合优化
物理模拟: 粒子输运、量子蒙特卡洛、分子动力学
机器学习: 贝叶斯推断(MCMC)、强化学习(蒙特卡洛树搜索)
运筹优化: 排队论、库存管理、交通流模拟

Tip

Further Reading

For deeper understanding of Monte Carlo methods: - Variance Reduction: Importance sampling, control variates - Advanced Sampling: Markov Chain Monte Carlo (MCMC), Gibbs sampling - Quasi-Monte Carlo: Low-discrepancy sequences for better convergence

感谢聆听! 欢迎提问与讨论。