Использование функции 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 (Вики-страница по дистрибутиву Gamma). Более того, если я поставлю 10 ^ 4 в качестве верхней границы (нижняя граница по умолчанию равна 0) для variance(), она вернет следующее:

variance(10^4)
## [1] 18

И интеграл от 10^4 до 10^5 будет:

variance(10^5) - variance(10^4)
## [1] -18

Кто-нибудь знает, почему variance(10^5) в этом случае дает бессмысленные результаты? Так же буду признателен за любые дополнительные комментарии по стилю поста.