如何找到蒙特卡洛算法中的百分比误差？

Question

我编写了一个蒙特卡洛程序来积分函数 f(x)。

我现在被要求计算百分比误差。

通过快速的文献搜索，我发现这可以用公式 %error = (sqrt(var[f(x)]/n))*100 给出，其中 n 是我使用的随机点数得出我的答案。

但是，当我运行我的集成代码时，我的百分比误差大于这个公式给出的百分比。

我有正确的公式吗？

任何帮助将不胜感激。谢谢x

Answer 1

这是一个简单的示例 - 使用 Monte-Carlo 估计区间 [0...1] 上的线性函数的积分。要估计误差，您必须收集第二个动量（值的平方），然后以原始单位和 % 计算方差、标准差和（假设 CLT）模拟误差

代码、Python 3.7、Anaconda、Win10 64x

import numpy as np

def f(x): # linear function to integrate
    return x

np.random.seed(312345)

N = 100000

x  = np.random.random(N)
q  = f(x)  # first momentum
q2 = q*q   # second momentum

mean = np.sum(q) / float(N) # compute mean explicitly, not using np.mean
var  = np.sum(q2) / float(N) - mean * mean # variance as E[X^2] - E[X]^2
sd   = np.sqrt(var) # std.deviation

print(mean) # should be 1/2
print(var)  # should be 1/12
print(sd)   # should be 0.5/sqrt(3)
print("-----------------------------------------------------")

sigma = sd / np.sqrt(float(N)) # assuming CLT, error estimation in original units

print("result = {0} with error +- {1}".format(mean, sigma))

err_pct = sigma / mean * 100.0 # error estimate in percents

print("result = {0} with error +- {1}%".format(mean, err_pct))

请注意，我们计算了一个 sigma 误差，并且（即使不是说它本身是随机值）真实结果仅在 68% 的运行中处于打印的均值 +-误差范围内。您可以打印 mean+-2*error，这意味着 95% 的情况下真实结果在该区域内，对于 99.7% 的运行，mean+-3*error 真实结果在该区域内，依此类推。

更新

对于抽样方差估计，有一个称为Bias in the estimator 的已知问题。基本上，我们低估了一点采样方差，应该应用适当的校正（贝塞尔校正）

var  = np.sum(q2) / float(N) - mean * mean # variance as E[X^2] - E[X]^2
var *= float(N)/float(N-1)

在许多情况下（和许多示例），它被省略了，因为 N 非常大，这使得校正几乎不可见 - fe，如果您有 1% 的统计误差但 N 以百万为单位，则校正没有实际用途。