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@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.
