Question on Bayesian Cognitive Modeling code

[Shravan Vasishth]

Hi all,

I was playing with the memory retention model published on the github repository stan-dev, for the Lee and Wagenmakers book:

https://github.com/stan-dev/example-models/tree/master/Bayesian_Cognitive_Modeling/CaseStudies/MemoryRetention

I ran the model Retention_2_Stan.R with a single chain (as the authors gave it), and then when I use print to look at the output, R-hat is printed out as 1. How can R-hat be computed for a single chain? 

thanks,

Shravan

[Ben Goodrich]

For Stan, each chain is split in half when the R-hat statistic is calculated. This allows it to be calculated (technically) when there is only one chain to split and in general is better at catching chains that drift in similar ways. The Stan manual mentions this in the chapter on MCMC. Still, I wouldn't run just one chain.

Ben

[Avraham Adler]

It's also in the text and footnote in BDA3, where the splitting is mentioned. I've e-mailed Dr. Plummer asking if he implements the change in coda. FWIW, I think it's still based on the original Gelman and Rubin paper.

Avi

[Bob Carpenter]

On Sep 13, 2014, at 8:31 PM, Avraham Adler <avraha...@gmail.com> wrote:

> It's also in the text and footnote in BDA3, where the splitting is mentioned. I've e-mailed Dr. Plummer asking if he implements the change in coda. FWIW, I think it's still based on the original Gelman and Rubin paper.

Martyn said he was going to update Coda to use the split R-hat metric.

And yes, it's still based on the original R-hat --- the only difference
is that each chain's split in two to diagnose a bunch of chains all slowly
drifting in the same direction. Or to measure one long chain against
itself, though as Ben said, it's not what we'd recommend except during model
debugging. (The case of models like mixtures that don't identify parameters
is an argument for just running one long chain.)

- Bob

>
> Avi
>
> On Saturday, September 13, 2014 9:35:03 AM UTC-4, Ben Goodrich wrote:
> For Stan, each chain is split in half when the R-hat statistic is calculated. This allows it to be calculated (technically) when there is only one chain to split and in general is better at catching chains that drift in similar ways. The Stan manual mentions this in the chapter on MCMC. Still, I wouldn't run just one chain.
>
> Ben
>
> On Saturday, September 13, 2014 7:37:09 AM UTC-4, Shravan Vasishth wrote:
> Hi all,
>
> I was playing with the memory retention model published on the github repository stan-dev, for the Lee and Wagenmakers book:
>
> https://github.com/stan-dev/example-models/tree/master/Bayesian_Cognitive_Modeling/CaseStudies/MemoryRetention
>
> I ran the model Retention_2_Stan.R with a single chain (as the authors gave it), and then when I use print to look at the output, R-hat is printed out as 1. How can R-hat be computed for a single chain?
>
> thanks,
>
> Shravan
>

> --
> 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.
> For more options, visit https://groups.google.com/d/optout.

[Andrew Gelman]

No, it’s better to run 2 or 4 long chains, not just 1.

[Bob Carpenter]

We want to run multiple chains to diagnose convergence better.
Ironically, evaluating R-hat on non-identified parameters like
mixture components is what's difficult!

For most posterior inferences, the unidentified components are
going to be marginalized out anyway. Should we just not bother
monitoring convergence on the mixture component parameters?

I know this brings up how we improve global convergence monitoring.
Would that somehow help with mixture components?

In other cases than mixtures, there really shouldn't be unidentified parameters
in a proper model --- we want to identify so that we can adapt and mix
well without devolving into long random walks in the unconstrained space.

The mixture case is the one we don't seem to be
able to parameterize away due to the arguments you made about
correlation of parameters and their interaction in things like ordering.
For instance, you might identify two mixture components with mu1 and
mu2 having posterior 90% intervals (-1,1) and (0,2) or something like
that. But if you use an ordered vector, you wind up with mu1 < mu2 in
all iterations.

- Bob

[Andrew Gelman]

I think we should monitor convergence on everything because if we have lack of convergnce on some intermediate parameters, this is still something we’d like to know.
A

[Bob Carpenter]

But that's the problem --- how do we monitor convergence of
the parameters involved in the mixture components?

p(y | lambda, mu1, sigma1, mu2, sigma2)

= lambda * normal(y|mu1,sigma1)

+ (1 - lambda) * normal(y|mu2,sigma2)

(mu1,sigma1) and (mu2,sigma2) are going to be swappable. And if
there's a narrow enough trough between the modes, you might even
swap within a single chain.

- Bob

[Andrew Gelman]

Yes, this comes up sometimes. I think the parameters need to be identified in the prior. Simplest would be just to constrain mu1 < mu2 but there are other approaches.

[Bob Carpenter]

I'm so confused.

I thought you (Andrew) were the one who convinced me we couldn't
just constrain mu1 < mu2 because even if
the posterior mean of mu1 < posterior mean of mu2, it doesn't mean every
iteration would have mu1 < mu2, and putting the constraint messes
up the posterior by insisting every smaple in the scatter plot is
above the diagonal (or below depending on how you draw it).

- Bob

[Andrew Gelman]

It depends on the application. If we’re in a purely predictive setting, constraining mu1 < mu2 should work in this setting. If the parameters are supposed to have a deeper meaning, then this constraint could cause some weird artiifacts in the interpretation of the parameters.

[Ben Goodrich]

On Sunday, September 14, 2014 6:15:21 PM UTC-4, Bob Carpenter wrote:

On Sep 13, 2014, at 8:31 PM, Avraham Adler <avraha...@gmail.com> wrote:

> It's also in the text and footnote in BDA3, where the splitting is mentioned. I've e-mailed Dr. Plummer asking if he implements the change in coda. FWIW, I think it's still based on the original Gelman and Rubin paper.

Martyn said he was going to update Coda to use the split R-hat metric.

And yes, it's still based on the original R-hat --- the only difference
is that each chain's split in two to diagnose a bunch of chains all slowly
drifting in the same direction.  

Not quite the only difference: coda also does a degrees of freedom correction that was advocated in the original Gelman and Rubin paper on Rhat that Andrew has stopped advocating but Stan should probably be doing.

Ben

[Avraham Adler]

Off topic, sorry, but out of curiosity, why isn't the DoF correction in BDA3?

Avi

[Michael Betancourt]

To be completely clear — chains properly sampling the entirety of a multimodal
distribution will yield an Rhat of 1. The problem here is when we have samplers
that cannot explore modes efficiently.

If the modes are not meaningful because of, for example, label switching degeneracies
then the easy answer is to just find a transformed parameter that doesn’t depend
on the labels and compute Rhat from that. But if the modes are meaningful then
we would want Rhat to diagnose this as a problem because a chain that can’t
sample from both modes at once will have biased inferences.

Breaking degeneracies through modeling may work, but it will be delicate. In the
mixture of Gaussians example ordering will work only if | mu1 - mu2 | >> 0 so that
the overlap Bob mentions is inconsequential.

[Ben Goodrich]

On Monday, September 15, 2014 2:45:39 AM UTC-4, Avraham Adler wrote:

Off topic, sorry, but out of curiosity, why isn't the DoF correction in BDA3?

As Andrew would say (and actually has said), he got tired of explaining it.

[Avraham Adler]

LOL.

Isn't the purpose of a textbook "EORM" (explain once, read {by students} many)?

Practically speaking, is the following accurate about the DoF correction?

1. If we use coda (and packages that use it like runjags) we should be fine
2. It will be implemented in Stan going forward
3. If we need to write our own, it would be better to implement it based on the original paper, but just using what is described in BDA3 will not cause the four horsemen of the Bayesian apocalypse to manifest

Avi

[Ben Goodrich]

On Monday, September 15, 2014 1:48:31 PM UTC-4, Avraham Adler wrote:

LOL.

Isn't the purpose of a textbook "EORM" (explain once, read {by students} many)?

Well, there are now 3 editions of BDA, so ...
 

Practically speaking, is the following accurate about the DoF correction?

Not completely

• If we use coda (and packages that use it like runjags) we should be fine

Coda does not use splitting currently, so you could be vulnerable to the problem of chains drifting in harmony.

• It will be implemented in Stan going forward

Maybe

• If we need to write our own, it would be better to implement it based on the original paper, but just using what is described in BDA3 will not cause the four horsemen of the Bayesian apocalypse to manifest

Splitting is better than not splitting. The degrees of freedom correction is better than the uncorrected estimate, but it tends to only be noticeable for heavy-tailed margins (in which case the Rhat on the original variable may not be appropriate anyway). Of the two, splitting is probably a bigger deal, particularly if you only do 1 or 2 chains for some reason.

Ben