蒙特卡洛法估算圆周率π时重复次数上升σ异常升高问题求助

问题背景 #

我正在使用Python完成习题3.1，题目如下：

实现代码 #

我编写的可正常运行的代码如下：

#NUMERICAL ESTIMATE OF PI
import numpy as np               # 数值计算库
import matplotlib.pyplot as plt  # 绘图库
from scipy.stats import norm     # 高斯拟合依赖库
  
#*******************************************************************************
M = 10**2                               # 估算π的重复次数
N = 10**4                               # 单次估算生成的随机点数量
mean_pi=[]                              # 存储每次估算π值的空列表
for i in range(M):
    x=np.random.uniform(-1,1,N)         # 生成[-1,1)区间均匀分布的x坐标数组
    y=np.random.uniform(-1,1,N)         # 生成[-1,1)区间均匀分布的y坐标数组
    x_sel=x[(x**2+y**2)<=1]             # 筛选落在单位圆内的x坐标
    y_sel=y[(x**2+y**2)<=1]             # 筛选落在单位圆内的y坐标
    mean_pi+=[4*len(x_sel)/len(x)]      # 计算本次估算的π值存入列表
#*******************************************************************************
plt.figure(figsize=(8,3))
# 绘制直方图，将分箱范围存入bins变量
_,bins,_=plt.hist(mean_pi,bins=int(np.sqrt(N)),density=True, color="skyblue")
mu, sigma = norm.fit(mean_pi)                                     # 拟合得到数据的均值和标准差
k = sigma*np.sqrt(N)                                              # 计算k参数
best_fit_line = norm.pdf(bins, mu, sigma)                         # 生成拟合曲线数据
print("\nTime of repetitions:", M, ". The mean of the distribution is: ", mu, ". The standard deviation is:", sigma, ". The k parameters is:", k ,". \n")
#*******************************************************************************
plt.plot(bins, best_fit_line, color="red")
plt.grid()
plt.xlabel('Bins',fontweight='bold')
plt.ylabel('Pi',fontweight='bold')
plt.title("Histogram for Pi vs. bins")
plt.show()
print("\n")
#*******************************************************************************
M = 10**3                               # 调整重复次数为1000
N = 10**4                               # 单次随机点数量保持不变
mean_pi=[]
for i in range(M):
    x=np.random.uniform(-1,1,N)
    y=np.random.uniform(-1,1,N)
    x_sel=x[(x**2+y**2)<=1]
    y_sel=y[(x**2+y**2)<=1]
    mean_pi+=[4*len(x_sel)/len(x)]
#*******************************************************************************
plt.figure(figsize=(8,3))
_,bins,_=plt.hist(mean_pi,bins=int(np.sqrt(N)),density=True, color="skyblue")
mu, sigma = norm.fit(mean_pi)
k = sigma*np.sqrt(N)
best_fit_line = norm.pdf(bins, mu, sigma)
print("Time of repetitions:", M, ". The mean of the distribution is: ", mu, ". The standard deviation is:", sigma, ". The k parameters is:", k ,". \n")
#*******************************************************************************
plt.plot(bins, best_fit_line, color="red")
plt.grid()
plt.xlabel('Bins',fontweight='bold')
plt.ylabel('Pi',fontweight='bold')
plt.title("Histogram for Pi vs. bins")
plt.show()
print("\n")
#*******************************************************************************
M = 5*10**3                             # 调整重复次数为5000
N = 10**4                               # 单次随机点数量保持不变
mean_pi=[]
for i in range(M):
    x=np.random.uniform(-1,1,N)
    y=np.random.uniform(-1,1,N)
    x_sel=x[(x**2+y**2)<=1]
    y_sel=y[(x**2+y**2)<=1]
    mean_pi+=[4*len(x_sel)/len(x)]
#*******************************************************************************
plt.figure(figsize=(8,3))
_,bins,_=plt.hist(mean_pi,bins=int(np.sqrt(N)),density=True, color="skyblue")
mu, sigma = norm.fit(mean_pi)
k = sigma*np.sqrt(N)
best_fit_line = norm.pdf(bins, mu, sigma)
print("Time of repetitions:", M, ". The mean of the distribution is: ", mu, ". The standard deviation is:", sigma, ". The k parameters is:", k ,". \n")
#*******************************************************************************
plt.plot(bins, best_fit_line, color="red")
plt.grid()
plt.xlabel('Bins',fontweight='bold')
plt.ylabel('Pi',fontweight='bold')
plt.title("Histogram for Pi vs. bins")
plt.show()
#*******************************************************************************
print("\n How many couples N you need to estimate pi at better than 0.0001? The number of couples N is:", (k**2)*10**8 ,".")
#*******************************************************************************

运行结果 #

运行输出结果如下图：

问题描述 #

运行结果和预期不符：随着估算π的重复次数M从100提升到1000、5000，拟合得到的π值分布标准差σ反而上升，按照预期重复次数越多σ应该越小。尝试增大单次运行生成的随机点数量N，结果也没有改善，需要排查错误原因。

排查结论 #

• 核心问题是统计概念混淆：你计算的σ是单次π估算值的分布标准差，这个值的理论值仅和单次估算的随机点数量N有关，和重复次数M没有关系。对于蒙特卡洛估算π的场景，单次估算的标准差理论值约为$\frac{1.64}{\sqrt{N}}$，当N固定为1e4时，这个值约为0.0164。你观测到的σ随M增大轻微上升只是小样本下的统计波动，当M足够大时，你计算得到的σ会收敛到这个固定的理论值，而不会下降。
• 你预期会随M增大下降的是M次估算结果均值的标准差（标准误），这个值的计算公式是$\frac{\sigma}{\sqrt{M}}$，确实会随M增大而减小，但它和你当前计算的单次分布标准差σ不是同一个统计量。
• 代码本身的计算逻辑没有问题，你可以尝试计算每组M次实验得到的mean_pi的均值，会发现这个均值的波动随M增大明显变小，符合你最初的预期。

内容的提问来源于stack exchange，提问作者J.Snowden