Calculating log_like for WAIC/LOO with user defined functions

Modeling
, , ,
1

Hi folks,

I’m very new to Stan, and I am currently using capture-recapture models from Kery & Schaub’s Bayesian Population Analysis to estimate population abundance. I’ve used many of these models for several years in JAGS and am quite comfortable with them, but I’m a complete novice when it comes to Stan syntax. With some much needed help from Hiroki Ito, I’ve been modifying his Stan conversions from the BPA book: https://github.com/stan-dev/example-models/blob/master/BPA/Ch.10/js_super_indran.stan . Now I’ve got four different model formulations that I’ve run on my dataset, and after about ten days of computer time, I realized for the first time that Stan does not automatically compute DIC or the metrics necessary for calculating WAIC, LOO, etc. for model comparison. I’m trying to figure out how to calculate a log_lik using a user defined function, and am hitting a wall (no doubt due to my Stan ignorance).

Following this example, I tried adding the following code to the generated quantities block:

  // Log Likelihood for WAIC:
{
  real log_lik;

  log_lik = js_super_lpn(y, first, last, p, phi, psi, nu, chi);

}

But I receive the following error:

SYNTAX ERROR, MESSAGE(S) FROM PARSER:

Sampling statements (~) and increment_log_prob() are
only allowed in the model block or lp functions.

The Stan model code is almost identical to the BPA-converted code above, minus individual effects and with the user defined function name changed to js_super_lpn (because I learned from a previous error message that functions ending in _lp cannot go in the generated quantities block). I’ve attached my code here as well.

Many thanks in advance for your input

js_super_ll.stan (8.6 KB)
.

Cheers,
Josh

2

Hi Josh,

I think you’ll need to modify the js_super_lp function so that instead of incrementing the log density with the likelihood (either via statements like 1 ~ bernoulli(psi) or target += log_sum_exp(lp) from within the function, the function instead increments elements of a vector with n_ind elements that contains the log likelihood of the capture histories (maybe call this vector log_lik), and then returns that vector at the end of the function. Each element in this vector contains the log likelihood for each individual.

If you add the vector log_lik to the list of parameters to collect via the par argument to rstan::stan(), then you will be able to use loo::extract_log_lik() etc.

Roughly this would look like js_super_indran2.stan (8.0 KB)

I have not checked the results from this model with results from the original implementation, but you probably should before investing another 10 days of compute time.

1 Like
3

Hi Max,

Thanks a ton for your quick reply and assistance. I think I understand what’s happening here—in the original model, it’s assessing the likelihood of each capture using 1 ~ bernoulli statements, whereas in your adjustment it is explicitly calculating the bernoulli_lmpf for each individual and aggregating that into log_lik individual by individual (correct me if I’m wrong?). That allows for the explicit calculation of a log likelihood, instead of it just running in the background and not being recorded.

Hiroki Ito made very similar changes when he was helping me add in a sex effect with a mix of observed and unknown sexes. Adding the sex effect roughly doubled the population abundance, and now I’m curious to see if that was truly just the addition of sex or if that had something to do with the way the likelihood was calculated for each individual. Hopefully it’s the former, and this will be a good test.

I modified your code to fit my specific case and it’s running well. I should know in a couple of days if the results are the same as the previous model, so I’ll report back then.

Many thanks again for your help!

Cheers,
Josh

4

The new version surprisingly ran faster than the original (…?) and produced virtually the same results. And I was able to calculate waic and looic using the loo package. Many thanks again for your assistance, Max!

Cheers,
Josh

1 Like
5

No problem! Glad it seems to be working. I don’t know why it might be more efficient to collect these into a vector and then increment the log probability once with the sum, but maybe one of the devs has an idea.

6

There’s a long arXiv paper explaining the autodiff, but the short answer is that we can create an expression graph that looks like (+ x1 x2 x3 ... xN) that has one node and N edges rather than (+ x1 (+ x2 (+ x3 + ... + xN))). This is much faster to traverse and only requires a single virtual function call (which then loops over the N) rather than a bunch of virtual functions, one for each addition. It’s also much tighter in space, which has all kinds of good side effects for modern computer speed.