蒙特卡洛方法估计圆周率π的收敛速度与方差缩减技术

蒙特卡洛方法估计圆周率π的收敛速度与方差缩减技术

1. 理解蒙特卡洛方法的核心思想

蒙特卡洛方法利用随机抽样和概率统计来求解数学问题。用它估计圆周率π的思路基于几何概率。

想象一个边长为2的正方形，内切一个半径为1的圆形。正方形的面积是4，圆形的面积是π。那么，随机在正方形内均匀地撒点，点落在圆形内的概率就等于圆面积除以正方形面积，即 π/4。

因此，通过统计大量随机点中落在圆内的比例，乘以4，就能得到π的估计值。点撒得越多，这个估计就越准。

2. 实现最基础的π估计

编写一段简单的程序来模拟这个过程。这里以Python为例。

import random
import math

def estimate_pi_basic(num_samples):
    """使用基础蒙特卡洛方法估计π。"""
    points_in_circle = 0
    for _ in range(num_samples):
        # 生成一个在正方形[-1,1] x [-1,1]内的均匀随机点
        x = random.uniform(-1.0, 1.0)
        y = random.uniform(-1.0, 1.0)
        # 判断点是否在单位圆内
        if math.sqrt(x**2 + y**2) <= 1.0:
            points_in_circle += 1
    # 计算π的估计值
    pi_estimate = 4 * points_in_circle / num_samples
    return pi_estimate

# 运行一次示例
print(f"估计值 (N=100000): {estimate_pi_basic(100000)}")

运行此代码，你会得到一个接近3.14的值。增加 num_samples 的值，结果会更接近真实的π。

3. 分析收敛速度：误差如何随样本量变化

基础方法的“收敛速度”描述了估计误差（估计值与真实π的差值）随着样本量 N 增大而减小的快慢。

根据概率论中的中心极限定理，基础蒙特卡洛估计量的标准误差与 $1/\sqrt{N}$ 成正比。这意味着：

• 要想将误差减小一半，你需要将样本量增加到原来的 4 倍。
• 要想将误差减小到十分之一，你需要将样本量增加到原来的 100 倍。

误差的这种 $1/\sqrt{N}$ 收敛速度是基础蒙特卡洛方法的一个主要缺点：为了获得高精度，需要消耗巨大的计算成本。

4. 引入第一种方差缩减技术：对偶变量法

方差缩减技术的目标是，在不增加样本量（或轻微增加计算量）的前提下，显著降低估计量的方差，从而加速收敛。

对偶变量法的核心是引入负相关的样本对。对于估计π，一个自然的对偶是：如果有一个随机点 (x, y)，那么 (-x, y) 就是其对偶点。

修改估算函数：

import numpy as np

def estimate_pi_antithetic(num_samples):
    """使用对偶变量法估计π。"""
    points_in_circle = 0
    for _ in range(num_samples // 2):  # 注意：只需生成一半的样本对
        # 生成一个随机角度和半径
        theta = np.random.uniform(0, 2 * np.pi)
        r = np.sqrt(np.random.uniform(0, 1))  # 保证在单位圆内均匀分布

        # 第一个点
        x1 = r * np.cos(theta)
        y1 = r * np.sin(theta)

        # 其对偶点：角度旋转180度（即π）
        x2 = -x1
        y2 = -y1  # 或者使用 y2 = -y1，根据分布定义调整

        # 这里使用更简单的几何对偶：关于y轴对称
        # 我们重新使用基于正方形内点的方法来演示对偶变量法
        x = np.random.uniform(-1.0, 1.0)
        y = np.random.uniform(-1.0, 1.0)

        # 判断第一个点
        if x**2 + y**2 <= 1.0:
            points_in_circle += 1

        # 判断其对偶点 (-x, y)
        if (-x)**2 + y**2 <= 1.0:
            points_in_circle += 1

    pi_estimate = 4 * points_in_circle / num_samples  # 分母仍是总样本数
    return pi_estimate

print(f"对偶变量法估计值 (N=100000): {estimate_pi_antithetic(100000)}")

为什么有效？ 当 (x, y) 落在圆外时，(-x, y) 也倾向于落在圆外；反之亦然。这种负相关性使得两个样本的平均值方差小于两个独立样本的平均值方差。用同样多的总点数，对偶变量法的方差通常小于基础方法。

5. 引入第二种方差缩减技术：重要性采样

重要性采样的思想是：不均匀地撒点。我们主动让随机点更多地落在我们关心的区域（即单位圆）附近，从而“浪费”更少的点在无效区域（正方形但不在圆内的部分）。

具体到估计π，我们可以改变采样分布，使点在靠近原点的地方更密集。一个简单的做法是使用极坐标采样，但调整径向 r 的分布。

改进采样过程，让 r 的分布使得点在单位圆内更密集：

import numpy as np

def estimate_pi_importance(num_samples):
    """使用重要性采样估计π。"""
    points_in_circle = 0
    total_weight = 0.0  # 用于计算加权平均

    for _ in range(num_samples):
        # 采用新的采样分布：r的分布改为 f(r) = 2r (0<=r<=1)，这样在圆内更密集
        # 这相当于在面积上均匀采样（dA = 2πr dr -> f(r) ∝ r）
        r = np.sqrt(np.random.uniform(0, 1))  # 注意这里用平方根，与之前一致
        theta = np.random.uniform(0, 2 * np.pi)

        x = r * np.cos(theta)
        y = r * np.sin(theta)

        # 所有点都生成在单位圆内，因此总是满足条件
        # 但我们需要计算似然比（重要性权重）来纠正采样偏差
        # 目标分布（均匀在正方形）的概率密度 g(x,y) = 1/4
        # 采样分布（我们的新分布）的概率密度 h(r,θ) = (1/(2π)) * (2r) * |Jacobian|^{-1}
        # 计算复杂，这里演示一个更简单的思路：直接统计并加权

        # 一个更简单的实现：我们仍判断点是否在圆内（总是为真），但用原始均匀分布来生成对应的“虚拟点”以计算权重
        # 为了简化，我们采用一个等效方法：直接生成原始均匀分布的点，但用重要性权重调整贡献
        # 这里我们用另一种直接做法：在极坐标下，径向按f(r)=2r采样，保证圆内均匀。我们只需统计即可。
        points_in_circle += 1  # 因为所有生成的点都在圆内

    # 由于采样分布已经改变，估计量不再是简单的比例。
    # 正确的估计量是： π = 4 * E[ 1_{点在圆内} * (g(X)/h(X)) ]
    # 在我们的例子中，g(X)/h(X) 是一个常数因子（与点的位置无关），可以事先计算。
    # 对于在正方形内均匀采样（基础方法）vs 在圆内均匀采样（本例），
    # 似然比 g/h 在圆内是一个常数： (1/4) / (1/π) = π/4。
    # 因此，估计量变为： π ≈ 4 * (π/4) * (1) = π， 这是一个无意义的恒等式。
    # 这说明直接将点限制在圆内并不能估计π，我们需要的是对比。

    # 真正有效的策略是：使用一个在正方形内非均匀但与被积函数相关性高的分布。
    # 一个更实际的做法是：使用 f(x,y) ∝ 1/(1 + x^2 + y^2) 或其他在圆附近增大的分布。
    # 由于实现较复杂，此处仅说明原理。

    # 一个更简洁的演示：使用坐标变换。令 u = x, v = y * sqrt(1-x^2)，但这不直观。

    # 回到最实用且简单的方差缩减：对偶变量法已被证明有效。
    # 重要性采样需要仔细设计采样分布和计算权重，在此简单问题中收益不明显。
    # 但它在积分维度高、被积函数集中时威力巨大。

    # 为了完整，我们展示一个基于重要性采样的正确但简单的估计量
    # 重新设计：我们从一个与圆相关的分布中采样，例如以原点为中心的二元正态分布
    from scipy.stats import multivariate_normal
    mean = [0, 0]
    cov = [[0.5, 0], [0, 0.5]]  # 协方差矩阵，使样本集中在原点附近
    samples = multivariate_normal.rvs(mean=mean, cov=cov, size=num_samples)

    # 计算每个样本在目标分布（正方形均匀）下的概率密度 g(x,y) = 1/4
    g = 1.0 / 4.0
    # 计算每个样本在采样分布（二元正态）下的概率密度 h(x,y)
    h = multivariate_normal.pdf(samples, mean=mean, cov=cov)

    # 判断点是否在圆内
    r_sq = samples[:, 0]**2 + samples[:, 1]**2
    in_circle = (r_sq <= 1.0).astype(float)

    # 计算加权平均
    weights = g / h
    weighted_in_circle = np.mean(in_circle * weights)

    pi_estimate = 4 * weighted_in_circle
    return pi_estimate

print(f"重要性采样估计值 (N=100000): {estimate_pi_importance(100000)}")

关键点：重要性采样的估计量是 $\pi = 4 \times \mathbb{E}\left[ \mathbf{1}_{\{\text{点在圆内}\}} \times \frac{g(X)}{h(X)} \right]$，其中 $g$ 是目标分布（均匀）的概率密度，$h$ 是采样分布的概率密度。选择一个与被积函数（这里就是指示函数）形状相似的 $h$，可以显著降低方差。当 $h$ 集中在函数值大的区域（即圆内）时，权重 $g/h$ 在圆内较小且稳定，在圆外较大但乘以的指示函数为0，从而减少了总体波动。

6. 总结技术路径

首先，掌握基于几何概率的基础蒙特卡洛方法。然后，认识到其 $1/\sqrt{N}$ 的慢收敛速度是瓶颈。接着，学习并应用对偶变量法，这是一种简单有效的方差缩减技术，通过引入负相关样本来降低方差。最后，理解重要性采样的更强大思想，通过改变采样分布来聚焦计算资源，尽管其实现更复杂，但在处理高维或复杂积分问题时是核心工具。实践中，从基础方法开始，逐步集成这些方差缩减技术，可以在相同计算预算下获得更精确的估计。