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
Different values for double integral
4 messages · Duncan Murdoch, Brian Ripley, Andreas Wittmann
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
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
Thank you all, for the very helpful advice.
i want to estimate the parameters omega and beta of the gamma-nhpp-model
with numerical integration.
so one step in order to do this is to compute the normalizing constant,
but as you see below i get different values
## some reliability data, theses are times between events like a
failures during software test.
`s` <-
c(8100, 4800, 900, 450, 450, 6000, 2400, 2100, 2100,
1260, 1200, 300, 9000, 600, 2400, 600, 15060, 120, 360,
1200, 300, 1200, 1200, 3300, 19800, 3000, 600, 9600,
8400, 8100, 15600, 1800, 5400, 3900, 2400, 1200, 79500,
9000, 48900)
## likelihood function of the NHPP-gamma reliability model
`likelihood` <- function(s, omega, beta, alpha=1)
{
me<-length(s)-1
te<-cumsum(s)[length(s)] ## endpoint of observation
omega^me*prod(dgamma(cumsum(s)[-length(s)], shape=alpha, rate=beta)) *
exp(-omega*pgamma(te, shape=alpha, rate=beta))
}
## normalizing constant first integrating omega then beta
`normConst1` <- function(s, alpha=1, lowBeta=-Inf, uppBeta=Inf,
lowOmega=-Inf, uppOmega=Inf)
{
"g" <- function(beta, ...)
{
integrate(function(omega) {likelihood(s=s, omega, beta,
alpha=alpha)}, lowOmega, uppOmega)$value
}
val <- integrate(Vectorize(function(beta) {g(beta, lowOmega=lowOmega,
uppOmega=uppOmega, s=s,
alpha=alpha)}), lowBeta,
uppBeta)$value
return(val)
}
## normalizing constant first integrating beta then omega
`normConst2` <- function(s, x, alpha=1, lowBeta=-Inf, uppBeta=Inf,
lowOmega=-Inf, uppOmega=Inf)
{
"g" <- function(omega, ...)
{
integrate(function(beta) {likelihood(s=s, omega, beta,
alpha=alpha)}, lowBeta, uppBeta)$value
}
val <- integrate(Vectorize(function(omega) {g(omega, lowBeta=lowBeta,
uppBeta=uppBeta, s=s,
alpha=alpha)}), lowOmega,
uppOmega)$value
return(val)
}
## integration limits were choosen by taking the quantiles of beta and
omega
## of a mcmc run
normConst1(s=s, lowBeta=2.8e-06, uppBeta=2.8e-05,
lowOmega=12.5, uppOmega=90)
[1] 5.131021e-162
normConst2(s=s, lowBeta=2.8e-06, uppBeta=2.8e-05,
lowOmega=12.5, uppOmega=90)
[1] 5.246008e-158
## i get normalizing constants with different values
best regards and many thanks
Andreas
Prof Brian Ripley wrote:
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.