MCMC diagnostics and finite-run limitations

Let $x_1,\dots ,x_S\in \mathbb{R}$ be a nonempty finite pooled collection, where $S\ge 1$. Give each observation its average rank $r_i$ among all $S$ values; equal values receive the average of the ranks they occupy. Let $\mathrm{\Phi }$ be the standard-normal CDF from Standard normal and normal laws. Define the Blom normal score

$z_i={\mathrm{\Phi }}^{-1}\text{\hspace{0.17em}⁣}\left(\frac{r_i-3\mathrm{/}8}{S+1\mathrm{/}4}\right)\mathrm{.}$

The standard-normal CDF is continuous, strictly increasing, and has limits $0$ and $1$ at the two infinities by The standard normal CDF is continuous and strictly increasing. The published quantile result The standard normal CDF has unique symmetric quantiles supplies the existence and uniqueness of a finite inverse for every probability in $(0,1)$. Since $1\le r_i\le S$, the displayed Blom fraction lies strictly between $0$ and $1$, so each score is well defined. This convention averages ties before applying the score transform.

Let $x_1,\dots ,x_S\in \mathbb{R}$ with $S\ge 1$, and let $x_{(1)}\le \cdots \le x_{(S)}$ be these values in nondecreasing order, retaining multiplicities. For $0<u\le 1$, the left empirical $u$-quantile is

$Q_S(u):=x_{(\lceil Su\rceil )}\mathrm{.}$

Thus ties are retained, and the empirical CDF is $F_S(t)=S^{-1}\mathrm{\#}\{i:x_i\le t\}$. In particular, for $0<u<1$, $Q_S(u)=\inf \{t:F_S(t)\ge u\}$. No interpolation between adjacent order statistics is performed.

Let $\lambda$ be Lebesgue measure on $\mathbb{R}$ and put $f(x)=1\mathrm{/}(\pi (1+x^2))$. Define the standard Cauchy law to be the measure ${\mu }_C=fd\lambda$ given by The measure with density $f$ relative to $\mu$. The proof below verifies that this measure has total mass one, so it is a probability law.

Let $(X_n{)}_{n\ge 0}$ be random elements of a measurable space $(E,\mathcal{E})$ on a common probability space, and let $h:E\to \mathbb{R}$ be measurable. Suppose $Y_n:=h(X_n)$ is a strictly stationary real process in the sense of Stationary process and canonical path shift, and suppose $0<{\sigma }_h^2:=\mathrm{Var}(Y_0)<\mathrm{\infty }$ using Moments, variance, and covariance on a probability space. For $k\ge 1$, put

${\rho }_h(k):=\frac{\mathrm{Cov}(Y_0,Y_k)}{{\sigma }_h^2}\mathrm{.}$

If ${\sum }_{k\ge 1}\mathrm{\mid }{\rho }_h(k)\mathrm{\mid }<\mathrm{\infty }$ and ${\tau }_h:=1+2{\sum }_{k\ge 1}{\rho }_h(k)>0$, define the theoretical functional-specific effective sample size for an integer $n\ge 1$ draws by ${\mathrm{ESS}}_h(n):=n\mathrm{/}{\tau }_h$. Here $Y_n$ is the scalar functional whose Monte Carlo average is being assessed. When a central limit theorem for that average has asymptotic variance ${\sigma }_h^2{\tau }_h$, this definition is the iid sample size with the same asymptotic variance. It does not assert that such a central limit theorem holds. If the covariance series is not absolutely summable or ${\tau }_h\le 0$, this definition does not assign a finite positive effective sample size.

Inputs: $M\ge 2$ scalar chains $x_{m,1},\dots ,x_{m,L}$ of equal length $L\ge 4$.

1. Let $N=\lfloor L\mathrm{/}2\rfloor$. Replace each original chain by two chains: its first $N$ observations and its last $N$ observations. If $L$ is odd, discard the middle observation. There are now $M^{\mathrm{\prime }}=2M$ chains of length $N\ge 2$.
2. For each split chain $j$, compute its mean ${\bar{x}}_j$ and sample variance $s_j^2=(N-1{)}^{-1}{\sum }_{i=1}^N(x_{j,i}-{\bar{x}}_j{)}^2$. Put ${\bar{x}}_{\cdot }=M^{\mathrm{\prime }-1}{\sum }_j{\bar{x}}_j$, $W=M^{\mathrm{\prime }-1}{\sum }_j{s}_j^2$, and $B=N(M^{\mathrm{\prime }}-1{)}^{-1}{\sum }_j({\bar{x}}_j-{\bar{x}}_{\cdot }{)}^2$.
3. If $W=0$, return not estimable. Otherwise return ${\widehat{R}}_{\mathrm{s}\mathrm{p}\mathrm{l}\mathrm{i}\mathrm{t}}=\sqrt{\frac{{\widehat{V}}^{+}}{W}},{\widehat{V}}^{+}=\frac{N-1}{N}W+\frac{B}{N}\mathrm{.}$

This finite-sample formula is not clamped below at $1$; identical split-chain means can therefore give a value slightly below $1$. Its interpretation as a convergence diagnostic requires more than the arithmetic procedure itself.

Inputs: $M\ge 2$ scalar chains $y_{m,1},\dots ,y_{m,L}$ of equal length $L\ge 4$. This estimator targets the functional-specific quantity in Functional-specific effective sample size from finite arrays; it does not test stationarity or a central limit theorem.

1. Split every input chain into its first and last $N=\lfloor L\mathrm{/}2\rfloor$ observations. If $L$ is odd, discard its middle observation. Denote the resulting $M^{\mathrm{\prime }}=2M$ chains by $z_{j,1},\dots ,z_{j,N}$ and set $S^{\mathrm{\prime }}=M^{\mathrm{\prime }}N$.
2. Compute their means ${\bar{z}}_j$, sample variances $s_j^2=(N-1{)}^{-1}{\sum }_{i=1}^N(z_{j,i}-{\bar{z}}_j{)}^2$, their average $W=M^{\mathrm{\prime }-1}{\sum }_j{s}_j^2$, grand mean ${\bar{z}}_{\cdot }=M^{\mathrm{\prime }-1}{\sum }_j{\bar{z}}_j$, between-chain term $B=N(M^{\mathrm{\prime }}-1{)}^{-1}{\sum }_j({\bar{z}}_j-{\bar{z}}_{\cdot }{)}^2$, and ${\widehat{V}}^{+}=(N-1)W\mathrm{/}N+B\mathrm{/}N$. If any $s_j^2=0$ or ${\widehat{V}}^{+}=0$, return not estimable.
3. For $0\le t<N$, compute the biased autocovariances ${\widehat{\gamma }}_{t,j}=N^{-1}{\sum }_{i=1}^{N-t}(z_{j,i}-{\bar{z}}_j)(z_{j,i+t}-{\bar{z}}_j)$, the within-chain autocorrelations ${\widehat{r}}_{t,j}={\widehat{\gamma }}_{t,j}\mathrm{/}{\widehat{\gamma }}_{0,j}$, and the between/within adjusted combined autocorrelation ${\widehat{\rho }}_t=1-\frac{W-M^{\mathrm{\prime }-1}{\sum }_j{s}_j^2{\widehat{r}}_{t,j}}{{\widehat{V}}^{+}}\mathrm{.}$
4. For $0\le k\le \lfloor (N-2)\mathrm{/}2\rfloor$, form $P_k={\widehat{\rho }}_{2k}+{\widehat{\rho }}_{2k+1}$. If $P_0\le 0$, return not estimable. Let $K$ be the largest $k$ for which $P_1,\dots ,P_k$ are all positive; $K=0$ if $P_1\le 0$ or no $P_1$ is available. Put ${\tilde{P}}_k={\min }_{0\le j\le k}{P}_j$ and ${\widehat{\tau }}_{\mathrm{o}\mathrm{d}\mathrm{d}}=-1+2{\sum }_{k=0}^K{\tilde{P}}_k$. If the next even lag $2K+2<N$ is available, put ${\widehat{\tau }}_{\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{n}}={\widehat{\tau }}_{\mathrm{o}\mathrm{d}\mathrm{d}}+2{\widehat{\rho }}_{2K+2}$; otherwise put ${\widehat{\tau }}_{\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{n}}={\widehat{\tau }}_{\mathrm{o}\mathrm{d}\mathrm{d}}$. Set $\widehat{\tau }=({\widehat{\tau }}_{\mathrm{o}\mathrm{d}\mathrm{d}}+{\widehat{\tau }}_{\mathrm{e}\mathrm{v}\mathrm{e}\mathrm{n}})\mathrm{/}2$.
5. If $\widehat{\tau }\le 0$, return not estimable. Otherwise return $\stackrel{}{\mathrm{ESS}}=\min \{\frac{S^{\mathrm{\prime }}}{\widehat{\tau }},S^{\mathrm{\prime }}{\log }_{10}{S}^{\mathrm{\prime }}\}\mathrm{.}$ The second term is the paper's finite software cap for very short, strongly antithetic chains.

Negative serial correlation can make the estimate exceed the $S^{\mathrm{\prime }}$ retained split draws, subject to that software cap. The output is an estimator, not a guarantee about the target law.

Inputs: $M\ge 2$ equal-length scalar chains with at least four observations each. Let $S$ be the total number of pooled draws.

1. Apply Blom rank-normal scores for pooled draws to all pooled draws, split the transformed chains as in Split R-hat from multiple chains, and calculate the split statistic $R_z$ on the transformed values.
2. Compute the pooled median $m$ by sorting all $S$ draws and averaging the two central order statistics when $S$ is even (using the central order statistic when $S$ is odd). Fold each original draw to $\mathrm{\mid }{x}_{m,n}-m\mathrm{\mid }$, apply the same pooled rank-normal transformation to these folded values, split each folded chain, and calculate the split statistic $R_f$.
3. If either split calculation is not estimable, return not estimable. Otherwise return $\widehat{R}=\max \{R_z,R_f\}$.

The rank-normalized component addresses location and trend differences in a finite output array; folding the absolute deviations makes the reported maximum sensitive to scale differences as well. It is a diagnostic procedure, not a proof of stationarity or convergence.

Inputs: $M\ge 2$ equal-length scalar chains with $N\ge 2$ observations each. Output: A finite bulk effective-sample-size estimate or not estimable.

1. If $N<4$, return not estimable. This procedure calls Split multi-chain autocorrelation effective sample size, whose input chains must have length at least $4$.
2. Pool all draws and apply Blom rank-normal scores for pooled draws. Apply Split multi-chain autocorrelation effective sample size to the resulting $M$ transformed chains and return its output.

Rank normalization gives finite transformed values for every finite real-valued draw array, including arrays whose target distribution might have infinite mean or variance. It is an efficiency diagnostic for the transformed values and is not, by itself, a Monte Carlo standard error for the original posterior mean.

Inputs: $M\ge 2$ equal-length chains of scalar draws with $N\ge 2$ observations each. Output: The minimum of the two tail ESS estimates or not estimable.

1. If $N<4$, return not estimable. This procedure calls Split multi-chain autocorrelation effective sample size, whose input chains must have length at least $4$.
2. Pool the draws and calculate $q_{.05}=Q_S(.05)$ and $q_{.95}=Q_S(.95)$ using Left empirical quantile of finite draws. For each $p\in \{.05,.95\}$, form the bounded indicator chains $I_{m,n}^{(p)}=\mathbf{1}\{x_{m,n}\le q_p\}$, pool their values, apply Blom rank-normal scores for pooled draws, and apply Split multi-chain autocorrelation effective sample size to the transformed indicator chains.
3. If either indicator ESS is not estimable, return not estimable. Otherwise return the minimum of the two values.

This is a tail-focused efficiency diagnostic; it does not establish that a quantile estimate has a specified error bound.

Inputs: $M\ge 2$ equal-length scalar chains with $N\ge 2$ observations each, a quantile level $\alpha \in (0,1)$, and the fixed levels $u_{-}=.16$ and $u_{+}=.84$ used for a one-standard-error approximation. Output: The approximate error interval and MCSE below, or not estimable.

1. If $N<4$, return not estimable. This procedure calls Split multi-chain autocorrelation effective sample size, whose input chains must have length at least $4$.
2. Pool the draws and compute ${\widehat{q}}_{\alpha }=Q_S(\alpha )$ by Left empirical quantile of finite draws. Form the indicator chains $I_{m,n}=\mathbf{1}\{x_{m,n}\le {\widehat{q}}_{\alpha }\}$, Blom-normalize their pooled values by Blom rank-normal scores for pooled draws, and compute $E_{\alpha }$ with Split multi-chain autocorrelation effective sample size. If the indicator chains are constant or $E_{\alpha }$ is not estimable, return not estimable.
3. Set $a=E_{\alpha }\alpha +1$ and $b=E_{\alpha }(1-\alpha )+1$. For $u\in (0,1)$, define the lower Beta quantile $Q_{a,b}(u)=\inf \{t\in [0,1]:F_{a,b}(t)\ge u\}$, where $F_{a,b}(t)=\mathrm{Beta}(a,b)((0,t])$ and Beta distribution supplies the probability measure. Let $v_{-}=Q_{a,b}(u_{-})$ and $v_{+}=Q_{a,b}(u_{+})$. For integers $n\ge 1$, the sets $E_n=(0,1\mathrm{/}n]$ in $(0,1)$ decrease to the empty set. By The real numbers, $\mathbb{R}$ is the Cauchy-sequence real field, and The Cauchy-sequence reals have the least-upper-bound property says it is a complete ordered field. Thus For every $\epsilon >0$ in a complete ordered field there is a natural $n\ge 1$ with $1\mathrm{/}n<\epsilon$ gives, for every $p\in (0,1)$, an $n$ with $1\mathrm{/}n<p$. The sets have finite Beta measure because the Beta law is a probability measure. Thus Continuity from above when one set has finite measure gives $F_{a,b}(1\mathrm{/}n)=\mathrm{Beta}(a,b)(E_n)\to 0$. Since $F_{a,b}$ is nondecreasing by Measures are monotone and $F_{a,b}(1)=1$, for each $u\in \{u_{-},u_{+}\}\subset (0,1)$ there is an $n$ with $F_{a,b}(1\mathrm{/}n)<u$; no $t\le 1\mathrm{/}n$ then belongs to the defining level set, while $t=1$ does. Hence $0<1\mathrm{/}n\le Q_{a,b}(u)\le 1$, so $v_{-},v_{+}\in (0,1]$ and both calls to $Q_S$ lie in its domain. Also $u_{-}<u_{+}$ makes the upper level set a subset of the lower one, so $v_{-}\le v_{+}$. The order-statistic definition makes $Q_S$ nondecreasing, hence $Q_S(v_{-})\le Q_S(v_{+})$.
4. Return the approximate Monte Carlo error interval $[Q_S(v_{-}),Q_S(v_{+})]$ and the estimate $\stackrel{}{\mathrm{MCSE}}({\widehat{q}}_{\alpha })=\frac{Q_S(v_{+})-Q_S(v_{-})}{2}\mathrm{.}$

All Beta shape parameters are positive because $E_{\alpha }>0$ and $0<\alpha <1$. This is a finite-run approximation, not a theorem that the interval has a prescribed coverage probability. Its statistical interpretation requires conditions for the relevant quantile and indicator-chain limit behavior; none are asserted here.