蒙特卡洛方法计算圆周率：积分公式选择的方差分析

发表于 更新于

在使用蒙特卡洛方法计算圆周率时，选择积分公式的关键在于方差的大小和实现的效率。本文通过数学推导和实验验证，对比两种常见积分公式的优劣。

1. 积分公式的数学等价性

两个积分的数学结果均为π，但被积函数特性不同：

  • 公式1：
    \[ 2 \int_0^1 \frac{1}{\sqrt{1 - x^2}} \, dx = \pi \]
    被积函数 \(f_1(x) = \frac{1}{\sqrt{1 - x^2}}\) 在 \(x \to 1\) 时发散，导致数值计算的不稳定性。

  • 公式2：
    \[ 2 \int_{-1}^1 \sqrt{1 - x^2} \, dx = \pi \]
    被积函数 \(f_2(x) = \sqrt{1 - x^2}\) 在区间 \([-1, 1]\) 上有界且连续，方差更小。

2. 蒙特卡洛方法的方差分析

2.1 公式1的方差问题

  • 发散特性：
    \(f_1(x)\) 在 \(x \to 1\) 处趋向无穷大，导致蒙特卡洛采样时方差显著增大。例如，当 \(x = 0.99\) 时，\(f_1(x) \approx 7\)，远高于平均值 \(\pi/2 \approx 1.57\)，极端值会严重干扰估计结果。

  • 实现挑战：
    需要均匀采样 \(x \in [0, 1]\)，但靠近 \(x = 1\) 的点对结果贡献极大，易引入噪声。

2.2 公式2的方差优势

  • 有界性：
    \(f_2(x)\) 的值域为 \([0, 1]\)，方差远小于公式1。例如，\(N=10^6\) 次采样时，公式2的估计误差通常比公式1低一个数量级。

  • 数值稳定性：
    无需处理函数发散问题，所有采样点的贡献均在可控范围内。

3. 实验验证：方差的量化对比

假设使用 \(N\) 次采样：

  • 公式1：
    由于 \(f_1(x)\) 的发散性，方差 \(D_1 \propto \int_0^1 \left( \frac{1}{\sqrt{1-x^2}} \right)^2 dx = \infty\)（发散）。

  • 公式2：
    方差 \(D_2 = \int_{-1}^1 \left( \sqrt{1-x^2} \right)^2 dx - \left( \frac{\pi}{2} \right)^2 = \frac{\pi}{2} - \frac{\pi^2}{4} \approx 0.43\)，远小于公式1。

4. 代码实现对比

4.1 公式2的C++实现（高效且低方差）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
import random
import math

def estimate_pi(n):
    count = 0
    for _ in range(n):
        x = random.uniform(-1, 1)       # 生成x ∈ [-1, 1]
        y = random.uniform(0, 1)        # 生成y ∈ [0, 1]
        if y <= math.sqrt(1 - x**2):    # 判断是否在半圆内
            count += 1
    # 计算公式：2 * (积分值) = 2 * (count/N * 矩形面积)
    return 2.0 * (count * 2.0 / n)      # 矩形面积为2（宽2×高1）

# 示例：使用100万个样本
print(f"Estimated π: {estimate_pi(1000000)}")

4.2 公式1的Python实现（高方差问题）

1
2
3
4
5
6
7
8
9
10
import numpy as np

def monte_carlo_formula1(n):
    x = np.random.uniform(0, 1, n)
    f = 1 / np.sqrt(1 - x**2)
    integral = np.mean(f)
    return 2 * integral

# 测试示例（方差极大）
print(monte_carlo_formula1(1000000))  # 结果可能波动明显

5. 结论

公式2（\(2\int_{-1}^{1} \sqrt{1 - x^2} \, dx\)）是更优选择，原因如下：

  1. 方差更小：被积函数 \(\sqrt{1 - x^2}\) 在区间 \([-1, 1]\) 上有界且连续，使得蒙特卡洛方法的收敛速度更快。
  2. 数值稳定性：无需处理函数发散问题（如公式1中 \(\frac{1}{\sqrt{1-x^2}}\) 在 \(x \to 1\) 处的发散），避免极端值对估计结果的干扰。
  3. 计算效率：在相同采样次数下，公式2的估计误差显著低于公式1，且计算过程更稳定。

选择积分公式时，应优先考虑被积函数的有界性和光滑性，以降低蒙特卡洛方法的方差，从而提升计算效率。