Improving speed for ecological model

[Cole Monnahan]

I've been comparing the speed of Stan vs JAGS on some ecological models. In reading the book 

Korner-Nievergelt et al (2015). Bayesian Data Analysis in Ecology Using Linear Models with R, BUGS, and Stan. (link to website with code)

I ran across example 14.5, which is a mark-recapture model that uses a "trick" to avoid discrete latent parameters by integrating them out with a while loop (the parameter chi, which represents the probability of never seeing the animal again). The authors found it was substantially more efficient than WinBUGS. However, when I coded it myself in JAGS I incorporated this trick instead of modeling the discrete states directly, making the parameterization identical. And to my surprise JAGS was more efficient (minESS/time), albeit not that much. To me this suggests that it was removing the need to track the states directly that was the major difference.

I was quite surprised by that, and wondering if there are any computational improvements that can be made to the Stan model.

I looked through the manual and thought maybe trying to use ragged arrays. The variable `last` represents the last time period the animal was seen, effectively determining the size of the loop over time for that animal. But it is not clear to me how to implement them, or even if that is something worth trying in this case. I also considered using the "ones trick" from BUGS so that the likelihood statements could be vectorized, but wasn't sure how to wrangle the code to accomplish this. 

Does anyone have any thoughts? Or is this just not the type of model that Stan would be expected to outperform JAGS?

[Andrew Gelman]

Hi, can you include the exact Stan and Jags models you used?

A

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Cole Monnahan]

Sure, if you want to run it you can use the attached model files, and data saved as an RDS file. I should note that to make the two models match I had to alter the originally linked Stan model. In particular I used the 'ones trick', and did some preprocessing of the data since JAGS doesn't have a transformed parameter section. This is mostly because JAGS can't do the while loop (at least as far as I could find), so I had to work around that. Also note that this requires JAGS 4.2.0 and rjags 4.4. The rjags version isn't up on CRAN but can be found here:  https://sourceforge.net/projects/mcmc-jags/files/rjags/4/ and installed manually.  I was really curious to see the differences on identical parameterizations, hence all the extra leg work. The Stan model runs about the same as the original.

Thanks for taking a look.

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/xVKzZGexHEU/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Bob Carpenter]

Thanks for doing the legwork and writing in.

How are you calculating n_eff? Stan is much more conservative
than JAGS (uses cross-chain information and also longer-term correlations)
and there can be huge differences in estimated effective sample size.

How big is the data set? You'll find that HMC has more
of an advantage with more dimensions.

Marginalizing out parameters almost always helps with efficiency. Especially
discrete ones that are difficult to mix with Gibbs sampling.

The Stan code isn't very efficient as written.

for(j in 1:4){
b0[j]~normal(0, 10);
b1[j]~normal(0, 10);
}

Vectorize to

b0 ~ normal(0, 10);

Same for the other distributions.

This

sigmayearphi~uniform(0,3);

doesn't match its declaration, which should include the upper bound:

real<lower=0, upper=3> sigmayearphi;

The model should define a finite log density for any parameter
meeting the constraints. With only lower=0, sigmayearphi = 4 satisfies
the constraint, but produces a -infinity log density.

You can pull the first line above the loop:

for (k in 1:last[i]) {
if(k>1) increment_log_prob(log(phi[i, k-1]));

and just use

increment_log_prob(log(phi[i, 1:(last[i] - 1)]));

Those were the things that jumped out at me, but there may be
more if you can find expressions that are repeated --- you don't
want to every apply log or inv_logit or any other expensive function
to the same argument twice.

- Bob

> On May 3, 2016, at 6:22 PM, Cole Monnahan <mon...@uw.edu> wrote:
>

[Andrew Gelman]

Let’s get rid of the bound and get rid of the uniform(0,3) prior. That’s just bad news, both statistically and possibly computationally.
A

[Andrew Gelman]

P.S. It will be good to resolve this one. If we _can_ make it faster with these steps, it will be a good little case study. I think a few case studies like this (including features such as getting rid of the uniform prior) might communicate some of these ideas better than a million replies on stan-users. Our case studies typically are so clean that people don’t reallize all the traps they can fall into.

> On May 3, 2016, at 10:35 PM, Bob Carpenter <ca...@alias-i.com> wrote:
>

[Cole Monnahan]

Hi Bob, thanks for the reply. Unfortunately I had forgotten how much I changed from the original model, and you'll find that I've vectorized where possible in the previous email which had model files, data, inits, and a script to setup the runs. Sorry for the confusion, please take a look at the attached.

I mentioned on here a few months ago that I was doing some testing of JAGS and Stan, initially out of curiosity and as a check on my understanding. But I think it would be an illustrative paper for for many ecologists: an intuitive illustration of HMC and some practical speed tests vs the status quo JAGS. If this forum is a good place to discuss parameterizations and such then I might post a few other models.

. 

A few inline thoughts  below.

On Tue, May 3, 2016 at 7:35 PM, Bob Carpenter <ca...@alias-i.com> wrote:

Thanks for doing the legwork and writing in.

How are you calculating n_eff?  Stan is much more conservative
than JAGS (uses cross-chain information and also longer-term correlations)
and there can be huge differences in estimated effective sample size.

I'm using rstan::monitor() for both Stan and JAGS output. I've done this for a bunch of models and compare it to CODA estimates, and indeed it can be very different and the new way is almost always lower ESS.  

How big is the data set?  You'll find that HMC has more
of an advantage with more dimensions.

Marginalizing out parameters almost always helps with efficiency.  Especially
discrete ones that are difficult to mix with Gibbs sampling.

JAGS isn't using Gibbs sampling here. It's using slice sampling (base::RealSlicer for all 177 parameters). I find that most of the ecology models I test typically use slice sampling, which leads to unexpected things like JAGS taking longer per iteration but having a higher ESS than Stan (the exact opposite of what is expected!), in some cases. So far I've discovered that it's hard to predict which platform is more efficient.

 

The Stan code isn't very efficient as written.

    for(j in 1:4){
      b0[j]~normal(0, 10);
      b1[j]~normal(0, 10);
      }

Vectorize to

  b0 ~ normal(0, 10); 

Same for the other distributions.

As above, my fault for linking to an outdated file. I've done this.

 

This

  sigmayearphi~uniform(0,3);

doesn't match its declaration, which should include the upper bound:

  real<lower=0, upper=3> sigmayearphi;

The model should define a finite log density for any parameter
meeting the constraints.  With only lower=0, sigmayearphi = 4 satisfies
the constraint, but produces a -infinity log density.

Yes I dropped that. I saw the discussion on eliminating the uniform prior density a while back, so took it out.  

You can pull the first line above the loop:

 for (k in 1:last[i]) {
        if(k>1) increment_log_prob(log(phi[i, k-1]));

and just use

  increment_log_prob(log(phi[i, 1:(last[i] - 1)]));

Those were the things that jumped out at me, but there may be
more if you can find expressions that are repeated --- you don't
want to every apply log or inv_logit or any other expensive function
to the same argument twice.

- Bob

> On May 3, 2016, at 6:22 PM, Cole Monnahan <mon...@uw.edu> wrote:
>
> I've been comparing the speed of Stan vs JAGS on some ecological models. In reading the book
>
> Korner-Nievergelt et al (2015). Bayesian Data Analysis in Ecology Using Linear Models with R, BUGS, and Stan. (link to website with code)
>
> I ran across example 14.5, which is a mark-recapture model that uses a "trick" to avoid discrete latent parameters by integrating them out with a while loop (the parameter chi, which represents the probability of never seeing the animal again). The authors found it was substantially more efficient than WinBUGS. However, when I coded it myself in JAGS I incorporated this trick instead of modeling the discrete states directly, making the parameterization identical. And to my surprise JAGS was more efficient (minESS/time), albeit not that much. To me this suggests that it was removing the need to track the states directly that was the major difference.
>
> I was quite surprised by that, and wondering if there are any computational improvements that can be made to the Stan model.
>
> I looked through the manual and thought maybe trying to use ragged arrays. The variable `last` represents the last time period the animal was seen, effectively determining the size of the loop over time for that animal. But it is not clear to me how to implement them, or even if that is something worth trying in this case. I also considered using the "ones trick" from BUGS so that the likelihood statements could be vectorized, but wasn't sure how to wrangle the code to accomplish this.
>
> Does anyone have any thoughts? Or is this just not the type of model that Stan would be expected to outperform JAGS?
>
> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

--

You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/xVKzZGexHEU/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Bob Carpenter]

I meant the general form that provided the "G" in "JAGS",
that is something that draws one parameter at a time conditional
on the others.

Did you provide a new prior for the scale parameters?

Would you mind sharing

* the Stan program
* the JAGS program
* the data
* the MCMC settings you used for comparison
* initialization if any

Thanks.

- Bob

[Cole Monnahan]

Bob: attached are the files that should get you going. Also note that this requires JAGS 4.2.0 and rjags 4.4. The rjags version isn't up on CRAN but can be found here:  https://sourceforge.net/projects/mcmc-jags/files/rjags/4/ and installed manually. I've just been using the defaults in Stan for now: 

adapt_delta=.8

max_treedepth=10

metric='diag_e'

warmup=N/2

and adaptation of stepsize turned on.

I've run them with Niter=200,000. Here's a little table of results.

model	Npar	max.cor	JAGS	Stan	Stan/JAGS
swallows_nc	177	0.800	0.421	0.098	0.230
swallows	177	0.889	0.570	0.063	0.110

where the top row is with non-centered random effects and the bottom row centered. 'max.cor' is just the maximum correlation from draws of an independent sampler. The JAGS and Stan columns are minESS/time, and the last is the ratio. Stan is about 5-10 times slower. This is averaged over 10 random seeds and initial values. I've stored summary outputs about adaptation if that helps.

In some earlier tests I experimented with setting adapt_delta across a gradient (.5 to .99) and found no major improvement in efficiency, but admittedly that was with fewer iterations.

Thanks,

Cole

[Michael Betancourt]

Did you compare the actual posterior distributions?  Where

they consistent across the different algorithms?

You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

<data.RDS><run_model.R><swallows_nc.jags><swallows_nc.stan><sims.ind.RDS>

[Cole Monnahan]

Michael: yes I did this by running them both with a high thinning rate, checking Rhat and ESS, and then plotting marginal QQ plots. Attached are those two plots, the second is a multipage PDF with groups of parameters.

It was enough to convince me that the two posteriors were the same, and that sub-sampling 10 samples from them would be effectively independent, for use in rerunning with default with those as initial values . Assuming convergence, of course.

I've setup some R code that does all this work, given two model files, data, and whatever inputs for Stan you want (metric, adapt_delta, etc.). My idea is to plug some representative ecological models in to see how they do.

[Daniel Lee]

Hi Cole,

Mind calculating ESS for Stan without thinning since you've already got the draws? And maybe for JAGS too. 

Daniel

<model_comparison_qqplots.pdf>

<model_comparison_convergence.png>

[Cole Monnahan]

Daniel: I forgot to mention earlier, there's no thinning for either JAGS or Stan in what I'm doing. 

[Cole Monnahan]

Daniel: Sorry I meant for the table I sent previously. I don't have the unthinned long chains since thinning was done by the functions... stan(..,. thin=500) etc. When I do the actual ESS comparisons I don't do any thinning, since each Stan draw contains more information. -Cole

[Bob Carpenter]

Great question, Daniel!

In general, thinning loses information (lower n_eff)
compared to not thinning.

The reason Daniel's asking here is that thin=500 will mask long
autocorrelation times. Suppose we have N iterations.
If one sampler has N / n_eff = 10 and the other has N / n_eff=100,
then thin = 500 is going to penalize the former sampler's
efficiency by a factor of 10.

So I'd be curious as to what the comparison looks like
for N=500, N=1000, and N=2000 without thinning. What we've
found in almost every model we've carefully compared is
that Stan's faster when measured on a time per n_eff basis.

- Bob

[Cole Monnahan]

Hi Bob (and Daniel), let me clarify a bit. I am only doing a thinning rate of 500 to verify that the JAGS and Stan models have the same posterior (Michael asked if I had done this, so I posted those plots). In this case I don't care at all about comparing efficiencies -- I'm simply trying to make sure it's valid to proceed with that stage of testing. Think of it as a debugging step. 

All of my efficiency comparisons don't use thinning at all: all parameters are saved. Your example is exactly why I don't thin: it will penalize the more efficient sampler.

I have tried doing this with fewer iterations, and then (following advice I saw somewhere on here) doubled the iterations until at least it appears that Stan is adapting nicely. If anything this is unfair to JAGS because it wouldn't need as long of a warmup, but to keep things simple I use the same iterations and warmup (N/2) periods for both platforms.

I realize that it is difficult to do such comparisons fairly, as there are many moving pieces. I hope what I'm currently doing is fairly representative of what an end user would experience, particularly a new one migrating over from a BUGS variant.

Has anyone thought about trying to vectorize the likelihood calculations using a ragged array? That seemed to be the most obvious place to look for this model...

[Bob Carpenter]

Thanks for the clarification. That's exactly how we've
been doing evaluations ourselves. I appreciate the amount
of effort you're putting in to do this right.

I'm surprised that JAGS converges in fewer iterations
than Stan does. Is the newer version very different from
older versions?

- Bob

[Cole Monnahan]

I'm not a JAGS expert, so I won't comment on anything too deep. However, when I tried to recreate the step of marginalizing out the latent parameters, I found that only the newest version had the indexing capabilities I needed. Hence the dependence on that newest version. This is not the only model that JAGS outperfms, so maybe something has changed internally to make it faster? -Cole

[Andrew Gelman]

Interesting, we should try to understand this one.
A