[R] Different values for double integral

[Andreas Wittmann]

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