Aki writes:
tl;dr
In the case of well mixing chains, the functions ess_basic and mcse_quantile in the posterior package are the most accurate methods. The ess_basic package computes sum of autocorrelations using the truncation rule of Geyer (1992), and mcse_quantile uses in addition an inverse approach for quantile Monte Carlo standard error. In the case of well separated modes, the posterior package also provides appropriately small effective sample size indicating that the chains are not mixing.
What should effective sample size estimates look like when we have chains that are not mixing? Here’s what Aki says:
If the modes are well separated so that between-chain variance is much bigger than within-chain variances, each chain is mostly indicating the loaction of the mode and not worth much more than one observation. The posterior package combines Geyer truncated autocorrelation estimates with multi-chain convergence diagnostic R-hat. If the chains are not mixing, then the draws within each chain are considered to be more correlated within the chain and the multi-chain ESS estimate gets smaller. If the modes are well separated, multi-chain ESS is close to the number of chains. Arguably even this is an overestimate as increasing the number of non-mixing chains does not provide reliable estimate of the relative masses of the masses.
This is a big deal. If you estimate effective sample size separately for each chain and then average the chains, then you’ll be much too optimistic when your chains aren’t mixing. That’s why we combine between and within-chain information as described in BDA but which is done more effectively in our recent paper, Rank-normalization, folding, and localization: An improved R-hat for assessing convergence of MCMC.
Finally, here’s the background from Aki:
The MCMC effective sample size (ESS) and Monte Carlo standard error (MCSE) estimated for one chain includes estimation of the correlation between the iterations, for example, using autocorrelation time (or spectral density at frequency zero). As a finite number of MCMC draws are available the autocorrelation estimates for bigger lags are noisier. In our paper, we wrote about autocorrelation estimators:
In our experiments, Geyer’s (1992) truncation had superior stability compared to flat-top (Doss et al., 2014) and slug-sail (Vats and Knudson, 2018) lag window approaches.
This notebook includes additional experiments . . .
Again, here’s the link to Aki’s notebook.
Given the large discrepancy between the estimated and actual ESS in the experiments, maybe the ESS itself (and more generally other diagnostics such as R hat) should come with some interval estimate.
Yuling:
For the original R-hat paper I came up with an interval estimate, but in practice I didn’t feel we were getting much out of it so I phased it out. One thing is that neither R-hat nor ESS is the goal in itself.
> Given the large discrepancy between the estimated and actual ESS in the experiments
That is an interesting conclusion from the notebook, as I thought the best estimates are very accurate except for very short and highly autocorrelated chains, and even for the less good methods I did mention that the discrepancies probably don’t have much effect in actual practical workflow.
> maybe the ESS itself (and more generally other diagnostics such as R hat) should come with some interval estimate.
As ESS is used for MCSE, this would be the same as adding error bars on error bars https://xkcd.com/2110/