В последующем я намеренно не выставляю случайное семя.
Как я уже упоминал в моих комментариях, есть по крайней мере два вводных вопроса и ответа по интеграции Монте-Карло в переполнении стека:
Оба объяснили, как получить оценку Монте-Карло, но забыли о стандартной ошибке оценки. Просто получается, что оценка Монте-Карло имеет чрезвычайно медленную скорость сходимости по вашей функции.
Общеизвестно, что интеграция Монте-Карло имеет скорость сходимости O(1 / sqrt(N)), где N - это размер выборки, а O() - это big O нотация . Однако константа, стоящая за этим большим O , может быть очень большой для некоторых функций, поэтому реалистическая скорость сходимости может быть намного медленнее.
Ваши функции могут быть определены как минимум двумя способами:
## direct definition
f <- function (x) x^4 * sin(x) * exp(-x/5)
## using gamma distribution; see ?rgamma
g <- function (x) sin(x) * 5 ^ 5 * gamma(5) * dgamma(x, 5, 1/5)
curve(f, from = 0, to = 100)
curve(g, add = TRUE, col = 2)
В первом разделе вопросов и ответов объяснялось, как вычислить интеграцию Монте-Карло с использованием равномерно распределенных выборок. Ваша функция f или g почти равна нулю за пределами x = 200, поэтому интеграция на [0, +Inf) эффективно на [0, 200]. Следующая функция вернет вам интеграцию и ее стандартную ошибку:
MCI1 <- function (n) {
x <- runif(n, 0, 200)
y <- 200 * f(x)
c(mean.default(y), sqrt(var(y) / n))
}
Другой способ - выборка по важности, как объяснено во 2-м вопросе и ответе. Здесь гамма-распределение используется в качестве распределения предложений (как предложил Бен Болкер).
MCI2 <- function (n) {
x <- rgamma(n, 5, 0.2)
y <- sin(x) * 75000
c(mean.default(y), sqrt(var(y) / n))
}
Теперь давайте проверим скорость сходимости.
n <- seq(1000, by = 5000, length = 100)
tail(n)
#[1] 471000 476000 481000 486000 491000 496000
b1 <- sapply(n, MCI1)
b2 <- sapply(n, MCI2)
Для равномерной выборки имеем
par(mfrow = c(1, 2))
plot(b1[1, ], main = "estimate")
plot(b1[2, ], main = "standard error")
b1[, (ncol(b1) - 5):ncol(b1)]
# [,1] [,2] [,3] [,4] [,5] [,6]
#[1,] 115.1243 239.9631 55.57149 -325.8631 -140.3745 78.61126
#[2,] 181.0025 179.9988 178.99367 178.2152 177.2193 175.31446
Для гамма-выборки мы имеем
par(mfrow = c(1, 2))
plot(b2[1, ], main = "estimate")
plot(b2[2, ], main = "standard error")
b2[, (ncol(b2) - 5):ncol(b2)]
# [,1] [,2] [,3] [,4] [,5] [,6]
#[1,] -100.70344 -150.71536 24.40841 -49.58032 169.85385 122.81731
#[2,] 77.22445 76.85013 76.53198 76.03692 75.69819 75.25755
Какой бы это ни был метод, обратите внимание, насколько велика стандартная ошибка (по сравнению с самой оценкой) и насколько медленно она уменьшается.
Гораздо проще использовать числовое интегрирование (что неудивительно для интеграции одномерных функций):
integrate(f, 0, 200)
#11.99356 with absolute error < 0.0012
## trapezoidal rule
200 * mean.default(f(seq(0, 200, length = 10000)))
#[1] 11.99236
В трапециевидном правиле, даже если взяты только 1e + 4 равномерно распределенных точки отбора, интегрирование достаточно близко к истине.
Примечание
Интеграция с Монте-Карло была бы менее трудной, если бы мы делали интеграцию в более ограниченной области. Из рисунка f или g мы видим, что это колебательная функция. И на самом деле, он пересекает ось X с периодом pi. Давайте рассмотрим интеграцию на [lower, upper].
MCI3 <- function (n, lower, upper) {
x <- runif(n, lower, upper)
y <- (upper - lower) * f(x)
c(mean.default(y), sqrt(var(y) / n))
}
a1 <- sapply(n, MCI3, lower = 0, upper = pi)
a2 <- sapply(n, MCI3, lower = pi, upper = 2 * pi)
a3 <- sapply(n, MCI3, lower = 2 * pi, upper = 3 * pi)
a4 <- sapply(n, MCI3, lower = 3 * pi, upper = 4 * pi)
a1[, (ncol(a1) - 5):ncol(a1)]
# [,1] [,2] [,3] [,4] [,5] [,6]
#[1,] 17.04658711 16.97935808 17.01094302 17.02117843 16.96935285 16.99552898
#[2,] 0.02407643 0.02390894 0.02379678 0.02368683 0.02354298 0.02342799
a2[, (ncol(a2) - 5):ncol(a2)]
# [,1] [,2] [,3] [,4] [,5]
#[1,] -406.5646843 -404.9633321 -405.4300941 -405.4799659 -405.8337416
#[2,] 0.3476975 0.3463621 0.3442497 0.3425202 0.3409073
# [,6]
#[1,] -405.8628741
#[2,] 0.3390045
a3[, (ncol(a3) - 5):ncol(a3)]
# [,1] [,2] [,3] [,4] [,5] [,6]
#[1,] 1591.539911 1592.280780 1594.307951 1591.375340 1593.171500 1591.648529
#[2,] 1.197469 1.190251 1.183095 1.177079 1.172049 1.165667
a4[, (ncol(a4) - 5):ncol(a4)]
# [,1] [,2] [,3] [,4] [,5]
#[1,] -3235.561677 -3239.147235 -3241.532097 -3238.421556 -3238.667702
#[2,] 2.336684 2.321283 2.311647 2.300856 2.286624
# [,6]
#[1,] -3237.043068
#[2,] 2.279032