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

发表于： 2025-03-18

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

title: “蒙特卡洛方法计算圆周率：积分公式选择的方差分析” description: keywords: “蒙特卡洛方法计算圆周率：积分公式选择的方差分析” date: 2025-03-18T00:22:27+08:00 lastmod: 2025-03-18T00:22:27+08:00 categories:

科研笔记

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低一个数量级。
数值稳定性：
无需处理函数发散问题，所有采样点的贡献均在可控范围内。

科研笔记

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$）是更优选择，原因如下：

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

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