为何蒙特卡洛法求π的Java程序最终结果精度低于过程值？

问题背景 #

使用蒙特卡洛法近似计算圆周率π时，发现π的估算精度并未随采样点增加逐步稳定提升：循环结束时的最终结果精度始终低于运行过程中出现的最优结果。最初使用double类型，怀疑是舍入误差导致，改用BigDecimal重写后问题依旧。即使将采样点总数MAX_PTS增大至5亿，最终误差仍收敛在约0.002%左右。

示例代码 #

import java.util.*;
import java.math.*;
import java.text.*;


public class PiMonte
{
  private static final int MAX_PTS = 500000;

  private static final BigDecimal BD100 = BigDecimal.valueOf(100);
  private static final BigDecimal BD10 = BigDecimal.valueOf(10);
  private static final BigDecimal BD4 = BigDecimal.valueOf(4);

  private static final BigDecimal PI = BigDecimal.valueOf(Math.PI);
  

  public static void main(String[] args)
  {
    Random rand = new Random();

    BigDecimal circlePts = BigDecimal.ZERO;
    BigDecimal totalPts = BigDecimal.ZERO;
    BigDecimal piEstimate = BigDecimal.ZERO;

    BigDecimal err = BigDecimal.ZERO;
    BigDecimal eps = BigDecimal.valueOf(100);  // percentage

    MathContext mc = new MathContext(20);

    System.out.println("Math.PI: \t" + Math.PI);
    System.out.println("PI Estimation");

    for (int i = 0; i < MAX_PTS; i++) {
      double x = rand.nextDouble()*2-1;  // in the range [-1,1]
      double y = rand.nextDouble()*2-1;

      if ((x*x + y*y) <= 1)  // inside circle
        circlePts = circlePts.add(BigDecimal.ONE);
      totalPts = totalPts.add(BigDecimal.ONE);

      piEstimate = circlePts.multiply(BD4).divide(totalPts, mc);
      err = piEstimate.subtract(PI).multiply(BD100).divide(PI, mc).abs();

      if (err.compareTo(eps) < 0) {
        System.out.println(" Step " + i + ": \t" + piEstimate + 
                               " < " + eps + "%; " + format(err) + "%");
        eps = eps.divide(BD10);
      }
    }

    System.out.println("\nEnd " + MAX_PTS + ": \t" + piEstimate + 
                                           " and " + format(err) + "%");

  } /* main */


  public static String format(BigDecimal num) 
  {
    DecimalFormat df = new java.text.DecimalFormat("0.##########");
    return df.format(num);
  }

}

程序运行示例（MAX_PTS=500000） #

> java PiMonte
Math.PI:        3.141592653589793
PI Estimation
 Step 2:        1.3333333333333333333 < 100%; 57.5586818422%
 Step 20:       2.8571428571428571429 < 10%; 9.0543182332%
 Step 107:      3.1111111111111111111 < 1%; 0.9702576317%
 Step 111:      3.1428571428571428571 < 0.1%; 0.0402499435%
 Step 190:      3.1413612565445026178 < 0.01%; 0.0073655967%
 Step 1327:     3.1415662650602409639 < 0.001%; 0.000839973%
 Step 24174:    3.1415925542916235781 < 0.0001%; 0.0000031608%
 Step 87878:    3.1415924168458903720 < 0.00001%; 0.0000075358%
 Step 186181:   3.1415926351634422232 < 0.000001%; 0.0000005865%

End 500000:     3.143968 and 0.0756096245%

现象原因解析 #

这个问题的核心是蒙特卡洛法的统计特性，和数值精度（double/BigDecimal）无关：

• 随机采样的波动性：蒙特卡洛法依赖随机点的均匀分布，估算值是围绕真实π值波动的随机变量。过程中出现的最优结果只是波动中的偶然低点，并非稳定收敛的状态。随着采样点增加，估算值的波动幅度会逐渐减小，但始终存在波动，最终结果恰好落在误差较大的区间是完全正常的。
• 误差收敛速率的限制：蒙特卡洛法的误差收敛速率是O(1/√N)（N为采样点数），也就是说，要将误差降低一个数量级，需要将采样点数增加100倍。5亿采样点对应的理论误差量级约为1/√(5e8) ≈ 0.00045%，但实际运行中由于随机波动，误差会在这个值附近上下浮动，你观察到的0.002%属于正常波动范围。
• 最优结果的偶然性：程序中记录的过程最优值，是在大量随机尝试中恰好出现的误差极小的情况，属于小概率事件。当采样点继续增加时，后续的随机点会拉低这个“偶然最优”的精度，让估算值回归到真实值附近的正常波动区间。

简单来说，过程最优值是随机波动带来的“运气值”，而最终结果才是当前采样点数下的典型收敛水平。多次运行程序会发现，最终结果的误差在理论范围内波动，过程中也会随机出现不同的最优精度点。

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