[R] Different values for double integral

[Christos Argyropoulos]

Each of the two integrals (g1, g2) seem to be divergent (or at least is considered to be so by R :) )
Try this:

z <- c(80, 20, 40, 30)

"f1" <- function(x, y, z) {dgamma(cumsum(z)[-length(z)], shape=x, rate=y)}

"g1" <- function(y, z) {integrate(function(x) {f1(x=x, y=y, z=z)}, 0.1,
0.5, rel.tol=1.5e-20)$value}
"g2" <- function(x, z) {integrate(function(y) {f1(x=x, y=y, z=z)}, 0.1,
0.5,rel.tol=1.5e-20)$value}



integrate(Vectorize(function(y) {g1(y=y, z=z)}), 0.1, 0.5)$value
Error in integrate(function(x) { : the integral is probably divergent

integrate(Vectorize(function(x) {g2(x=x, z=z)}), 0.1, 0.5)$value
Error in integrate(function(x) { : the integral is probably divergent

Are you sure:
a) you are not integrating over a singularity (in that case, you'd have to calculate the integral as a Cauchy principal value)
b) you are not over/underflowing during your "manual integration" of the Gamma density (would R's integrate function report this? - I don't know)
c) there is a  loss of precision during numerical integration inside the g1/g2 functions

To see that c) is indeed possible try your code with z set to

1) c(8,2,4,3)*0.5 :
> integrate(Vectorize(function(y) {g1(y=y, z=z)}), 0.1, 0.5)$value
[1] 0.002982671
> integrate(Vectorize(function(x) {g2(x=x, z=z)}), 0.1, 0.5)$value
[1] 0.002985362

2)   c(8,2,4,3) :
> integrate(Vectorize(function(y) {g1(y=y, z=z)}), 0.1, 0.5)$value
[1] 0.0006453548
> integrate(Vectorize(function(x) {g2(x=x, z=z)}), 0.1, 0.5)$value
[1] 0.0006466886

3)  c(8,2,4,3)*5:
> integrate(Vectorize(function(y) {g1(y=y, z=z)}), 0.1, 0.5)$value
[1] 1.198051e-06
> integrate(Vectorize(function(x) {g2(x=x, z=z)}), 0.1, 0.5)$value
[1] 1.253991e-06

4) c(8,2,4,3)*20:
> integrate(Vectorize(function(y) {g1(y=y, z=z)}), 0.1, 0.5)$value
[1] 5.832236e-13
> integrate(Vectorize(function(x) {g2(x=x, z=z)}), 0.1, 0.5)$value
[1] 3.230487e-13

to get an idea of how this would play out.

You could try to tweak integrate (e.g. increasing the maximum number of subdivisions, changing the absolute / relative precision) but
your integrand (? are you trying to marginalize over nuisance parameters in a Bayesian setting) is numerically ill-behaved that I do not think
you may be able to integrate it in a straightforward manner.
You may try want to try:
a) adapt (from library adapt) to carry out the integration over the shape and scale simultaneously
b) evaluate one or (even both!) integrations symbolically if possible (by hand, tables, or CAS e.g. Mathematica or Maxima) and recode accordingly

Click on the following link to  symbolically evaluate the indefinite numerical integral of the gamma density w.r.t to its scale parameter:
http://integrals.wolfram.com/index.jsp?expr=y^(a-1)*Exp[-y%2Fx]%2F(x^s*Gamma[a])&random=false

This will give you an idea of the numerical nastiness that the poor integrate function has to deal with .. :)




_________________________________________________________________
Show them the way! Add maps and directions to your party invites.