Как найти процентную ошибку в алгоритме Монте-Карло? - PullRequest
Новая
       24

Как найти процентную ошибку в алгоритме Монте-Карло?

0 голосов
/ 28 апреля 2019

Я написал программу Монте-Карло для интеграции функции f (x).

Теперь меня попросили вычислить процентную ошибку.

Сделав быстрый поиск литературы, я обнаружил, что это можно сделать с помощью уравнения% error = (sqrt (var [f (x)] / n)) * 100, где n - число случайных чисел.баллы, которые я использовал, чтобы получить свой ответ.

Однако, когда я запускаю свой код интеграции, моя процентная ошибка больше, чем указано в этой формуле.

У меня есть правильная формула?

Любая помощь будет принята с благодарностью.Спасибо х

1 Ответ

1 голос
/ 01 мая 2019

Вот краткий пример - оценка интеграла линейной функции на интервале [0 ... 1] с использованием метода Монте-Карло. Чтобы оценить ошибку, вы должны собрать второй импульс (значения в квадрате), а затем вычислить дисперсию, стандартное отклонение и (при условии CLT), ошибку моделирования в исходных единицах, а также в%

Код, Python 3.7, Анаконда, 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))

Имейте в виду, что мы вычислили одну сигма-ошибку и (даже не говоря о том, что она сама является случайным значением), истинный результат находится в пределах напечатанного среднего + -ошибки только для 68% прогонов. Вы можете вывести среднюю ошибку + -2 *, и это будет означать, что истинный результат находится внутри этой области для 95% случаев, средний + + 3 * истинный результат находится внутри этой области для 99,7% прогонов и т. Д. И т. Д.

UPDATE

Для оценки дисперсии выборки существует известная проблема, называемая Смещение в оценщике . По сути, мы недооцениваем небольшую дисперсию выборки, должна применяться правильная коррекция (коррекция Бесселя) 

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 в миллионах, коррекция не имеет практического применения.

...