为何While循环估算π比固定次数向量化方法用更少样本更精准？

问题背景 #

我尝试用随机数估算π值：最初采用固定生成1000000个随机数的向量化方法，统计落入单位圆的点数来估算π；但使用While循环方法时，仅用更少的抽样次数就得到了精度更高的π估算结果。

初始固定100万次抽样的向量化代码 #

import numpy as np
import matplotlib.pyplot as plt

def pi_estimation1(n):
    # 创建随机数生成器
    rng = np.random.default_rng()

    # 同时生成x和y坐标
    coordinates = rng.random((n, 2))

    # 用向量化操作筛选圆内/外的点
    x_in_circle = coordinates[coordinates[:, 0]**2 + coordinates[:, 1]**2 <= 1]
    x_out_circle = coordinates[coordinates[:, 0]**2 + coordinates[:, 1]**2 > 1]

    # 根据圆内点比例计算π值
    pi = 4 * len(x_in_circle) / n
    return pi, x_in_circle, x_out_circle,n

# 使用示例
estimated_pi, x_in_circle, x_out_circle,n = pi_estimation1(1000000)
# 绘制圆内/外的点
plt.figure(figsize=(5, 5))
plt.scatter(x_in_circle[:, 0], x_in_circle[:, 1], color='blue')
plt.scatter(x_out_circle[:, 0], x_out_circle[:, 1], color='red')
plt.show()

print("估算的π值为 ", estimated_pi)
print("精确π值为 ", np.pi)
print("相对误差为 ", abs(estimated_pi - np.pi) / np.pi)
print("抽样点数 ",n)

运行结果

估算的π值为  3.140236
精确π值为  3.141592653589793
相对误差为  0.0004318362497579398
抽样点数  1000000

While循环实现的π估算代码 #

import numpy as np

def pi_estimation_while_loop(deviation):
    
    # 初始化变量
    pi = 0
    n = 0
    dev=1
    in_circle = 0
    
    rng=np.random.default_rng()
    
    # 循环直到偏差小于设定阈值
    while deviation < dev :
        coordinate = rng.random((1, 2))
        if coordinate[0,0]**2 + coordinate[0,1]**2 <= 1:
            in_circle += 1
        
    
        n += 1
        pi = 4 * in_circle / n
        dev=(abs(pi-np.pi)/np.pi)
        
    return pi, n, dev

pi, n, dev = pi_estimation_while_loop(1e-9)
print("估算的π值为 ", pi)
print("精确π值为 ", np.pi)
print("相对误差为 ", abs(pi - np.pi) / np.pi)
print("迭代次数为 ", n)

运行结果

估算的π值为  3.141592654644592
精确π值为  3.141592653589793
相对误差为  3.3575292254192973e-10
迭代次数为  499581

精度差异的核心原因 #

1. 终止逻辑本质不同
   向量化方法是固定抽样次数：不管当前估算精度是否达标，生成100万样本后就停止，这次运行刚好得到了误差约4e-4的结果。而While循环是精度驱动的动态终止：持续抽样直到相对误差低于设定的1e-9阈值，相当于“按需抽样”，刚好在满足精度要求时停止，所以抽样次数更少但精度更高。

2. 随机样本的偶然性
   蒙特卡洛方法本身依赖随机样本，存在天然的波动。向量化那次的100万样本刚好偏差稍大，而While循环的样本序列刚好收敛到更高精度——但这是次要因素，核心还是终止条件的差异。如果多次运行向量化方法，你会发现误差在一定范围内波动，而While循环每次都会停在满足精度的次数上。

用最少抽样次数优化π估算的方案 #

1. 优化动态终止条件 #

不要直接用“当前估算值与真实π的误差”（因为你其实不知道真实π的精确值），改用置信区间判断：
蒙特卡洛估算的π的标准差为 $\sigma = \sqrt{\frac{\pi(4-\pi)}{n}}$，95%置信区间为 $\hat{\pi} \pm 1.96\sigma$。当置信区间的半宽小于你需要的精度时停止抽样，这样更科学，避免依赖对真实π的提前认知。

2. 采用低方差蒙特卡洛变种 #

这些方法能在相同抽样次数下降低误差，或者用更少次数达到相同精度：

• 分层抽样：把[0,1]×[0,1]区域分成若干小格子，每个格子内抽样相同次数，确保样本分布更均匀，减少随机波动。
• 重要性抽样：在单位圆附近增加抽样密度，因为这些点对π估算的贡献更大，能有效降低方差。
• 对称性利用：只在第一象限抽样（x∈[0,1], y∈[0,1]），计算圆内点比例后乘以4，这样抽样次数减半，同时方差也会降低。

3. 向量化+动态终止结合 #

单样本While循环速度较慢，可改为批量向量化抽样：每次生成一批样本（比如10000个），更新当前估算值和误差，未达到精度要求就继续生成下一批，兼顾向量化的计算效率和动态终止的抽样经济性。

示例代码（批量向量化+置信区间终止）：

import numpy as np

def pi_estimation_optimal(target_precision):
    rng = np.random.default_rng()
    n = 0
    in_circle = 0
    batch_size = 10000
    
    while True:
        # 批量生成样本
        coords = rng.random((batch_size, 2))
        in_batch = np.sum(coords[:,0]**2 + coords[:,1]**2 <= 1)
        in_circle += in_batch
        n += batch_size
        
        # 计算当前估算值和置信区间半宽
        pi_hat = 4 * in_circle / n
        sigma = np.sqrt((pi_hat * (4 - pi_hat)) / n)
        ci_half_width = 1.96 * sigma
        
        # 当置信区间半宽小于目标精度时停止
        if ci_half_width < target_precision:
            break
    
    return pi_hat, n, ci_half_width

pi_hat, n, ci_width = pi_estimation_optimal(1e-5)
print(f"估算π值: {pi_hat}")
print(f"抽样次数: {n}")
print(f"95%置信区间半宽: {ci_width}")

内容的提问来源于stack exchange，提问作者PortorogasDS