Интеграция Монте-Карло - Как найти ошибку - PullRequest
Новая
       67

Интеграция Монте-Карло - Как найти ошибку

2 голосов
/ 25 марта 2020

Мне нужно применить интеграцию Монте-Карло для оценки следующей функции в [0,1]: 

f(x) = exp(-ax) cos(bx)

, где a=0.3060734 и b=0.11221230. Но мне нужно использовать четыре различных варианта Монте-Карло: хит или промах, грубый, выборка важности и контроль дисперсии. Все они на странице 392 из этой книги: https://edisciplinas.usp.br/pluginfile.php/5168099/mod_resource/content/1/Julio%20Stern.pdf. Однако после оценки интеграла от f (x) в каждом варианте мне нужно вычислить относительную ошибку ( | g* - g | / g ) < 1%, где g - это реальное значение интеграла (неизвестно), а g* - это оценочное значение. Как рассчитать ошибку? Я думал об использовании дисперсии каждого варианта, но я не уверен, как это сделать и сделать его <1%. У меня уже есть код, вычисляющий каждую оценку и дисперсию. OBS: я должен использовать Python или R. </p>

1 Ответ

0 голосов
/ 25 марта 2020

Вот пример (с R) для симуляции Монте-Карло для интеграции 

a <- 0.3060734
b <- 0.11221230
f <- function(x) exp(-a*x)*cos(b*x)

MCsim <- function(n) {
  gs <- sapply(n,function(k) mean(runif(k) <= f(runif(k))))
  g <- integrate(f,0,1)$value
  rel_err <- abs(gs-g)/g
  data.frame(n,gs,g,rel_err)
}

n <- 10**(1:7)
res <- MCsim(n)

такой, что 

> res
      n        gs         g      rel_err
1 1e+01 0.8000000 0.8597815 6.953106e-02
2 1e+02 0.8700000 0.8597815 1.188497e-02
3 1e+03 0.8590000 0.8597815 9.089768e-04
4 1e+04 0.8570000 0.8597815 3.235149e-03
5 1e+05 0.8606900 0.8597815 1.056639e-03
6 1e+06 0.8596650 0.8597815 1.355245e-04
7 1e+07 0.8597802 0.8597815 1.536980e-06
...