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