R help - Jan 2009 - Different values for double integral

Dear R useRs,

i have the function f1(x, y, z) which i want to integrate for x and y. 
On the one hand i do this by first integrating for x and then for y, on 
the other hand i do this the other way round and i wondering why i 
doesn't get the same result each way?

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)$value}
"g2" <- function(x, z) {integrate(function(y) {f1(x=x, y=y, z=z)},
0.1,
0.5)$value}

integrate(Vectorize(function(y) {g1(y=y, z=z)}), 0.1, 0.5)$value
[1] 5.909943e-09
integrate(Vectorize(function(x) {g2(x=x, z=z)}), 0.1, 0.5)$value
[1] 5.978334e-09

If you have any advice what is wrong here, i would be very thankful.


best regards
Andreas

On 24/01/2009 5:23 AM, Andreas Wittmann wrote:
> Dear R useRs,
> 
> i have the function f1(x, y, z) which i want to integrate for x and y. 
> On the one hand i do this by first integrating for x and then for y, on 
> the other hand i do this the other way round and i wondering why i 
> doesn't get the same result each way?
> 
> 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)$value}
> "g2" <- function(x, z) {integrate(function(y) {f1(x=x, y=y,
z=z)}, 0.1,
> 0.5)$value}
> 
> integrate(Vectorize(function(y) {g1(y=y, z=z)}), 0.1, 0.5)$value
> [1] 5.909943e-09
> integrate(Vectorize(function(x) {g2(x=x, z=z)}), 0.1, 0.5)$value
> [1] 5.978334e-09
> 
> If you have any advice what is wrong here, i would be very thankful.
It looks as though your f1 returns a vector result whose length will be 
the max of length(z)-1, length(x), and length(y):  that's not good, when 
you don't have control over the lengths of x and y.  I'd guess
that's
your problem.  I don't know what your intention was, but if the lengths 
of x and y were 1, I think it should return a length 1 result.

More geneally, integrate() does approximate integration, and it may be 
that you won't get identical results from switching the order even if 
you fix the problems above.

Finally, if this is the real problem you're working on, you can use the 
pgamma function to do your inner integral:  it will be faster and more 
accurate than integrate.

Duncan Murdoch

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.