1

This is my first post here and I hope I'll follow all the rules of the community.

I'm trying to calculate variance of gamma distribution with shape parameter 2 and scale parameter 3 in R using function antiD from mosaic package. The R code I use is the following

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

The problem is that the number I get doesn't make sense as the variance for Gamma distribution with those parameters should be equal to 2*3^2 = 18 (Wiki page on Gamma distribution). Moreover if I put 10^4 as an upper bound (the default lower bound is 0) for variance() it will return the following:

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

And the integral from 10^4 to 10^5 will be:

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

Does anyone know why variance(10^5) produce nonsensical results in this case? I also will be grateful for any additional comments on the style of the post.

Nicholas Horton
  • 19
  • 1
Georgiy Syunyaev
  • 11
  • 1
  • Don't know the answer, although I think someone will find it in the code for `mosaic::numerical.integration`. Your question is clear and complete. – IRTFM Sep 21 '14 at 09:37
  • `antiD()` is a complex wrapper (through `numerical_integration()`) for the `integrate()` function in base R. You can see the same results much more simply with `integrate(function(x) dgamma(x,shape=2),lower=0,upper=10^5)` (or `upper=10^4`). The problem is that adaptive numerical integration doesn't work well if the range is so large that the initial mesh doesn't "see" any non-trivial values. Interestingly, `upper=Inf` actually works better. – Ben Bolker Sep 21 '14 at 15:59
  • closely related: http://stackoverflow.com/questions/20665689/integrate-a-very-peaked-function-in-r – Ben Bolker Sep 21 '14 at 16:02
  • @BenBolker Thank you very much, I think using (upper=Inf) solves the issue for me. – Georgiy Syunyaev Sep 21 '14 at 16:39

0 Answers0