[R] Different values for double integral

[Duncan Murdoch]

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