using nimIntegrate to integrate over exponential density with time-varying rate

[gesta...@gmail.com]

Awhile back, in this thread, Perry helped me get started on using the nimIntegrate function. Recently I needed to do a double integral, and got that working too, with all constant parameters.

Now I would like to integrate over an exponential density with a time-varying rate parameter, but I keep getting stuck on several points.  I want the density to be defined by an exponentially decreasing log hazard rate.

For example:

Let the log hazard rate be defined by an intercept and a coefficient controlling exponential decay, something like

b0 <- -5.5
b1 = 4
loghaz <- b0 + b1*exp(-(0:100)/10)
plot(0:100, loghaz)
lam <- exp(loghaz)

Then I want to calculate the density. I understand that the density function for time varying hazard is

f(x) = lam(x)exp(-LAM(x)), 

where LAM(x) is the integral of lam(s) over s = 0 to x.

Lastly, I want to integrate this density function from x = 0 to a. I tried the following, but I don't think it is right, and it doesn't totally run either.  

# inner function
InnerFunc <- nimbleFunction(
run = function(
x = double(1),
theta = double(1)){

f <- exp(theta[1] + theta[2]* exp(-x/10))
    return(dexp(x,f))
    returnType(double(1))
  })
 
# outer integral  
OuterFunc <- nimbleFunction(
run = function(
y = double(1),
theta = double(1)){
return(nimIntegrate(InnerFunc, 0, theta[3], param = theta)[1])
    returnType(double(1))
  })


a = 100
theta = c(-5.5, 4, a)

#this works
nimIntegrate(InnerFunc, 0, theta[3], param = theta)
# this doesn't and produces an error about result of wrong length
nimIntegrate(f = OuterFunc, lower = 0, upper = a, param = theta)

Error in integrate(f, lower, upper, param, subdivisions = subdivisions,  :
  evaluation of function gave a result of wrong length

The first integration works, but I don’t think it is correct because it is not a density with an integrated hazard.

I also tried the following, but it doesn’t quite work either. I am not quite sure how to define the integration limits for the time-varying hazard:

# time varying hazard function
InnerFunc1 <- nimbleFunction(
run = function(
x = double(1),
theta = double(1)){
f2 = exp(theta[1] + theta[2]* exp(-x/10))
    return(f2)
    returnType(double(1))
  })

# integration of hazard function to specify distribution function
# (I'm not sure how to define the upper limit) – seems like it
# should be y, but that generates an error

InnerFunc2 <- nimbleFunction(
run = function(
y = double(1),
theta = double(1)){
lam = nimIntegrate(InnerFunc1, 0, theta[3], param = theta)[1]
return(dexp(x,lam))
    returnType(double(1))
  })

 
# outer integral  
OuterFunc <- nimbleFunction(
run = function(
z = double(1),
theta = double(1)){
return(nimIntegrate(InnerFunc2, 0, theta[3], param = theta)[1])
    returnType(double(1))
  })

# this works, I think
nimIntegrate(InnerFunc1, 0, theta[3], param = theta)[1]
# this doesn't work either
nimIntegrate(f = OuterFunc, lower = 0, upper = age, param = theta)

 Is it possible to do what I am trying? I could do unit cumulative hazard approximations, but execution is much slower than using integrals. Appreciate if anyone can set me straight.

 

[Perry de Valpine]

This reply is only looking at nimble software mechanics, not whether the math you're trying to do is correct.

For nimIntegrate(f, a, b, param), according to help(nimIntegrate), f must be a function (nimbleFunction) that accepts a vector first argument and returns a corresponding vector as output. (It also takes a vector second argument, theta as you have it.) In your first example, it looks like your innerFunc satisfies this protocol. Vector x is input and the result of the calculations you've made will be to return a vector output. However, it doesn't look like outerFunc satisfies this protocol. Vector y is input. But then it calls nimIntegrate and returns the value of that (first element of its return vector), which is a scalar. Hence outerFunc looks like it takes a vector input and returns a scalar, so it does not follow the protocol to be able to called within nimIntegrate. I am wondering if you need to write a for loop to iterate over values of y and call nimIntegrate of innerFunc for each value of y. I'm not following because y is not used at all in outerFunc. Nevertheless I hope this points you in the right direction.

It does seem that if you repeatedly call innerFunc, you will be repeatedly integrating over parts of the same thing, hence computationally it seems possible that you'd be better off writing this manually (including by a discretization method) rather than using nimIntegrate, in order to avoid redundant computations. But I am not sure.

HTH

Perry

--
You received this message because you are subscribed to the Google Groups "nimble-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to nimble-users...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/nimble-users/0092b556-934a-442a-83c9-c9dfc33a21efn%40googlegroups.com.

[gesta...@gmail.com]

Thanks Perry. I am having a really hard time wrapping my head around this, so I am going to turn to a discretized version which I do know how to implement. If I ever come back to this and solve it, I'll post a solution here. 

[Colton Padilla]

Hey gesta,

Would you be willing to share how you got that double integral working? I am actually working on a very similar problem that is going to involve a double integral and I am trying to figure it out.