can't get nimIntegrate to work as expected

[gesta...@gmail.com]

I have an R function that I am experimenting with to integrate survival over possible infection times, and am interested in using in a nimble model. It seems to work as expected as an R function, but when I try to do the same thing using nimble functions, I get an error about a missing argument even though I think I am providing a value for the argument. I am likely missing something very simple, but I have puzzled over this for awhile and would appreciate some additional insight.

Below is the R code that works, and further below the nimble functions that I can’t get to work. 

################
# R function for integrating over survival and infection
################
prob.func <- function(lambdasc, lambdaic, lambdafc, e){
   out <- rep(0, length(e))
   integrand <- function(x,lambdaf, lambdas, lambdai, a, b){
    out <- dexp(x-a,lambdaf)*(1-pexp(x-a,lambdas))*((1-pexp(b-x,lambdai)))
    return(as.vector(out))
  }
 for(i in 1:length(e)){
      #Prep values
      sum.lambdasc <- lambdasc+lambdafc
     
      #Calculate Probability
      tempa <- integrate(integrand, lower=0, upper=e[i], lambdaf = lambdafc, lambdai = lambdaic,
                         lambdas = lambdasc, a = 0, b = e[i])
     
      #return log probability for likelihood  
      out[i] <-log(tempa$value + (1-pexp(e[i], sum.lambdasc)))
  }
  return(out)
}
################
# Test
################
Ssurv <- 0.95
FOI <- 0.25  
Isurv <- 0.5  

#Weekly survival & infection hazards
lambdasc <- -log(Ssurv)/52
lambdaic <- -log(Isurv)/52
lambdafc <- -log(1-FOI)/52

# vector of entry times
e <- round(runif(1000,1,nyears*52))

test2 <- prob.func(lambdasc, lambdaic, lambdafc, e)
sum(test2)

#######################################################
###  nimble functions
#######################################################
  integrand <- nimble::nimbleFunction(
run = function(
x = double(1),
lambdasc = double(1),
lambdaic = double(1),
lambdafc = double(1),
a = double(1),
b = double(1)){

   return(dexp(x-a[1],lambdafc[1])*(1-pexp(x-a[1],lambdasc[1]))*((1-pexp(b[1]-x,lambdaic[1]))))
   returnType = double(1)
  })

prob.func.nimble <- nimble::nimbleFunction(
run = function(
lambdasc = double(0),
lambdaic = double(0),
lambdafc = double(0),
lower = double(1),
    upper = integer(1)){

 
  out <- nimNumeric(length(upper), init = FALSE)
 
 for(i in 1:length(e)){
      #Prep values
      sum.lambdasc <- lambdasc+lambdafc
  a = lower
  b = upper[i]
      param=c(lambdasc,lambdaic,lambdafc, a, b)
     
  #Calculate Probability
      tempa <- integrate(f = integrand, lower = a, upper = b, param = param)
     
      #return log probability for likelihood  
      out[i] <-log(tempa[1] + (1-pexp(b[i] - a[1], sum.lambdasc)))

  }
  return(out)
  returnType = double(1)
})

test2 <-prob.func.nimble(lambdasc, lambdaic, lambdafc, lower = 0, upper = e)
# I get the following error
"Error in dexp(x - a[1], lambdafc[1]) :
  argument "a" is missing, with no default"

[Perry de Valpine]

Thanks for the question.

integrand must be defined to take two arguments. Yours takes 6 arguments. The second argument should be a vector that you will need to unpack internally into the 5 pieces you need (i.e. to invert the way you packed them up as param=c(lambdasc,lambdaic,lambdafc, a, b), which makes sense).

The reason for the decreased flexibility of nimble's integrate (aka nimInterate) vs R's integrate is that it is not so easy in generated C++ (from compileNimble) to have functions of arbitrary numbers of arguments. (It could be done, but the current system is simpler than that.) So, instead, there is a single additional argument to nimIntegrate and the integrand function must be coded to unpack it internally. If you need lengths of components to unpack them (e.g. the lengths of a or b, which may not be determined easily when they are concatenated in a vector), you can include those lengths as elements in the vector as well.

I hope that helps. There is an attempt at explaining this in help(nimIntegrate). (integrate in a nimbleFunction becomes nimIntegrate.)

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 on the web visit https://groups.google.com/d/msgid/nimble-users/29e27c15-4a07-4785-bb97-bbfd84e5da6cn%40googlegroups.com.

[gesta...@gmail.com]

Thanks Perry. I got it now, and it is working as expected. I figured I was missing something easy.

To complete the thread,  I changed he integrand to the following, then defined theta each time in the loop where I calculated the integral.

  integrand <- nimble::nimbleFunction(
run = function(
x = double(1),

theta = double(1)){

   return(dexp(x-theta[4],theta[3])*(1-pexp(x-theta[4],theta[1]))*((1-pexp(theta[5]-x,theta[2]))))
   returnType = double(1)
  })