[R] Different values for double integral

[Prof Brian Ripley]

More generally, if you want to do a two-dimensional integral, you will 
do better to us a 2D integration algorithm, such as those in package 
'adapt'.

Also, these routines are somewhat sensitive to scaling, so if the 
correct answer is around 5e-9, you ought to rescale.   You seem to be 
in the far right tail of your gamma distribution (but as you are not 
integrating over z, pgamma is not appropriate).

More specifically, cumsum(z)[-length(z)] is constant ( c(80,100, 140)
) are can be done once.  But without knowing the intention of f1, it 
is impossible to show you better code.

On Sat, 24 Jan 2009, Duncan Murdoch wrote:

> 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
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595