Functions to debug and profile stan models in R

I’ve found that sometimes I have simple errors in my code that can be hard to debug from the error thrown in the sampling command. Additionally, it can be hard to profile the code using sampling since there is a lot going on with the HMC sampler.

It’s not really clear how to get the log prob value from the stan model, but I was able to come up with this. I was following the rstan thread, and wanted to share this here in case anyone wants to use it, and also in hopes that this sort of thing would be easier in the new rstan version.

checkStanModel <- function(stan.model){
    x <- new(stan.model@mk_cppmodule(stan.model),out.data,rstan:::grab_cxxfun(stan.model@dso))
    x$log_prob(rep(.1,x$num_pars_unconstrained()),adjust_transform=TRUE,gradient=FALSE)
}
profileStanModel <- function(stan.model,n=1000){
    x <- new(stan.model@mk_cppmodule(stan.model),out.data,rstan:::grab_cxxfun(stan.model@dso))   
    return(system.time(replicate(n,x$log_prob(runif(x$num_pars_unconstrained(),-2,2),
                                              adjust_transform=TRUE,gradient=TRUE))))
}
2 Likes

There’s a built in for evaluating log density from a compiled model with data. It’s log_prob and grad_log_prob if you need gradient. That’s on the unconstrained scale, so you may need to transform from constrained to unconstrained.

It felt like there must be a better way… How do I get the compiled model with the data? I couldn’t figure out a way to do it without calling sampling.

I don’t think it is possible, though you could ask @bgoodri or @jonah. One of my main goals for RStan and PyStan 3 is to fix this so that there is an object for compiled program plus data (basically, the posterior density object).

A quick update on this since I just linked here I realized the code does not work on the latest version of RStan. This is what I am using now.

getCallableStan <- function(stan.model){
  if(packageVersion("rstan")>=2.16)
    new(stan.model@mk_cppmodule(stan.model),out.data,0L,rstan:::grab_cxxfun(stan.model@dso))
  else
    new(stan.model@mk_cppmodule(stan.model),out.data,rstan:::grab_cxxfun(stan.model@dso))    
}
checkStanModel <- function(stan.model){
  x <- getCallableStan(stan.model)
  x$log_prob(rep(.1,x$num_pars_unconstrained()),adjust_transform=TRUE,gradient=FALSE)
}
profileStanModel <- function(stan.model,n=1000){
  set.seed(1)
  x <- getCallableStan(stan.model)  
  system.time(replicate(n,x$log_prob(runif(x$num_pars_unconstrained(),-2,2),
                                     adjust_transform=TRUE,gradient=TRUE)))
}

Yeah there’s currently no single object with both the model and data. Definitely something to change for the next major release.