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

Question

我编写了一个蒙特卡洛程序来集成函数f（x）。

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

做了一个快速的文献检索，我发现这可以用方程式给出％error =（sqrt（var [f（x）] / n））* 100，其中n是我用来导出我的随机点的数量回答。

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

我有正确的配方吗？

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

Answer 1

这是一个快速示例 - 使用蒙特卡罗估计区间[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错误并且（即使不是在谈论它是随机值本身）真实结果在打印平均值+ -error内只有68％的运行。您可以打印平均值+ -2 *错误，这意味着真实结果在95％的情况下在该区域内，平均值+ -3 *错误真实结果在该区域内99.7％的运行，依此类推。

UPDATE

对于抽样方差估计，存在称为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非常大，这使得校正几乎不可见 - 例如，如果你有1％的统计误差但N是百万，那么校正没有实际用途。