Passing a function to an RTMB function

[Paul vdb]

What is the most efficient way to allow a user to pass a custom function to the likelihood?

Here is one way I have done it.

extractValue <- function(x, y){
  if(!is.null(x)) return(x)
  return(y)
}

mycars_bayeslm <- function(sd_prior = NULL, beta_priors = NULL){
  dprior_sd <- extractValue(sd_prior, function(x){dgamma(1, 1, log = TRUE)})
  dprior_beta <- extractValue(sd_prior, function(x){sum(dnorm(0, 10, log = TRUE))})
 
  parinit <- list(beta = c(0,0), sd = 5)
  negll <- function(pars){
    getAll(cars, pars)
    -sum(dnorm(dist, beta[1] + beta[2] * speed, sd, log = TRUE)) - dprior_sd(sd) - dprior_beta(beta)
  }
  obj <- MakeADFun(negll, parinit)
  return(obj)
}

obj1 <- mycars_bayeslm()
obj2 <- mycars_bayeslm(sd_prior = function(x){dnorm(x[1], 1, 10, log = TRUE)})

If I am not explicit about returning a scalar, then it does not work, which is minorly annoying.

obj3 <- mycars_bayeslm(sd_prior = function(x){dnorm(x, 1, 10, log = TRUE)})

Does anyone have a better solution? Any suggestions? Thanks so much.

[James Thorson]

Sorry that I haven't looked closely at the code, but I do something similar in `package:ecostate`, search `log_prior` argument here.

--
To post to this group, send email to us...@tmb-project.org. Before posting, please check the wiki and issuetracker at https://github.com/kaskr/adcomp/. Please try to create a simple repeatable example to go with your question (e.g issues 154, 134, 51). Use the issuetracker to report bugs.
---
You received this message because you are subscribed to the Google Groups "TMB Users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to tmb-users+...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/tmb-users/df2b9797-946e-4c05-8aa7-cf508d3a742cn%40googlegroups.com.

[Paul vdb]

For future improvements or others stuck, here is my working solution for including a user defined function in an R package so that it deals with AD overload and ensures that the function returns a scalar. Thanks Jim for the hints, especially passing the AD overload via environment. If someone has a better solution, please let me know.

addUserFunction <- function(fn){
  if(!is.function(fn)){
    return(\(...){
    "c" <- ADoverload("c")
    "[<-" <- ADoverload("[<-")
    0
  })
  }else{
    env_fn = local({
      "c" <- ADoverload("c")
      "[<-" <- ADoverload("[<-")
      environment()
    })
    environment(fn) <- env_fn
    fsum <- function(x, ...){sum(fn(x,...))}
    return(fsum)
  }
}

mycars_bayeslm <- function(sd_prior = NULL, beta_prior = NULL){
  dprior_sd <- addUserFunction(sd_prior)
  dprior_beta <- addUserFunction(beta_prior)
 
  parinit <- list(beta = c(0,0), logsd = log(5))

  negll <- function(pars){
    getAll(cars, pars)

    sd <- exp(logsd)

    -sum(dnorm(dist, beta[1] + beta[2] * speed, sd, log = TRUE)) - dprior_sd(sd) - dprior_beta(beta)
  }
  obj <- MakeADFun(negll, parinit)
  return(obj)
}

obj <- mycars_bayeslm(sd_prior = function(x){dgamma(x, 0.5, 0.5, log = TRUE) + log(x)}, beta_prior = function(x){dnorm(x, 0, 1, log = TRUE)})
obj$fn()
fit <- nlminb(obj$par, obj$fn, obj$gr)
fit$par