Использование функции antiD для дисперсии гамма-распределения
Это мой первый пост здесь, и я надеюсь, что буду следовать всем правилам сообщества.
Я пытаюсь вычислить дисперсию гамма-распределения с параметром формы 2 и параметром масштаба 3 в R, используя функцию antiD из пакета мозаики. Код R, который я использую, следующий
stopifnot(require(mosaic))
f <- function(y) {
dgamma(y, shape = 2, scale = 3)
}
mean_integral <- antiD( z*f(z) ~ z )
mn <- mean_integral(10^4)
g <- function(y) {
(y - mn)^2
}
variance <- antiD(f(x)*g(x) ~ x)
variance(10^5)
## [1] 7.115334e-09
Проблема в том, что число, которое я получаю, не имеет смысла, так как дисперсия для гамма-распределения с этими параметрами должна быть равна 2*3^2 = 18 ( вики-страница о гамма-распределении). Более того, если я положу 10^4 в качестве верхней границы (по умолчанию нижняя граница равна 0) для variance(), он вернет следующее:
variance(10^4)
## [1] 18
И интеграл от 10^4 до 10^5 будет:
variance(10^5) - variance(10^4)
## [1] -18
Кто-нибудь знает почему variance(10^5) привести к бессмысленным результатам в этом случае? Буду также признателен за любые дополнительные комментарии по стилю поста.