Translation of examples in Hill, Gelman books into Stan code

[Ben Lambert]

Hi,

I am new to this group, so forgive me if this is the incorrect place to post this. I wanted to help by translating some existing examples in the aforementioned books (and maybe some others, for example some in Lee's book on 'Bayesian cognitive modelling') into working Stan code. I am newish to Stan; having used BUGS previously, but have been using it to successfully fit some hierarchical models that I previously struggled with, over the last couple of months. It may take me a bit longer than more experienced Stan developers to code some examples, but am hopeful I can contribute!

I am aware that there is already considerable effort to 'translate'/make working Stan code for the books I mentioned, and so wanted to ask - has work begun on writing Stan code for Gelman, Hill and other books? I can see that the majority of BUGS book examples are done already. Also, where would be most helpful for me to begin?

Best,

Ben Lambert

[Daniel Lee]

Hi Ben,

This is a good place to start the discussion. We've started coding up models here:
  https://github.com/stan-dev/example-models
Some of the models from Gelman and Hill have been implemented under the folder ARM.

More below.

On Sun, May 10, 2015 at 2:11 PM, Ben Lambert <ben.c....@googlemail.com> wrote:

Hi,

I am new to this group, so forgive me if this is the incorrect place to post this. I wanted to help by translating some existing examples in the aforementioned books (and maybe some others, for example some in Lee's book on 'Bayesian cognitive modelling') into working Stan code. I am newish to Stan; having used BUGS previously, but have been using it to successfully fit some hierarchical models that I previously struggled with, over the last couple of months. It may take me a bit longer than more experienced Stan developers to code some examples, but am hopeful I can contribute!

I'm sure you'll be contributing in short order!

I am aware that there is already considerable effort to 'translate'/make working Stan code for the books I mentioned, and so wanted to ask - has work begun on writing Stan code for Gelman, Hill and other books? I can see that the majority of BUGS book examples are done already. Also, where would be most helpful for me to begin?

Some background:

The example-models repository has just been a place where we have contributed models, without checking for quality. I'm sure not everything there even compiles under Stan. That said, some of them implementations might be the best way to do it. We just haven't systematically verified anything.

Here are some things I think would be helpful:

* discuss how users would use the example-models repository and how we can make that experience better for them. I think one of the things we can do is to provide text descriptions about the models in a standard format.

* create a standard layout for each model; this should include the data and how to run the model

* consistent formatting, perhaps syntax highlighting in the description?

* apply best practices for model building

* place a little bit of process on the repository so everything that gets in meets the standards that we set

Daniel

[Ben Lambert]

Hi Dan,

Thanks for that - much appreciated.

Perhaps a good place to start in which case would be for me to go through the various examples that are done, and firstly check that they compile in Stan. I can then have a think about how perhaps we can present them in a unified manner. What do you think?

Best,

Ben

On 10 May 2015, at 19:35, Daniel Lee <bea...@alum.mit.edu> wrote:

Hi Ben,

This is a good place to start the discussion. We've started coding up models here:
  https://github.com/stan-dev/example-models
Some of the models from Gelman and Hill have been implemented under the folder ARM.

More below.

On Sun, May 10, 2015 at 2:11 PM, Ben Lambert <ben.c....@googlemail.com> wrote:

Hi,

I am new to this group, so forgive me if this is the incorrect place to post this. I wanted to help by translating some existing examples in the aforementioned books (and maybe some others, for example some in Lee's book on 'Bayesian cognitive modelling') into working Stan code. I am newish to Stan; having used BUGS previously, but have been using it to successfully fit some hierarchical models that I previously struggled with, over the last couple of months. It may take me a bit longer than more experienced Stan developers to code some examples, but am hopeful I can contribute!

I'm sure you'll be contributing in short order!

I am aware that there is already considerable effort to 'translate'/make working Stan code for the books I mentioned, and so wanted to ask - has work begun on writing Stan code for Gelman, Hill and other books? I can see that the majority of BUGS book examples are done already. Also, where would be most helpful for me to begin?

Some background:

The example-models repository has just been a place where we have contributed models, without checking for quality. I'm sure not everything there even compiles under Stan. That said, some of them implementations might be the best way to do it. We just haven't systematically verified anything.

Here are some things I think would be helpful:

* discuss how users would use the example-models repository and how we can make that experience better for them. I think one of the things we can do is to provide text descriptions about the models in a standard format.

* create a standard layout for each model; this should include the data and how to run the model

* consistent formatting, perhaps syntax highlighting in the description?

* apply best practices for model building

* place a little bit of process on the repository so everything that gets in meets the standards that we set

Daniel

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

[Bob Carpenter]

What we want is someone to take a pass over all the
models and make sure they're up to our current thinking
in terms of priors, coding efficiency, and coding style.

We'd ideally like to have for each BUGS model:

1) direct translation of the BUGS models, preserving priors

2) model with our preferred priors

So while they're mostly there, they're in various stages
of having both versions for each model and (2) being up to our
current coding practices (vectorization, priors, etc.)

For example, the very first example in the vol1 needs to be
updated:

https://github.com/stan-dev/example-models/tree/master/bugs_examples/vol1/blocker

I'm assuming this is a reasonable direct translation (though maybe using the regular
binomial instead of binomial-logit would be more faithful, we're more concerned
about having the same model than the same coding, so algebraic differences like
binomial(inv_logit(alpha)) vs. binomial_logit(alpha) shouldn't matter).

What we want is a "best practices" version with weakly informative priors as
recommended in the regression chapter of the manual.

We'd also love to get simulations going for all these models, but that's
a secondary concern.

For the ARM models, they all need to be tested --- I think many of them
are broken.

- Bob

[Andrew Gelman]

In my reboot of ARM, I’d also like to be adding informative priors, but no need to worry about that right now.

[Michael Betancourt]

I’d actually argue that simulations are paramount!

Recall that we have two problems: is our fit parsimonious with the
given model and is our model parsimonious with the “truth”?
Simulated data and our diagnostics checks the former, PPCs
check the latter.

I think it’d be great if for every example we had:

Step 1:
Simulate data from the model
Fit the model to the simulated data
Check diagnostics (n_divergent is zero for all iterations, tree_depth < max_depth for all iterations, R_hat < 1.1)
— if diagnostics fail improve model with a note explaining the fix
Check fit to exact results
— if fit is poor then improve model with a note explaining the fix

Step 2:
Fit the model to the true data
Check diagnostics (n_divergent is zero for all iterations, tree_depth < max_depth for all iterations, R_hat < 1.1)
— if diagnostics fail improve model with a note explaining the fix
Check posterior predictive checks
— if PPCs are bad improve model with a note explaining the fix

Basically, we want to encourage people to always see if they can recover the results from
simulated data before moving on to evaluating their model with PPCs.

[Bob Carpenter]

Yes, this is how we should be rolling.

I added this to the notes.

- Bob

[Ben Lambert]

Hello all,

I have started to have a go at looking at the Stan translation of the BUGS example 1 (https://github.com/stan-dev/example-models/tree/master/bugs_examples/vol1/blocker), that is in the repository. I have tried to follow Bob's list of example points (https://github.com/stan-dev/example-models/blob/master/README.md), by first simulating data from the model, and fitting the model to the fake data, then checking the diagnostics. Two  of the checks appear ok: the value of Rhat < 1.1 for all parameters, and tree_depth<max_depth for all iterations. However, for the "n_divergent" parameter, I find that for about 10/2000 iterations this parameter is 1. 

Forgive me; this is probably due to my misunderstanding of the problem. I understand from reading the literature that this parameter indicates that the integrator is diverging during warm up (the first 1000 iterations). I have tried various re-parameterisations, as well as changing the priors (I have changed the scale parameter to be a cauchy(0,2.5) and tried using a normal rather than t prior for delta). I have also tried manually specifying a stepsize, to prevent the integrator stepping too far, but I cannot seem to prevent this divergence from occurring. 

Looking at the indices where the divergence occurs, it appears to happen frequently at iterations around: 2, 3, 4, 9, 102, 152, 252, 452, 952. 

I attach the RStan script file, which first generates simulated data, then fits the model to the simulated data. I also attach the Stan file, which has only been changed a bit vs the version on git. A finally attach the real data file.

Any hints/tips would be most welcome. As I have said, I am not an expert at HMC!

Best,

Ben

On Mon, May 11, 2015 at 5:41 PM, Bob Carpenter <ca...@alias-i.com> wrote:

Yes, this is how we should be rolling.

I added this to the notes.

- Bob

> On May 10, 2015, at 6:33 PM, Michael Betancourt <betan...@gmail.com> wrote:
>

[Michael Betancourt]

Firstly, divergences during warmup are ignorable as they usually typically

denote poor steps in the adaptation instead of pathologies in the sampler.

So if the divergences are really occurring only in the first 1000 iterations

then you’re fine.

It’s only when you see divergences in the sampling phase that you have

to worry.  In this case manually specifying a step size will not be sufficient

because the adaptation will rapidly change the initial value (although

setting a lower step size can sometimes speed up adaptation if larger

step sizes lead to awkward parameter values and expensive likelihood

computations).  What you have to do is increase the target acceptance

probability (change “delta” parameter from the default 0.8 to 0.9 or

higher as necessary for removing the divergences).

For more see http://arxiv.org/abs/1312.0906.

<blockerSimulatedData.R><blockerUpdated.stan><blocker.data.R>

[Ben Lambert]

Great - thanks for that; much appreciated.

On 12 May 2015, at 09:01, Michael Betancourt <betan...@gmail.com> wrote:

Firstly, divergences during warmup are ignorable as they usually typically

denote poor steps in the adaptation instead of pathologies in the sampler.

So if the divergences are really occurring only in the first 1000 iterations

then you’re fine.

It’s only when you see divergences in the sampling phase that you have

to worry.  In this case manually specifying a step size will not be sufficient

because the adaptation will rapidly change the initial value (although

setting a lower step size can sometimes speed up adaptation if larger

step sizes lead to awkward parameter values and expensive likelihood

computations).  What you have to do is increase the target acceptance

probability (change “delta” parameter from the default 0.8 to 0.9 or

higher as necessary for removing the divergences).

For more see http://arxiv.org/abs/1312.0906.

On May 12, 2015, at 1:49 AM, Ben Lambert <ben.l...@some.ox.ac.uk> wrote:

<blockerSimulatedData.R><blockerUpdated.stan><blocker.data.R>

[Ben Lambert]

Hello all,

I have finally got some time to get through writing up the first example of the BUGS book. Please find it attached. I think I have carried out all the steps which Bob wrote up on the github page. Before moving on to other examples, I wanted to get your feedback on the attached, as I am sure you will have preferences regarding content and style, as well as what I call the 'preferred' Stan model, which uses the most up-to-date thinking on priors, and coding style. I stopped short of building a beta-binomial model, as I allude to at the end, as I though this would be overkill, however, if you think this would be worthwhile, then I can look into it.

Once we have an agreed style in place, I don't think it will take as long to go through the various other examples, I hope.

Look forward to hearing your feedback - I'm sure there'll be plenty!

Best,

Ben

On Tue, May 12, 2015 at 12:41 PM, Ben Lambert <ben.l...@some.ox.ac.uk> wrote:

Great - thanks for that; much appreciated.

On 12 May 2015, at 09:01, Michael Betancourt <betan...@gmail.com> wrote:

Firstly, divergences during warmup are ignorable as they usually typically

denote poor steps in the adaptation instead of pathologies in the sampler.

So if the divergences are really occurring only in the first 1000 iterations

then you’re fine.

It’s only when you see divergences in the sampling phase that you have

to worry.  In this case manually specifying a step size will not be sufficient

because the adaptation will rapidly change the initial value (although

setting a lower step size can sometimes speed up adaptation if larger

step sizes lead to awkward parameter values and expensive likelihood

computations).  What you have to do is increase the target acceptance

probability (change “delta” parameter from the default 0.8 to 0.9 or

higher as necessary for removing the divergences).

For more see http://arxiv.org/abs/1312.0906.

On May 12, 2015, at 1:49 AM, Ben Lambert <ben.l...@some.ox.ac.uk> wrote:

<blockerSimulatedData.R><blockerUpdated.stan><blocker.data.R>

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

[Bob Carpenter]

Thanks --- that looks awesome. Much more than we were
expecting. Although it'd be great, the key thing is getting
the code for both the BUGS model and "best practices" model right.

Your updated model looks great!

As to the direct translation, it's fine as is, but it doesn't have
to be literally line for line. I think it's OK to translate the model,
not the implementation. So it's OK to use, for example, binomial_logit
in the translation. This is really a judgment call, but I'd rather
err on the side of efficient Stan idioms as long as we literally
translate the model as written to one with the same posterior.
Maybe some bread crumbs in the text on the translation?

Then we want to provide the model we would've used.

Some specific comments:

* Num -> N

* keep the opening curly brace on the same line as the for statement
(we basically follow the Google code style in Stan as well as in
C++).

* spaces after commas (e.g., in gamma(0.001, 0.001))

* space between operators ("1 / sqrt(tau)" rather than "1.0/sqrt(tau)").

* Indentation in code needs to be right. No tabs and everything
indented 2 spaces. You missed the indent on the binomial lines
and all of the generated quantities block. There's an emacs mode
that'll clean this up for you automatically if you get everything
on the right lines. We're not entirely consistent yet in the manual
and elsewhere, but we're aiming for it.

* I think instead of "1 / 1E-6" as your scale, you just want "1e3",
that is, sd = 1 / sqrt(tau).

* I'd just define sigma <- 1 / sqrt(tau) as a transformed parameter
and then use it in the student_t (sampling statement in model and
rng in generated quantities), and eliminate sigma from generated quantities

* The fake data script is for the original model; I'm not sure if
we want it for both models or just our form?

* It'd be nice if we could do the simulation in Stan, too --- the R
interface isn't really set up to do it very easily, but it's possible
to run for a single iteration with generated quantities and no params and
then pull out structured values --- just some awkward index fiddling.
But then there's also something to be said for standalone fake data
generation.

* always a space after commas, even in R

* Table 3 should be a plot with actual value on the x axis and fitted
value with some interval on the y axis (geom_errorbars or something like
that in ggplot2); it could be faceted by top-level param name

* diagnostic discussion is great, but it should probably be 4.1, not 4.0.3.

* n_{\mathrm{eff}} and n_{\mathrm{leapfrog}} and \hat{R} in LaTeX

* I didn't understand the point of table 4 --- they're comparing different
models and most of the BUGS entries are missing anyway

* The code after table 4 needs indentation

* Instead of


if (simTreatMort[i] > rt[i])
{

indicatorTreat[i] <- 1;
}
else
{

indicatorTreat[i] <- 0;
}

indent as

if (simTreatMort[i] > rt[i]) {
indicatorTreat[i] <- 1;
} else {
indicatorTreat[i] <- 0;
}

but what you really want to do is this:

indicatorTreat[i] <- (simTreatMort[i] > rt[i]);

* $p$-values, I think --- isn't the p usually italicized? I don't
know for sure as I don't do much with them.

* If you cite BDA, it should be BDA3

* Your bib form is messing up capitalization of journals for some reason

* Also, you need our fonts to get the code blocks right --- LaTeX messes
up ~ as a matter of course in most fonts, conflating it with the diacritic
tilde.

- Bob

> <Example1.pdf><blockerSimulatedData.R>

[Andrew Gelman]

Bob:
Tell me the truth here. Did you just include this as a test to see if I’m reading the messages on this list?
A

[Bob Carpenter]

Nope --- it's just that we decided to exactly translate
the BUGS models to help those nervous about our getting
different answers from BUGS sleep more easily. We are also
writing them the way we think they should be written, at least
until your thinking priors evolves further.

The thing to comment on is the original PDF with the
writeup of the model --- that's exactly up your alley.
I already channeled for you and for every table said "graph"!
(Oh, I forgot for the first one, but that should obviously be
a graph, too. It just looked so much like the way the BUGS writeups
start it slipped under my radar.)

- Bob

[Daniel Lee]

If you want to push Bob's buttons, you'd do this:

    gamma (0.001,0.001
)

where the indentation is not aligned with anything above and it uses spaces.

Bob, I think we should seriously think about a Stan linter or formatter.

[Bob Carpenter]

> On Jun 18, 2015, at 8:50 PM, Daniel Lee <bea...@alum.mit.edu> wrote:
>
> If you want to push Bob's buttons, you'd do this:
> gamma (0.001,0.001
> )
>
> where the indentation is not aligned with anything above and it uses spaces.
>
>
> Bob, I think we should seriously think about a Stan linter or formatter.

I could probably do that from the parsed model, but it
might be easier to do in an ad-hoc way from the outside.

- Bob

[Ben Lambert]

Hello all again,

I am now attaching a reworked write up of example 1, with graphs replacing tables etc. 

I have managed to, I hope, do everything that you mentioned before Bob, with the exception of the tildes within the code blocks, and the Bibliography. Which latex package do you use to write code in the Stan manuals? I am also attaching the .tex file, so that you can edit if you wish.

Let me know your feedback, as it would be good to get the style down before we move on to the other examples.

Best,

Ben

On Thu, Jun 18, 2015 at 9:07 PM, Bob Carpenter <ca...@alias-i.com> wrote:

Thanks --- that looks awesome.  Much more than we were
expecting.  Although it'd be great, the key thing is getting
the code for both the BUGS model and "best practices" model right.

Your updated model looks great!

As to the direct translation, it's fine as is, but it doesn't have
to be literally line for line.  I think it's OK to translate the model,
not the implementation.  So it's OK to use, for example, binomial_logit
in the translation.  This is really a judgment call, but I'd rather
err on the side of efficient Stan idioms as long as we literally
translate the model as written to one with the same posterior.
Maybe some bread crumbs in the text on the translation?

Then we want to provide the model we would've used.

Some specific comments:

* Num -> N

* keep the opening curly brace on the same line as the for statement
  (we basically follow the Google code style in Stan as well as in
   C++).

* spaces after commas (e.g., in gamma(0.001, 0.001))

> <Example1.pdf><blockerSimulatedData.R>

[Bob Carpenter]

Lucida, to go with the Lucida fonts. The only problem is
that they're proprietary and Bigelow and Holmes charge
for them and we don't have a site license to share the licenses.

I don't know if other fonts manage to fix some of the issues
with LaTeX that Lucida fix.

- Bob

> <Example1.pdf><Example1.tex>

[Bob Carpenter]

Thanks --- that looks fantastic. Are you generating the LaTeX
through knitr?

I like the plot with all the parameter values, but this style's going
to get problematic when the parameters aren't all on the same scale.
If I were doing it, I'd facet the different named parameters through
ggplot2. But this is already way more than we need in this context,
so I'm not suggesting redoing it!

In general, I think 90% intervals are OK --- the problem is that
the larger the interval, the noisier the endpoint estimates are.

It's great that you plotted that we get the same result as BUGS.
Could you indicate if it's OpenBUGS or WinBUGS and which version?

The overlplotting of orange and blue (great color choice, by the
way) in the BUGS vs. Stan is a bit hard to read, but I don't think
it's a big deal. What is an issue is that you probably want to be
comparing means +/- MCMC std error, as those are the things that should
converge to the same value if everything's working.

And to compare apples and apples, you need to do both with Stan's
method of computing n_eff. I have some code to convert R2jags
output to something Stan can process.

The bibliography's still messed up in terms of capitalization.
(And "2x 2" rather than "2 x 2".)

You want to use "a---b" for em-dashes to separate text in
American English typesetting (you get an extra space in the UK,
but you also want to use the em-dash).

The values 1 and 0 shouldn't be quoted. A quote usually
indicates a character rather than the value.

Spaces, so it's not [0.025,0.975], but [0.025, 0.975].

And remember, 95% intervals are going to be noisier than,
say 80% intervals. In both cases, you can calibrate by eye.

Sorry for all the typesetting comments --- I'm a nut about it,
but it's all minor stuff compared to the content.

- Bob

P.S. Some R to fit a JAGS model (for an IRT example) and then
pull out the effective sample sizes. It's not complete, as
there are dependencies on things like num_params, but I hope it's
clear what's going on (the relevant bit is the monitor and n_eff
extraction, which Ben showed me how to do):

start_time_jags <- proc.time();

fit_jags <- jags(model.file = "irt.jags",
data = list(I=I, J=J, y=y),
init = init_list,
parameters.to.save = c("alpha","beta","delta"),
n.chains = num_chains,
n.iter = 10000,
jags.seed = init_seed);

elapsed_time_jags = proc.time() - start_time_jags;
print(c("jags", elapsed_time_jags));

fit_jags_mcmc <- as.mcmc(fit_jags);
arr <- as.array(fit_jags_mcmc);
arr <- aperm(arr, c(1,3,2));
mon <- monitor(arr, inc_warmup = FALSE);

n_eff_jags <- mon[1:num_params, "n_eff"];
sec_jags <- elapsed_time_jags[3];
n_eff_per_sec_jags <- n_eff_jags / sec_jags;

> <Example1.pdf><Example1.tex>

[Ben Lambert]

Thanks - will respond to your points below.

On Tue, Jun 23, 2015 at 6:23 PM, Bob Carpenter <ca...@alias-i.com> wrote:

Thanks --- that looks fantastic.  Are you generating the LaTeX
through knitr?

No - I haven't been. I've never used this package. I have been making some of the figures outside of R, so that may be a snag?

I like the plot with all the parameter values, but this style's going
to get problematic when the parameters aren't all on the same scale.
If I were doing it, I'd facet the different named parameters through
ggplot2.    But this is already way more than we need in this context,
so I'm not suggesting redoing it!

I agree. I looked at plotting actual vs estimated on a scatterplot, but this looks really messy, and decided that the approach I've taken looks better. However, if you an example from ggplot2 that you could show me, then I can have a look at replicating it...

 

In general, I think 90% intervals are OK --- the problem is that
the larger the interval, the noisier the endpoint estimates are.

Yep, I will redo these with those intervals.

 

It's great that you plotted that we get the same result as BUGS.
Could you indicate if it's OpenBUGS or WinBUGS and which version?

Yes, I will include this information, FYI it's OpenBUGS. version 3.2.3. rev 1012. 

The overlplotting of orange and blue (great color choice, by the

way) in the BUGS vs. Stan is a bit hard to read, but I don't think

it's a big deal.  What is an issue is that you probably want to be
comparing means +/- MCMC std error, as those are the things that should
converge to the same value if everything's working.

Ok. I will redo with the plot with these. 

And to compare apples and apples, you need to do both with Stan's
method of computing n_eff.  I have some code to convert R2jags
output to something Stan can process.

I'm not sure I follow you. Your code makes sense for extracting n_eff from JAGS output, but how will I use this as an input for Stan? Sorry, I think this is somewhat due to my poor knowledge of n_eff; I know that it is an effective sample size essentially measures the number of independent iterations from mcmc, but unsure how to specify this a priori for Stan.

The bibliography's still messed up in terms of capitalization.
(And "2x 2" rather than "2 x 2".

Ok, is there a style that you recommend? At the moment I am just using 'plain'.

 

You want to use "a---b" for em-dashes to separate text in
American English typesetting (you get an extra space in the UK,
but you also want to use the em-dash).

The values 1 and 0 shouldn't be quoted.  A quote usually
indicates a character rather than the value.

Ok. 

Spaces, so it's not [0.025,0.975], but [0.025, 0.975].

Yep, my mistake. Will change.

 

> <Example1.pdf><Example1.tex>

[Bob Carpenter]

> On Jun 23, 2015, at 1:47 PM, Ben Lambert <ben.l...@some.ox.ac.uk> wrote:
>
> Thanks - will respond to your points below.
>
> On Tue, Jun 23, 2015 at 6:23 PM, Bob Carpenter <ca...@alias-i.com> wrote:
> Thanks --- that looks fantastic. Are you generating the LaTeX
> through knitr?
> No - I haven't been. I've never used this package. I have been making some of the figures outside of R, so that may be a snag?

It doesn't have to be R, but it'd be nice to have all the
code used to generate the examples available.

> I like the plot with all the parameter values, but this style's going
> to get problematic when the parameters aren't all on the same scale.
> If I were doing it, I'd facet the different named parameters through
> ggplot2. But this is already way more than we need in this context,
> so I'm not suggesting redoing it!
>
> I agree. I looked at plotting actual vs estimated on a scatterplot, but this looks really messy, and decided that the approach I've taken looks better.

You're using the usual one (though Andrew argues for plotting the other
way --- swapping x and y --- in some cases for reasons I still don't
fully understand having to do with viewing something as a regression).

> However, if you an example from ggplot2 that you could show me, then I can have a look at replicating it...
>
> In general, I think 90% intervals are OK --- the problem is that
> the larger the interval, the noisier the endpoint estimates are.
>
> Yep, I will redo these with those intervals.
>
> It's great that you plotted that we get the same result as BUGS.
> Could you indicate if it's OpenBUGS or WinBUGS and which version?
>
> Yes, I will include this information, FYI it's OpenBUGS. version 3.2.3. rev 1012.

Thanks. I should try to intsall OpenBUGS and see how it compares
to JAGS in terms of speed. So far, all we know is this:

https://martynplummer.wordpress.com/2010/09/20/how-fast-is-jags/

Martyn rescales by dividing time by the mean effective sample size.
We're not sure what the right way to compare is.

> The overlplotting of orange and blue (great color choice, by the
> way) in the BUGS vs. Stan is a bit hard to read, but I don't think
> it's a big deal. What is an issue is that you probably want to be
> comparing means +/- MCMC std error, as those are the things that should
> converge to the same value if everything's working.
>
> Ok. I will redo with the plot with these.
>
> And to compare apples and apples, you need to do both with Stan's
> method of computing n_eff. I have some code to convert R2jags
> output to something Stan can process.
>
> I'm not sure I follow you. Your code makes sense for extracting n_eff from JAGS output, but how will I use this as an input for Stan?

Some of those commands are for rstan --- that was a complete
example of pulling n_eff out for JAGS output --- presumably the
same thing would work for BUGS if they have Coda-compatible output.

> Sorry, I think this is somewhat due to my poor knowledge of n_eff; I know that it is an effective sample size essentially measures the number of independent iterations from mcmc, but unsure how to specify this a priori for Stan.

The point is that we do n_eff estimates more conservatively in
Stan than in packages like Coda, which the R packages for BUGS
and JAGS use. That leads to bigger MCMC std error estimates,
because they're just

mcmc_std_err = posterior_sd / sqrt(n_eff)

But this isn't critical here, so you could just use the BUGS
versions of mcmc std err (if they have them).

Overall, I'm still not sure we should be comparing to BUGS as
part of these writeups. One, we don't want to give the impression
that BUGS is somehow the "reference" answer --- we trust Stan more
to do the right thing in the tails and for hard models.

On the other hand, part of the motivation is to reassure people
who trust BUGS that Stan can be trusted, too.

> The bibliography's still messed up in terms of capitalization.
> (And "2x 2" rather than "2 x 2".
>
> Ok, is there a style that you recommend? At the moment I am just using 'plain'.

I don't care as long as it capitalizes properly! For
the manual, we use:

\bibliographystyle{apalike}
\makeatletter
\renewcommand\@biblabel[1]{}
\makeatother
\bibliography{../bibtex/all}

I don't have the exact example coded up in ggplot2. What you want
to do is something like facet_wrap(~ parameter) in a data frame with
all the params. Here's an example I did of plotting Stan vs. JAGs
times for all paramters (as opposed to just averaging):

df_stan <- data.frame(n_eff_per_sec = n_eff_per_sec_stan,
parameter = param_names,
mcmc = rep("stan", num_params));

df_jags <- data.frame(n_eff_per_sec = n_eff_per_sec_jags,
parameter = param_names,
mcmc = rep("jags", num_params));

df <- rbind(df_stan, df_jags);

library(ggplot2);
plot <- ggplot(df, aes(x = n_eff_per_sec, fill = mcmc)) +
geom_histogram() +
scale_x_log10() +
facet_wrap(~ parameter) +
ggtitle("Stan vs. JAGS: IRT 2PL, 20 questions (I), 200 students (J)") +
xlab("n_eff / sec") +
ylab("# parameters");

This produces:

[Sebastian Weber]

Hi!

First of all I want to say that this is really great! We need this type really bad and what you produced is really good stuff I think. Some comments:

1. I think it would be good to make the BUGS engine exchangeable, i.e. the user should be able to switch between WinBUGS (the mother of all BUGS), OpenBUGS and JAGS. There should be R packages out there which are of help here.

2. Before you move on to other examples I suggest to think about how all of the example can later on be processed in an automatic fashion, i.e. I am thinking of something similar (but more sophisticated) like the stan_demo function. Basically I would appreciate if one can get all these nice example models in a programmatic way as other people prefer different comparisons.

3. The way I once came up with to compare Stan vs WinBUGS output was the following:
- I declared WinBUGS as reference which is what you want to do as this type of work in anyway only for WinBUGS followers
- Then I standardize each quantity of interest from the models with respect to the median and robust standard deviation (inter-quantile-range divided by 1.34) such that the WinBUGS output and the Stan output becomes roughly standard normal
- Last, I computed the 2.5%, 5% ... 95%, 97.5% quantiles from each program and then checked that the absolute difference is below 0.1

4. The simulation of the fake data should be done somehow with the fixed_param algo from rstan as Bob points out.

Could we place this stuff in some git repo? This makes it easier to track the code & interact. My third point is something which really does not need to be part of this, but I would really appreciate if I could do this type of stuff with this work given that point 2 is taken care of.

Really cool!

Best,
Sebastian

On Tuesday, June 23, 2015 at 8:43:31 PM UTC+2, Bob Carpenter wrote:
> > On Jun 23, 2015, at 1:47 PM, Ben Lambert
> >

> > Thanks - will respond to your points below.
> >
> > On Tue, Jun 23, 2015 at 6:23 PM, Bob Carpenter

[Andrew Gelman]

> On Tuesday, June 23, 2015 at 8:43:31 PM UTC+2, Bob Carpenter wrote:
>>>
>>>
>>> I agree. I looked at plotting actual vs estimated on a scatterplot, but this looks really messy, and decided that the approach I've taken looks better.
>>
>> You're using the usual one (though Andrew argues for plotting the other
>> way --- swapping x and y --- in some cases for reasons I still don't
>> fully understand having to do with viewing something as a regression).

See Figures 8.1 and 8.2 of ARM.

[Bob Carpenter]

> On Jun 23, 2015, at 3:33 PM, Sebastian Weber <sdw....@gmail.com> wrote:
>
> Hi!
>
> First of all I want to say that this is really great! We need this type really bad and what you produced is really good stuff I think. Some comments:
>
> 1. I think it would be good to make the BUGS engine exchangeable, i.e. the user should be able to switch between WinBUGS (the mother of all BUGS), OpenBUGS and JAGS. There should be R packages out there which are of help here.

The languages are not the same, so it would mean
rewriting the models.

> 2. Before you move on to other examples I suggest to think about how all of the example can later on be processed in an automatic fashion, i.e. I am thinking of something similar (but more sophisticated) like the stan_demo function. Basically I would appreciate if one can get all these nice example models in a programmatic way as other people prefer different comparisons.

It wasn't clear to me what you're asking for here. Could you provide
an example?

> 3. The way I once came up with to compare Stan vs WinBUGS output was the following:
> - I declared WinBUGS as reference which is what you want to do as this type of work in anyway only for WinBUGS followers

I'd prefer not to describe WinBUGS as a reference in anything we
write.

> - Then I standardize each quantity of interest from the models with respect to the median and robust standard deviation (inter-quantile-range divided by 1.34) such that the WinBUGS output and the Stan output becomes roughly standard normal

You mean posterior median of a parameter? And what's the inter-quantile
range that you're dividing by 1.34?

I don't think this makes sense from an MCMC perspective, which
gives you guarantees on the means and error from the MCMC CLT.

> - Last, I computed the 2.5%, 5% ... 95%, 97.5% quantiles from each program and then checked that the absolute difference is below 0.1

My own preference is to make the comparison with WinBUGS output
a completely separate project.

You need to run a lot of iterations to get good 2.5% intervals.
And this is likely to be where we disagree with BUGS by getting
more into tails. And this is also why Gibbs shouldn't be the gold
standard. See the Betancourt and Girolami paper on hierarchical
models.

I would also strongly prefer to do synthetic code generation for
tests so that we know the answer rather than relying on comparisons
with BUGS.

But all this can be in a separate project. What we really need
are the models translated as a first priority.

> 4. The simulation of the fake data should be done somehow with the fixed_param algo from rstan as Bob points out.
>
> Could we place this stuff in some git repo?

Agreed. Ben --- we can give you push permission if you give us your
GitHub handle. For us, nothing exists that's not in GitHub.

- Bob

[Andrew Gelman]

Just a quick comment that I don’t think we need to insist on precise 2.5% intervals. That’s old-style thinking, in which the goal of a statistical analysis is to check whether the 95% interval includes zero.

In practice, the interval [-1.9, 2.1] is as good as [-2, 2] or [-2, 1.9] or whatever.

[Ben Lambert]

Hi Bob,

Have just set up a github handle for this, it's: 'benlambertStan'. I will acquaint myself with your rules for pushing and pulling etc. and will push the relevant files onto github after I have privileges.

Best,

Ben

On Tue, Jun 23, 2015 at 11:54 PM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

Just a quick comment that I don’t think we need to insist on precise 2.5% intervals.  That’s old-style thinking, in which the goal of a statistical analysis is to check whether the 95% interval includes zero.

In practice, the interval [-1.9, 2.1] is as good as [-2, 2] or [-2, 1.9] or whatever.

> On Jun 23, 2015, at 6:41 PM, Bob Carpenter <ca...@alias-i.com> wrote:
>
> You need to run a lot of iterations to get good 2.5% intervals.
> And this is likely to be where we disagree with BUGS by getting
> more into tails.  And this is also why Gibbs shouldn't be the gold
> standard.  See the Betancourt and Girolami paper on hierarchical
> models.

[Bob Carpenter]

We're much more laid back on the stan-dev/example-models repo,
which is where this should go.

The best way to get things reviewed going forward is through pull
requests. Then it's easy to leave notes in the code.

There are already BUGS models there and we want to remove
the old ones and insert the new ones.

My main goal for this is to

a) show people how to directly translate the BUGS models,

b) make sure they know we don't condone the BUGS priors,
and

c) provide an example of how we think the models should be
written.

Of these, (c) is most important for us (that is the Stan dev
team). But I understand lots of people are interested in (a)
for both speed comparisons (why you want to make the BUGS models
good implementations) and to satisfy themselves we're getting
the same answers.

After that, I think it'd be great if we could figure out a
way to unit test all these. But that's something we need to
figure out in general, not for this immediate project.

- Bob

P.S. I sent you an invite to give you permission to push to
example-models.

[Michael Betancourt]

Regarding the diagnostics:

The n_leapfrog isn’t technically what you want to be looking for, and is
mostly redundant with treedepth anyways, so I’d just look at divergences
and treedepth. I actually like showing a histogram of treedepth here, with
a line at 10, but this count is fine.

Also you have to be really careful with the language — the diagnostics
passing do not indicate convergence/mixing but rather show no signs of
non-convergence/non-mixing. This is a subtle but important point as
we can never guarantee convergence.

Additionally you don’t quote how many chains you run. I’d make this
explicit in the call and emphasize the robustness given by running
multiple chains.

> The overlplotting of orange and blue (great color choice, by the
> way) in the BUGS vs. Stan is a bit hard to read, but I don't think
> it's a big deal. What is an issue is that you probably want to be
> comparing means +/- MCMC std error, as those are the things that should
> converge to the same value if everything's working.

As Bob noted we really want to be treating Stan as the baseline and
(hoping) that we’re geometrically ergodic so the MCMC std error means
something. Assuming that this is true then what I’d recommend is

(a) rescale each parameter by subtracting the Stan mean and then
dividing by the Stan MCMC std error.

(b) The Stan parameters are all now N(0, 1) and can be represented
by a solid band from -2 to 2.

(c) Now you can plot the BUGS results as points + bars representing
the standardized means +/- 2 standardized deviations. Everything will
be on a consistent scale and deviations will be easy to spot.

You can do a similar resealing with the simulated data, only given that
the simulated data is exact you won’t need a band (just a line at 0) and
you’d want to use posterior quantiles (as you already do in Figure 1).
It would be great to see this particular plot for both Stan and BUGS.

Additionally Bob noted that performance is important. I recommend
adding a quick “Performance” section and plotting histograms of
n_eff / time for all parameters, one histogram for BUGS and one for
Stan. This can be done by overlapping transparent histograms on
one plot or just having two separate figures.

[Sebastian Weber]

Hi Bob!

To your points:

1. If you write the BUGS models with a bit of care, then they can be run with any of the three BUGS flavors. Just checkout the BUGSExamples package from Martyn Plummer. Only a handful of models need a single line of code different to run on JAGS.

2. I would like to be able to do something like this:

bugs_demo("winbugs_example_xy", engine="Stan") # this runs the "best practices" example in Stan; of course, the Stan engine will be the default
bugs_demo("winbugs_example_xy", engine="JAGS") # this runs the "best practices" example in JAGS or OpenBUGS or WinBUGS
bugs_demo("winbugs_example_xy", engine="JAGS", data="fake") # this runs the "best practices" example in JAGS or OpenBUGS or WinBUGS on fake data (which is also returned)
bugs_demo("winbugs_example_xy_vanilla", engine="JAGS") # this runs the vanilla example in JAGS or OpenBUGS or WinBUGS

3. This is really something which I do not ask for inclusion in this exercise here, but rather would like to be able to program it myself given that step 2 is done. Maybe a better wording is "with respect to the WinBUGS" result? The basic idea I am proposing is to first standardize the samples using robust estimators of the location and dispersion. This is essentially the same what Michael is suggesting and it effectively makes the 2-sample problem a 1-sample problem which is ok in this case, I think. And yes, going down to 2.5% and up to 97.5% requires a large sample size - using normal approximations one can get a handle on this.

Of course, the 2.5% steps are an arbitrary choice (one could even randomize this z-range which is compared), but then, the 95% is the magic number which is old, yes, but certainly used very often such that it is still up to date - at least for the people who are looking into this.

Best,
Sebastian

On Wednesday, June 24, 2015 at 12:42:36 AM UTC+2, Bob Carpenter wrote:
> > On Jun 23, 2015, at 3:33 PM,
> >

[Bob Carpenter]

> n_eff / time for all parameters, one histogram forBUGS and one for

> Stan. This can be done by overlapping transparent histograms on
> one plot or just having two separate figures.

Andrew objected to histogams (like Peter Li produced when comparing
the ensemble samplers) and wanted scatterplots. Then I suggested
faceting on parameter name or printing the parameters by color. The best
I could come up with was this, which shows a plot with equal axes and
lines coresponding to factors of 10 speedup off the diagonal, and is
followed by a zoom.

[Michael Betancourt]

> Andrew objected to histogams (like Peter Li produced when comparing
> the ensemble samplers) and wanted scatterplots. Then I suggested
> faceting on parameter name or printing the parameters by color. The best
> I could come up with was this, which shows a plot with equal axes and
> lines coresponding to factors of 10 speedup off the diagonal, and is
> followed by a zoom.

Well, I’m happy to argue for them as global comparisons.

I think the scatter plot is also good, although I don’t think the parameters
need to be individually identified (a ShinyStan plot where the parameter
names pop-up on rollover would be great). Also, coloring the bands
(imagine a light red band between 1x and 10x faster, dark red band
between 10x and 100x, a light blue band between 1x and 10x slower,
etc) would be awesome and allow us to drop the vertical and horizontal
lines (which aren’t super meaningful here) and the tick labels.

[Bob Carpenter]

I worry tremendously about feature creep leading to relatively
simple and useful projects never getting finished because they
won't hit the ultimate goal.

Project 1
-------------
It would be very useful to us to just have the Stan versions
of the BUGS models rewritten --- direct translation and
the model we think is better. And by direct translation, I mean
no-holds-barred best way we can code it in Stan to get
the same posterior. And I mean just the models. Nothing else.
No LaTeX, no comparisons with BUGs, etc. etc.

Project 1.5
---------------
Writing up nice doc for all these models explaining why the Stan
models are the way they are to help BUGS users program in Stan.
Note that this does not need comparison with BUGS output.

Project 2
--------------
Generating an R package that runs models in Stan, WinBUGS,
OpenBUGS, JAGS, Laplace's demon, MCMCpack, and anything else people
want to throw at it. Hopefully the dependencies can be optional for that.

Project 3
-------------
Comparison with BUGS (JAGS, OpenBUGs, ...?). There are really
two related things here.

We (the Stan devs, that is), would much prefer to do this on our
own turf writing models the way we want them for all the reasons
Michael has laid out computationally and Andrew has laid out
statistically.

Sebastian I believe wants to compare the BUGS models as written
by BUGS (or perhaps as translated by Martyn Plummer to JAGS --- I
was surprised about the flexiblity here --- I'd think the original
WinBUGS models would be a requirement).

- Bob

[Bob Carpenter]

I think we're aiming for something that can be printed
in a paper and submitted somewhere.

I think it's helpful to identify the parameter groups, but
they're only a few keystrokes to remove.

I agree that bands would be better --- that I don't know how
to do offhand. We need a color analyst for this project. I
hate having to choose colors. And in some models, we're going to
see much more than 100x speedup, so we need a general scheme
with extendible colours. Maybe one color shaded?

- Bob

[Michael Betancourt]

Naw, we’d need two. There are lots of nice divergent
color palettes out there. ColorBrewer has some nice ones,
http://colorbrewer2.org (click on 7 classes and then select
“diverging”).