updates on rhat and ess

Author: avehtari

src/cmdstan-guide/stansummary.qmd

@@ -12,22 +12,20 @@ diagnostic statistics on the sampler chains, reported in the following order:
 
 - Mean - sample mean
 - MCSE - Monte Carlo Standard Error, a measure of the amount of noise in the sample
-- StdDev - sample standard deviation - the variance around the sample mean.
-- MAD - Median Absolute Deviation  - the variance around the sample median.
+- StdDev - sample standard deviation - the standard deviation around the sample mean.
+- MAD - Median Absolute Deviation  - the median absolute deviation around the sample median.
 - Quantiles - default 5%, 50%, 95%
 - ESS_bulk
 - ESS_tail
-- R_hat - $\hat{R}$ statistic, a measure of chain equilibrium, must be within $0.05$ of  $1.0$.
+- R_hat - $\hat{R}$ statistic, a MCMC convergence diagnostic
 
 When reviewing the `stansummary` output, it is important to check the final three
-output columns first - these are the diagnostic statistics on chain convergence and
-number of independent draws in the sample.
-A $\hat{R}$ statistic of greater than $1.01$ indicates that the chain has not converged and
-therefore the sample is not drawn from the posterior, thus the estimates of the mean and
-all other summary statistics are invalid.
+output columns first - these are the diagnostic statistics on MCMC convergence and
+effective sample size.
+A $\hat{R}$ statistic of greater than $1$ indicates potential convergence problems and that the sample is not presentative of the target posterior, thus the estimates of the mean and all other summary statistics are likely to be invalid. A value $1.01$ can be used as generic threshold to decide whether more iterations or further convergence analysis is needed, but other thresholds can be used depending on the specific use case.
 
 Estimation by sampling produces an approximate value for the model parameters;
-the MCSE statistic indicates the amount of noise in the estimate.
+the MCSE statistic indicates the amount of uncertainty in the estimate.
 Therefore MCSE column is placed next to the sample mean column,
 in order to make it easy to compare this sample with others.
 
@@ -37,8 +35,9 @@ chapter of the Stan Reference Manual which describes both the theory and practic
 estimation techniques.
 
 The statistics - Mean, StdDev, MAD, and Quantiles - are computed directly from all draws across all chains.
-The diagnostic statistics - MCSE, ESS_bulk, ESS_tail, and R_hat are computed from the rank-normalized,
-folded chains according to the definitions in @Vehtari+etal:2021:Rhat.
+The diagnostic statistics - ESS_bulk, ESS_tail, and R_hat are computed from the rank-normalized,
+folded, and splitted chains according to the definitions by @Vehtari+etal:2021:Rhat.
+the MCSE statistic is computed using split chain R_hat and autocorrelations.
 The summary statistics and the algorithms used to compute them are described in sections
 [Notation for samples](https://mc-stan.org/docs/reference-manual/analysis.html#notation-for-samples-chains-and-draws)
 and

src/reference-manual/analysis.qmd

@@ -5,7 +5,7 @@ pagetitle: Posterior Analysis
 # Posterior Analysis {#analysis.chapter}
 
 Stan uses Markov chain Monte Carlo (MCMC) techniques to generate
-samples from the posterior distribution for full Bayesian inference.
+draws from the posterior distribution for full Bayesian inference.
 Markov chain Monte Carlo (MCMC) methods were developed for situations
 in which it is not straightforward to make independent draws
 @Metropolis:1953.
@@ -17,7 +17,7 @@ despite the fact that it is not actually a Markov chain.
 Stan's Laplace algorithm produces a sample from a normal approximation
 centered at the mode of a distribution in the unconstrained space.
 If the mode is a maximum a posteriori (MAP) estimate,
-the samples provide an estimate of the mean and standard deviation
+the sample provides an estimate of the mean and standard deviation
 of the posterior distribution.
 If the mode is a maximum likelihood estimate (MLE),
 the sample provides an estimate of the standard error of the likelihood.
@@ -66,15 +66,15 @@ applies.  These problems are addressed in the next two sections.
 
 Stan's posterior analysis tools compute a number of summary
 statistics, estimates, and diagnostics for Markov chain Monte Carlo
-(MCMC) samples.  Stan's estimators and diagnostics are more robust in
+(MCMC) sample.  Stan's estimators and diagnostics are more robust in
 the face of non-convergence, antithetical sampling, and long-term
 Markov chain correlations than most of the other tools available.  The
 algorithms Stan uses to achieve this are described in this chapter.
 
 
 ## Convergence
 
-By definition, a Markov chain generates samples from the target
+By definition, a Markov chain samples from the target
 distribution only after it has converged to equilibrium (i.e.,
 equilibrium is defined as being achieved when $p(\theta^{(n)})$ is the
 target density).  The following point cannot be expressed strongly
@@ -105,21 +105,23 @@ One way to monitor whether a chain has converged to the equilibrium
 distribution is to compare its behavior to other randomly initialized
 chains.  This is the motivation for the @GelmanRubin:1992 potential
 scale reduction statistic, $\hat{R}$.  The $\hat{R}$ statistic
-measures the ratio of the average variance of samples within each
-chain to the variance of the pooled samples across chains; if all
+measures the ratio of the average variance of drawss within each
+chain to the variance of the pooled draws across chains; if all
 chains are at equilibrium, these will be the same and $\hat{R}$ will
 be one.  If the chains have not converged to a common distribution,
 the $\hat{R}$ statistic will be greater than one.
 
 Gelman and Rubin's recommendation is that the independent Markov
 chains be initialized with diffuse starting values for the parameters
-and sampled until all values for $\hat{R}$ are below 1.1.  Stan
+and sampled until all values for $\hat{R}$ are below some threshold.
+@Vehtari+etal:2021:Rhat suggest in general to use a threshold $1.01$, but
+othe thresholds can be used depending on the use case. Stan
 allows users to specify initial values for parameters and it is also
 able to draw diffuse random initializations automatically satisfying
 the declared parameter constraints.
 
 The $\hat{R}$ statistic is defined for a set of $M$ Markov chains,
-$\theta_m$, each of which has $N$ samples $\theta^{(n)}_m$.  The
+$\theta_m$, each of which has $N$ draws $\theta^{(n)}_m$.  The
 *between-chain variance* estimate is
 
 $$
@@ -190,12 +192,12 @@ because the first half of each chain has not mixed with the second
 half.
 
 
-### Rank-normalization helps when there are heavy tails {-}
+### Rank normalization helps when there are heavy tails {-}
 
 Split R-hat and the effective sample size (ESS) are well defined only if
 the marginal posteriors have finite mean and variance.
-Therefore, following @Vehtari+etal:2021:Rhat, we compute the rank-normalized
-paramter values and then feed them into the formulas for split R-hat and ESS.
+Therefore, following @Vehtari+etal:2021:Rhat, we compute the rank normalized
+parameter values and then feed them into the formulas for split R-hat and ESS.
 
 Rank normalization proceeds as follows:
 
@@ -210,7 +212,31 @@ $$
 z_{(nm)} = \Phi^{-1} \left( \frac{r_{(nm)} - 3/8}{S - 1/4} \right)
 $$
 
+To further improve sensitivity to chains having different scales, 
+
+rank normalized R-hat is computed also for the 
+for the corresponding *folded*
+draws $\zeta^{(mn)}$, absolute deviations from the median,
+$$
+\label{zeta}
+\zeta^{(mn)} = \left|\theta^{(nm)}-{\rm median}(\theta)\right|.
+$$
+The rank normalized split-$\widehat{R}$ measure computed on the
+ $\zeta^{(mn)}$ values  is called \emph{folded-split}-$\widehat{R}$.
+  This measures convergence in the
+tails rather than in the bulk of the distribution. 
+
+To obtain a single conservative $\widehat{R}$ estimate, we propose
+to report the maximum of rank normalized split-$\widehat{R}$ and
+rank normalized folded-split-$\widehat{R}$ for each parameter.
+
+Bulk-ESS is defined as ESS for rank normalized split chains.  Tail-ESS
+is defined as the minimum ESS for the 5% and 95% quantiles.  See
+[Effective Sample Size section](#effective-sample-size.section) for
+details on how ESS is estimated.
+
 ### Convergence is global {-}
+
 A question that often arises is whether it is acceptable to monitor
 convergence of only a subset of the parameters or generated
 quantities.  The short answer is "no," but this is elaborated
@@ -275,7 +301,7 @@ analytically.
 So what we do in practice is hope that the finite number of draws is
 large enough for the expectations to be reasonably accurate. Removing
 warmup iterations improves the accuracy of the expectations but there
-is no guarantee that removing any finite number of samples will be
+is no guarantee that removing any finite number of draws will be
 enough.
 
 
@@ -297,7 +323,7 @@ distribution for the log posterior density (`lp__`) almost
 never looks Gaussian, instead it features long tails that can lead to
 large $\hat{R}$ even in the large $N$ limit.  Tweaks to $\hat{R}$,
 such as using quantiles in place of raw values, have the flavor of
-making the samples of interest more Gaussian and hence the $\hat{R}$
+making the sample of interest more Gaussian and hence the $\hat{R}$
 statistic more accurate.
 
 
@@ -316,7 +342,7 @@ and can apply the standard tests.
 ## Effective sample size {#effective-sample-size.section}
 
 The second technical difficulty posed by MCMC methods is that the
-samples will typically be autocorrelated (or anticorrelated) within a
+draws will typically be autocorrelated (or anticorrelated) within a
 chain.  This increases (or reduces) the uncertainty of the estimation of posterior
 quantities of interest, such as means, variances, or quantiles; see
 @Geyer:2011.
@@ -371,7 +397,7 @@ $$
    \, d\theta - \frac{\mu^2}{\sigma^2}.
 $$
 
-The effective sample size of $N$ samples generated by a process with
+The effective sample size of $N$ draws generated by a process with
 autocorrelations $\rho_t$ is defined by
 $$
 N_{\mathrm{eff}}
@@ -397,7 +423,7 @@ correlation.
 In practice, the probability function in question cannot be tractably
 integrated and thus the autocorrelation cannot be calculated, nor the
 effective sample size.  Instead, these quantities must be estimated
-from the samples themselves.  The rest of this section describes a
+from the draws themselves.  The rest of this section describes a
 autocorrelations and split-$\hat{R}$ based effective sample
 size estimator, based on multiple chains. As before, each chain
 $\theta_m$ will be assumed to be of length $N$.