Hi all,
I’m using the rstan package and noticed that it outputs rhat even if I run only one chain. To my understanding, calculating R-hat requires at least two chains. I wonder how r hat is being calculated in rstan under this condition.
Thank you!
I’m pretty sure that Stan uses split-R-hat, so if there are m chains, you’ll get 2m pieces which can be used for R-hat. I will always run multiple chains, but if you have just 1 chain, it is still split in two.
3 Likes
Thanks for the information! If I run 1 chain with 2000 iterations, does that mean the 2000 iterations will be split into two sets of 1000 iterations each, or will I get 2000 samples in total?
If you run 1 chain with 1000 warmup and 1000 saved iterations, then the 1000 saved iterations will be saved into two sets when computing split R-hat. But they are put back together when the simulations are sent out of Stan. In general, if you run m chains for n saved iterations each, Stan will compute split R-hat by splitting them into 2m chains, each of length n/2, and the Stan will put them back and return m chains each of length n. I guess that Stan will also return the warmup iterations if you ask it, but usually I don’t do anything with them.
Yes, Andrew’s right—Stan uses split-\widehat{R}.
You can get the precise definition we use in both the Stan Reference Manual and in Gelman et al.'s Bayesian Data Analysis (free pdf on the book’s home page). Here’s the relevant section of the Reference Manual:
https://mc-stan.org/docs/reference-manual/analysis.html
The more elaborate version we use now in RStan and ArviZ is based on this paper:
This will be coming soon to CmdStan and hence CmdStanPy and CmdStanR.
It’s actually already in CmdStanR because we use the posterior package to compute diagnostics. So fit$summary() will give you the new Rhat.
To be more precise, Bayesian Data Analysis does not describe the version which has been used for years in Stan. The paper Rank-Normalization, Folding, and Localization: An Improved Rˆ for Assessing Convergence of MCMC (with Discussion) does describe in addition to the rank-normalized version also the Rhat version which has been in used for years in Stan and mentions the differences to BDA3 version which is not really used anywhere.
1 Like
I guess that depends on what you call “Stan”. It’s still being used in CmdStan and CmdStanPy, and I only use CmdStanPy these days. I have no idea how far CmdStanR has drifted from mirroring CmdStanPy. [Edit: Oops—just saw @Jonah’s post that says they’re using the convergence monitoring from posterior in cmdstanr.]
Of course, you can drop in a different posterior analysis package like ArviZ, but I really can’t stand ArviZ (specifically, its dependencies on analysis packages and its god-object design [the latter it shares with all the stats packages in R like lm/glm]).
@avehtari took the expedient of updating the convergence monitoring for some of our interfaces, but not the core C++ packages. As a result, some of our interfaces remain out of sync. @mitzimorris is now taking the time to build performant implementations in C++ and will be wrapping them in CmdStan so they’ll be available to CmdStanPy. My guess is that CmdStanR will continue to forge its own path and ArviZ will continue to use its own implementations. Nobody, including me, is prioritizing keeping our interfaces in sync.
The improved R-hat diagnostic was added by aleksgorica (Aleks) · GitHub in Improved rhat diagnostic by aleksgorica · Pull Request #3266 · stan-dev/stan · GitHub and split rank-normalized ESS was added last week Feature/3299 improved ess rhat by mitzimorris · Pull Request #3312 · stan-dev/stan · GitHub. once the performant refactoring that Bob mentioned has been reviewed - Feature/3299 chainset by mitzimorris · Pull Request #3313 · stan-dev/stan · GitHub - CmdStan’s bin/stansummary will use the same algorithms as posterior and it should be available in the next release.
2 Likes
Thankfully the additions that @mitzimorris mentioned should bring the interfaces more into alignment on diagnostics. It’s not so much that we wanted to forge our own path with CmdStanR (although you’re right that we did), it’s just that we had the improved diagnostics available in R so it seemed like we should go ahead and use them. It wasn’t clear at the time when they would end up making it into CmdStan. It’s great that we’ll have them in stansummary going forward!
1 Like
Sorry, I just realize my mistake was not what I call “Stan”, but thinking also about the effective sample size (ESS). BDA3 Section 11.5 describes a version of ESS (which includes the Rhat computation) which is not used in CmdStan.
My bad—I thought the implementation in CmdStan was the one from BDA3!