Assessing Markov Chain Convergence :: SAS/STAT(R) 13.1 User's Guide

Assessing Markov Chain Convergence

Subsections:

Visual Analysis via Trace Plots
Statistical Diagnostic Tests

Simulation-based Bayesian inference requires using simulated draws to summarize the posterior distribution or calculate any relevant quantities of interest. You need to treat the simulation draws with care. There are usually two issues. First, you have to decide whether the Markov chain has reached its stationary, or the desired posterior, distribution. Second, you have to determine the number of iterations to keep after the Markov chain has reached stationarity. Convergence diagnostics help to resolve these issues. Note that many diagnostic tools are designed to verify a necessary but not sufficient condition for convergence. There are no conclusive tests that can tell you when the Markov chain has converged to its stationary distribution. You should proceed with caution. Also, note that you should check the convergence of all parameters, and not just those of interest, before proceeding to make any inference. With some models, certain parameters can appear to have very good convergence behavior, but that could be misleading due to the slow convergence of other parameters. If some of the parameters have bad mixing, you cannot get accurate posterior inference for parameters that appear to have good mixing. See Cowles and Carlin (1996) and Brooks and Roberts (1998) for discussions about convergence diagnostics.

Visual Analysis via Trace Plots

Trace plots of samples versus the simulation index can be very useful in assessing convergence. The trace tells you if the chain has not yet converged to its stationary distribution—that is, if it needs a longer burn-in period. A trace can also tell you whether the chain is mixing well. A chain might have reached stationarity if the distribution of points is not changing as the chain progresses. The aspects of stationarity that are most recognizable from a trace plot are a relatively constant mean and variance. A chain that mixes well traverses its posterior space rapidly, and it can jump from one remote region of the posterior to another in relatively few steps. Figure 7.1 through Figure 7.4 display some typical features that you might see in trace plots. The traces are for a parameter called $\gamma$ .

Figure 7.1: Essentially Perfect Trace for $\gamma$

Figure 7.1 displays a “perfect” trace plot. Note that the center of the chain appears to be around the value 3, with very small fluctuations. This indicates that the chain could have reached the right distribution. The chain is mixing well; it is exploring the distribution by traversing to areas where its density is very low. You can conclude that the mixing is quite good here.

Figure 7.2: Initial Samples of $\gamma$ Need to be Discarded

Figure 7.2 displays a trace plot for a chain that starts at a very remote initial value and makes its way to the targeting distribution. The first few hundred observations should be discarded. This chain appears to be mixing very well locally. It travels relatively quickly to the target distribution, reaching it in a few hundred iterations. If you have a chain that looks like this, you would want to increase the burn-in sample size. If you need to use this sample to make inferences, you would want to use only the samples toward the end of the chain.

Figure 7.3: Marginal Mixing for $\gamma$

Figure 7.3 demonstrates marginal mixing. The chain is taking only small steps and does not traverse its distribution quickly. This type of trace plot is typically associated with high autocorrelation among the samples. To obtain a few thousand independent samples, you need to run the chain for much longer.

Figure 7.4: Bad Mixing, Nonconvergence of $\gamma$

The trace plot shown in Figure 7.4 depicts a chain with serious problems. It is mixing very slowly, and it offers no evidence of convergence. You would want to try to improve the mixing of this chain. For example, you might consider reparameterizing your model on the log scale. Run the Markov chain for a long time to see where it goes. This type of chain is entirely unsuitable for making parameter inferences.

Statistical Diagnostic Tests

The Bayesian procedures include several statistical diagnostic tests that can help you assess Markov chain convergence. For a detailed description of each of the diagnostic tests, see the following subsections. Table 7.1 provides a summary of the diagnostic tests and their interpretations.

Table 7.1: Convergence Diagnostic Tests Available in the Bayesian Procedures

Name	Description	Interpretation of the Test
Gelman-Rubin	Uses parallel chains with dispersed initial values to test whether they all converge to the same target distribution. Failure could indicate the presence of a multi-mode posterior distribution (different chains converge to different local modes) or the need to run a longer chain (burn-in is yet to be completed).	One-sided test based on a variance ratio test statistic. Large $\widehat{\mbox{R}}_ c$ values indicate rejection.
Geweke	Tests whether the mean estimates have converged by comparing means from the early and latter part of the Markov chain.	Two-sided test based on a z-score statistic. Large absolute z values indicate rejection.
Heidelberger-Welch (stationarity test)	Tests whether the Markov chain is a covariance (or weakly) stationary process. Failure could indicate that a longer Markov chain is needed.	One-sided test based on a Cramer–von Mises statistic. Small p-values indicate rejection.
Heidelberger-Welch (half-width test)	Reports whether the sample size is adequate to meet the required accuracy for the mean estimate. Failure could indicate that a longer Markov chain is needed.	If a relative half-width statistic is greater than a predetermined accuracy measure, this indicates rejection.
Raftery-Lewis	Evaluates the accuracy of the estimated (desired) percentiles by reporting the number of samples needed to reach the desired accuracy of the percentiles. Failure could indicate that a longer Markov chain is needed.	If the total samples needed are fewer than the Markov chain sample, this indicates rejection.
autocorrelation	Measures dependency among Markov chain samples.	High correlations between long lags indicate poor mixing.
effective sample size	Relates to autocorrelation; measures mixing of the Markov chain.	Large discrepancy between the effective sample size and the simulation sample size indicates poor mixing.

Gelman and Rubin Diagnostics

Gelman and Rubin diagnostics (Gelman and Rubin, 1992; Brooks and Gelman, 1997) are based on analyzing multiple simulated MCMC chains by comparing the variances within each chain and the variance between chains. Large deviation between these two variances indicates nonconvergence.

Define $\{ \theta ^ t\}$ , where $t = 1,\ldots ,n$ , to be the collection of a single Markov chain output. The parameter $\theta ^ t$ is the tth sample of the Markov chain. For notational simplicity, $\theta$ is assumed to be single dimensional in this section.

Suppose you have M parallel MCMC chains that were initialized from various parts of the target distribution. Each chain is of length n (after discarding the burn-in). For each $\theta ^ t$ , the simulations are labeled as $\theta ^ t_ m, \; \mbox{where} \; t=1,\dots ,n$ and $m=1,\dots , M$ . The between-chain variance B and the within-chain variance W are calculated as

$\begin{eqnarray*} B & =& \frac{n}{M-1} \sum _{m=1}^{M}(\bar{\theta }_ m^{\cdot } - \bar{\theta }_{\cdot }^{\cdot })^2, \; \; \mbox{where} \; \; \bar{\theta }_{m}^{\cdot } = \frac{1}{n} \sum _{t=1}^{n}\theta _ m^ t, \; \; \bar{\theta }_{\cdot }^{\cdot } = \frac{1}{M} \sum _{m=1}^{M}\bar{\theta }_ m^{\cdot } \\ W & =& \frac{1}{M}\sum _{m=1}^ M s_ m^2, \; \; \mbox{where} \; \; s_ m^2 = \frac{1}{n-1}\sum _{t=1}^ n (\theta _ m^ t - \bar{\theta }_{m}^{\cdot } )^2 \end{eqnarray*}$

The posterior marginal variance, $\mr {var}(\theta |{\mb {y}})$ , is a weighted average of W and B. The estimate of the variance is

$\widehat{V} = \frac{n-1}{n}W + \frac{M+1}{nM}B$

If all M chains have reached the target distribution, this posterior variance estimate should be very close to the within-chain variance W. Therefore, you would expect to see the ratio $\widehat{V} / W$ be close to 1. The square root of this ratio is referred to as the potential scale reduction factor (PSRF). A large PSRF indicates that the between-chain variance is substantially greater than the within-chain variance, so that longer simulation is needed. If the PSRF is close to 1, you can conclude that each of the M chains has stabilized, and they are likely to have reached the target distribution.

A refined version of PSRF is calculated, as suggested by Brooks and Gelman (1997), as

$\widehat{\mbox{R}}_ c = \sqrt {\frac{\hat{d}+3}{\hat{d}+1} \cdot \frac{\widehat{V}}{W}} = \sqrt {\frac{\hat{d}+3}{\hat{d}+1} \left( \frac{n-1}{n} + \frac{M+1}{nM} \frac{B}{W} \right)}$

where

$\hat{d} = \frac{2 \widehat{V}^2}{\widehat{\mr {Var}} (\widehat{V})}$

and

$\begin{eqnarray*} \widehat{\textrm{Var}} (\widehat{V}) & =& \left( \frac{n-1}{n} \right)^2 \frac{1}{M} \widehat{\textrm{Var}}(s_ m^2) + \left( \frac{M+1}{nM} \right)^2 \frac{2}{M-1} B^2 \\ & & + ~ 2\frac{(M+1)(n-1)}{n^2M} \frac{n}{M} \left( \widehat{\textrm{cov}}(s_ m^2, (\bar{\theta }_ m^{\cdot })^2) - 2 \bar{\theta }_{\cdot }^{\cdot } \widehat{\textrm{cov}}(s_ m^2, \bar{\theta }_ m^{\cdot }) \right) \end{eqnarray*}$

All the Bayesian procedures also produce an upper $100(1-\alpha /2)\%$ confidence limit of $\widehat{\mbox{R}}_ c$ . Gelman and Rubin (1992) showed that the ratio $B/W$ in $\widehat{\mbox{R}}_ c$ has an F distribution with degrees of freedom $M-1$ and $2W^2M/\widehat{\mr {Var}}(s_ m^2)$ . Because you are concerned only if the scale is large, not small, only the upper $100(1-\alpha /2)\%$ confidence limit is reported. This is written as

$\sqrt {\left( \frac{n-1}{n} + \frac{M+1}{nM} \cdot F_{1-\alpha /2} \left(M-1, \frac{2W^2}{\widehat{\mr {Var}}(s_ m^2)/M} \right) \right) \cdot \frac{\hat{d}+3}{\hat{d}+1}}$

In the Bayesian procedures, you can specify the number of chains that you want to run. Typically three chains are sufficient. The first chain is used for posterior inference, such as mean and standard deviation; the other $M-1$ chains are used for computing the diagnostics and are discarded afterward. This test can be computationally costly, because it prolongs the simulation M-fold.

It is best to choose different initial values for all M chains. The initial values should be as dispersed from each other as possible so that the Markov chains can fully explore different parts of the distribution before they converge to the target. Similar initial values can be risky because all of the chains can get stuck in a local maximum; that is something this convergence test cannot detect. If you do not supply initial values for all the different chains, the procedures generate them for you.

Geweke Diagnostics

The Geweke test (Geweke, 1992) compares values in the early part of the Markov chain to those in the latter part of the chain in order to detect failure of convergence. The statistic is constructed as follows. Two subsequences of the Markov chain $\{ \theta ^ t\}$ are taken out, with $\{ \theta _1^ t : t = 1, \dots , n_1\}$ and $\{ \theta _2^ t: t = n_ a, \dots , n \}$ , where $1<n_1<n_ a<n$ . Let $n_2 = n - n_ a + 1$ , and define

$\bar{\theta }_1 = \frac{1}{n_1} \sum _{t=1}^{n_1} \theta ^ t \; \; \mbox{and} \; \; \bar{\theta }_2 = \frac{1}{n_2} \sum _{t=n_ a}^{n} \theta ^ t$

Let $\hat{s}_1(0)$ and $\hat{s}_2(0)$ denote consistent spectral density estimates at zero frequency (see the subsection Spectral Density Estimate at Zero Frequency for estimation details) for the two MCMC chains, respectively. If the ratios $n_1/n$ and $n_2/n$ are fixed, $(n_1+n_2)/n < 1$ , and the chain is stationary, then the following statistic converges to a standard normal distribution as $n \rightarrow \infty$ :

$Z_ n = \frac{\bar{\theta }_1 - \bar{\theta _2}}{\sqrt {\frac{\hat{s}_1(0)}{n_1} + \frac{\hat{s}_2(0)}{n_2}}}$

This is a two-sided test, and large absolute z-scores indicate rejection.

Spectral Density Estimate at Zero Frequency

For one sequence of the Markov chain $\{ \theta _ t\}$ , the relationship between the h-lag covariance sequence of a time series and the spectral density, f, is

$s_{h} = \frac{1}{2\pi } \int _{-\pi }^{\pi } \exp (\mbox{\Variable{i}}\omega h) f(\omega ) d\omega$

where i indicates that $\omega h$ is the complex argument. Inverting this Fourier integral,

$f(\omega ) =\sum _{h = -\infty }^{\infty } s_{h} \exp (-\mbox{\Variable{i}} \omega h) = s_{0} \left( 1 + 2 \sum _{h=1}^{\infty } \rho _{h} \cos ( \omega h) \right)$

It follows that

$f(0) = \sigma ^{2} \left( 1+2\sum _{h=1}^{\infty }\rho _{h} \right)$

which gives an autocorrelation adjusted estimate of the variance. In this equation, $\sigma ^{2}$ is the naive variance estimate of the sequence $\{ \theta _{t}\}$ and $\rho _{h}$ is the lag h autocorrelation. Due to obvious computational difficulties, such as calculation of autocorrelation at infinity, you cannot effectively estimate $f(0)$ by using the preceding formula. The usual route is to first obtain the periodogram $p(\omega )$ of the sequence, and then estimate $f(0)$ by smoothing the estimated periodogram. The periodogram is defined to be

$p(\omega ) = \frac{1}{n} \left[ \left( \sum _{t=1}^{n} \theta _{t} \sin (\omega t) \right)^2 + \left( \sum _{t=1}^{n} \theta _{t} \cos (\omega t)\right)^2 \right]$

The procedures use the following way to estimate $\hat{f}(0)$ from p (Heidelberger and Welch, 1981). In $p(\omega )$ , let $\omega = \omega _{k} = 2 \pi k / n$ and $k=1,\ldots ,[\frac{n}{2}]$ .[3] A smooth spectral density in the domain of $(0,\pi ]$ is obtained by fitting a gamma model with the log link function, using $p(\omega _ k)$ as response and $x_1(\omega _ k)=\sqrt {3}( 4 \omega _ k / (2 \pi ) - 1)$ as the only regressor. The predicted value $\hat{f}(0)$ is given by

$\hat{f}(0) = \exp ({\hat{\beta }_0 - \sqrt {3} \hat{\beta }_1})$

where $\hat{\beta }_0$ and $\hat{\beta }_1$ are the estimates of the intercept and slope parameters, respectively.

Heidelberger and Welch Diagnostics

The Heidelberger and Welch test (Heidelberger and Welch, 1981, 1983) consists of two parts: a stationary portion test and a half-width test. The stationarity test assesses the stationarity of a Markov chain by testing the hypothesis that the chain comes from a covariance stationary process. The half-width test checks whether the Markov chain sample size is adequate to estimate the mean values accurately.

Given $\{ \theta ^ t \}$ , set $S_0 = 0$ , $S_ n = \sum _{t=1}^ n \theta ^ t$ , and $\bar{\theta } = (1/n) \sum _{t=1}^ n \theta ^{t}$ . You can construct the following sequence with s coordinates on values from $\frac{1}{n}, \frac{2}{n}, \cdots , {1}$ :

$B_ n(s) = (S_{[ns]} - [ns] \bar{\theta }) / (n\hat{p}(0)) ^{1/2}$

where $\left[ ~ \right]$ is the rounding operator, and $\hat{p}(0)$ is an estimate of the spectral density at zero frequency that uses the second half of the sequence (see the section Spectral Density Estimate at Zero Frequency for estimation details). For large n, $B_ n$ converges in distribution to a Brownian bridge (Billingsley, 1986). So you can construct a test statistic by using $B_ n$ . The statistic used in these procedures is the Cramer–von Mises statistic [4]; that is $\int _0^1 B_ n(s)^2 ds = \mr {CVM}(B_ n)$ . As $n \rightarrow \infty$ , the statistic converges in distribution to a standard Cramer–von Mises distribution. The integral $\int _0^1 B_ n(s)^2 ds$ is numerically approximated using Simpson’s rule.

Let $y_ i = B_ n(s)^2$ , where $s = 0, \frac{1}{n}, \cdots , \frac{n-1}{n}, 1$ , and $i = ns = 0, 1, \cdots , n$ . If n is even, let $m = n/2$ ; otherwise, let $m = (n-1)/2$ . The Simpson’s approximation to the integral is

$\int _0^1 B_ n(s)^2 ds \approx \frac{1}{3n} \left[ y_0 + 4(y_1 + \cdots + y_{2m-1}) + 2(y_2 + \cdots + y_{2m-2} ) + y_{2m} \right]$

Note that Simpson’s rule requires an even number of intervals. When n is odd, $y_ n$ is set to be 0 and the value does not contribute to the approximation.

This test can be performed repeatedly on the same chain, and it helps you identify a time t when the chain has reached stationarity. The whole chain, $\{ \theta ^ t\}$ , is first used to construct the Cramer–von Mises statistic. If it passes the test, you can conclude that the entire chain is stationary. If it fails the test, you drop the initial $10\%$ of the chain and redo the test by using the remaining $90\%$ . This process is repeated until either a time t is selected or it reaches a point where there are not enough data remaining to construct a confidence interval (the cutoff proportion is set to be $50\%$ ).

The part of the chain that is deemed stationary is put through a half-width test, which reports whether the sample size is adequate to meet certain accuracy requirements for the mean estimates. Running the simulation less than this length of time would not meet the requirement, while running it longer would not provide any additional information that is needed. The statistic calculated here is the relative half-width (RHW) of the confidence interval. The RHW for a confidence interval of level $1-\alpha$ is

$\mr {RHW}=\frac{z_{(1-\alpha /2)} \cdot \left( \hat{s}_{n}/n\right) ^{1/2}}{ \hat{\theta }}$

where $z_{(1-\alpha /2)}$ is the z-score of the $100(1-\alpha /2)$ percentile (for example, $z_{(1-\alpha /2)} = 1.96$ if $\alpha = 0.05$ ), $\hat{s}_{n}$ is the variance of the chain estimated using the spectral density method (see explanation in the section Spectral Density Estimate at Zero Frequency), n is the length, and $\hat{\theta }$ is the estimated mean. The RHW quantifies accuracy of the $1-\alpha$ level confidence interval of the mean estimate by measuring the ratio between the sample standard error of the mean and the mean itself. In other words, you can stop the Markov chain if the variability of the mean stabilizes with respect to the mean. An implicit assumption is that large means are often accompanied by large variances. If this assumption is not met, then this test can produce false rejections (such as a small mean around 0 and large standard deviation) or false acceptance (such as a very large mean with relative small variance). As with any other convergence diagnostics, you might want to exercise caution in interpreting the results.

The stationarity test is one-sided; rejection occurs when the p-value is greater than $1-\alpha$ . To perform the half-width test, you need to select an $\alpha$ level (the default of which is 0.05) and a predetermined tolerance value $\epsilon$ (the default of which is 0.1). If the calculated RHW is greater than $\epsilon$ , you conclude that there are not enough data to accurately estimate the mean with $1-\alpha$ confidence under tolerance of $\epsilon$ .

Raftery and Lewis Diagnostics

If your interest lies in posterior percentiles, you want a diagnostic test that evaluates the accuracy of the estimated percentiles. The Raftery-Lewis test (Raftery and Lewis, 1992, 1995) is designed for this purpose. Notation and deductions here closely resemble those in Raftery and Lewis (1995).

Suppose you are interested in a quantity $\theta _ q$ such that $P(\theta \leq \theta _ q | {\mb {y}}) = q$ , where q can be an arbitrary cumulative probability, such as 0.025. This $\theta _ q$ can be empirically estimated by finding the $[n \cdot 100 \cdot q]$ th number of the sorted $\{ \theta ^ t\}$ . Let $\hat{\theta }_ q$ denote the estimand, which corresponds to an estimated probability $P(\theta \leq {\hat{\theta }_ q}) = \hat{P}_ q$ . Because the simulated posterior distribution converges to the true distribution as the simulation sample size grows, $\hat{\theta }_ q$ can achieve any degree of accuracy if the simulator is run for a very long time. However, running too long a simulation can be wasteful. Alternatively, you can use coverage probability to measure accuracy and stop the chain when a certain accuracy is reached.

A stopping criterion is reached when the estimated probability is within $\pm r$ of the true cumulative probability q, with probability s, such as $P(\hat{P}_ q \in (q-r, q+r)) = s$ . For example, suppose you want the coverage probability s to be 0.95 and the amount of tolerance r to be 0.005. This corresponds to requiring that the estimate of the cumulative distribution function of the 2.5th percentile be estimated to within $\pm 0.5$ percentage points with probability 0.95.

The Raftery-Lewis diagnostics test finds the number of iterations, M, that need to be discarded (burn-ins) and the number of iterations needed, N, to achieve a desired precision. Given a predefined cumulative probability q, these procedures first find $\hat{\theta }_ q$ , and then they construct a binary $0-1$ process $\{ Z_ t\}$ by setting $Z_ t = 1$ if $\theta ^ t \leq \hat{\theta }_ q$ and 0 otherwise for all t. The sequence $\{ Z_ t\}$ is itself not a Markov chain, but you can construct a subsequence of $\{ Z_ t\}$ that is approximately Markovian if it is sufficiently k-thinned. When k becomes reasonably large, $\{ Z_ t^{(k)}\}$ starts to behave like a Markov chain.

Next, the procedures find this thinning parameter k. The number k is estimated by comparing the Bayesian information criterion (BIC) between two Markov models: a first-order and a second-order Markov model. A jth-order Markov model is one in which the current value of $\{ Z_ t^{(k)}\}$ depends on the previous j values. For example, in a second-order Markov model,

$\begin{eqnarray*} \lefteqn{p \left(Z_ t^{(k)} = z_ t | Z_{t-1}^{(k)} = z_{t-1}, Z_{t-2}^{(k)} = z_{t-2}, \cdots , Z_{0}^{(k)} = z_{0} \right)} \\ & =& p \left(Z_ t^{(k)} = z_{t} | Z_{t-1}^{(k)} = z_{t-1}, Z_{t-2}^{(k)} = z_{t-2} \right) \end{eqnarray*}$

where $z_ i = \{ 0, 1 \} , i = 0, \cdots , t$ . Given $\{ Z_{t}^{(k)}\}$ , you can construct two transition count matrices for a second-order Markov model:

	$z_ t = 0$			$z_ t = 1$
	$z_{t-1} = 0$	$z_{t-1} = 1$		$z_{t-1} = 0$	$z_{t-1} = 1$
$z_{t-2} = 0$	$w_{000}$	$w_{010}$	$z_{t-2} = 0$	$w_{001}$	$w_{011}$
$z_{t-2} = 1$	$w_{100}$	$w_{110}$	$z_{t-2} = 1$	$w_{101}$	$w_{111}$

For each k, the procedures calculate the BIC that compares the two Markov models. The BIC is based on a likelihood ratio test statistic that is defined as

$G_{k}^2=2\sum _{i = 0}^{1} \sum _{j = 0}^{1}\sum _{l = 0}^{1} w_{ijl} \log \frac{ w_{ijl}}{\hat{w}_{ijl}}$

where $\hat{w}_{ijl}$ is the expected cell count of $w_{ijl}$ under the null model, the first-order Markov model, where the assumption $(Z_{t}^{(k)} \; \bot \; Z_{t-2}^{(k)} ) | Z_{t-1}^{(k)}$ holds. The formula for the expected cell count is

$\hat{w}_{ijl}=\frac{\sum _{i}w_{ijl}\cdot \sum _{l}w_{ijl}}{ \sum _{i}\sum _{l}w_{ijl}}$

The BIC is $G_{k}^2 - 2\log (n_ k -2)$ , where $n_ k$ is the k-thinned sample size (every kth sample starting with the first), with the last two data points discarded due to the construction of the second-order Markov model. The thinning parameter k is the smallest k for which the BIC is negative. When k is found, you can estimate a transition probability matrix between state 0 and state 1 for $\{ Z_ t^{(k)}\}$ :

$Q = \left( \begin{array}{cc} 1-\alpha & \alpha \\ \beta & 1-\beta \end{array} \right)$

Because $\{ Z_ t^{(k)}\}$ is a Markov chain, its equilibrium distribution exists and is estimated by

$\pi = (\pi _0, \pi _1) = \frac{(\beta , \alpha )}{\alpha + \beta }$

where $\pi _0 = P(\theta \leq \theta _ q | {\mb {y}})$ and $\pi _1 = 1 - \pi _0$ . The goal is to find an iteration number m such that after m steps, the estimated transition probability $P(Z_ m^{(k)} = i | Z_0^{(k)} = j)$ is within $\epsilon$ of equilibrium $\pi _ i$ for $i, j = 0, 1$ . Let $e_0 = (1, 0)$ and $e_1 = 1-e_0$ . The estimated transition probability after step m is

$P(Z_ m^{(k)} = i | Z_0^{(k)} = j) = e_ j \left[ \left( \begin{array}{cc} \pi _0 & \pi _1 \\ \pi _0 & \pi _1 \end{array} \right) + \frac{(1-\alpha -\beta )^ m}{\alpha +\beta } \left( \begin{array}{cc} \alpha & -\alpha \\ -\beta & \beta \end{array} \right) \right] e_ j’$

which holds when

$m = \frac{\log \left( \frac{(\alpha +\beta )\epsilon }{\max (\alpha , \beta )} \right) }{\log (1-\alpha -\beta )}$

assuming $1-\alpha -\beta > 0$ .

Therefore, by time m, $\{ Z_ t^{(k)}\}$ is sufficiently close to its equilibrium distribution, and you know that a total size of $M = mk$ should be discarded as the burn-in.

Next, the procedures estimate N, the number of simulations needed to achieve desired accuracy on percentile estimation. The estimate of $P(\theta \leq \theta _ q | {\mb {y}})$ is $\bar{Z}_ n^{(k)} = \frac{1}{n} \sum _{t=1}^ n Z_ t^{(k)}$ . For large n, $\bar{Z}_ n^{(k)}$ is normally distributed with mean q, the true cumulative probability, and variance

$\frac{1}{n} \frac{(2-\alpha -\beta ) \alpha \beta }{(\alpha + \beta )^3}$

$P(q-r \leq \bar{Z}_ n^{(k)} \leq q+r) = s$ is satisfied if

$n = \frac{(2-\alpha -\beta ) \alpha \beta }{(\alpha + \beta )^3} \left\{ \frac{\Phi ^{-1} \left( \frac{s+1}{2}\right)}{r}\right\} ^2$

Therefore, $N = nk$ .

By using similar reasoning, the procedures first calculate the minimal number of iterations needed to achieve the desired accuracy, assuming the samples are independent:

$N_{\mathit{min}} = \left\{ \Phi ^{-1} \left( \frac{s+1}{2} \right) \right\} ^2 \frac{q(1-q)}{r^2}$

If $\{ \theta ^ t\}$ does not have that required sample size, the Raftery-Lewis test is not carried out. If you still want to carry out the test, increase the number of Markov chain iterations.

The ratio $N / N_{\mathit{min}}$ is sometimes referred to as the dependence factor. It measures deviation from posterior sample independence: the closer it is to 1, the less correlated are the samples. There are a few things to keep in mind when you use this test. This diagnostic tool is specifically designed for the percentile of interest and does not provide information about convergence of the chain as a whole (Brooks and Roberts, 1999). In addition, the test can be very sensitive to small changes. Both N and $N_{\mathit{min}}$ are inversely proportional to $r^2$ , so you can expect to see large variations in these numbers with small changes to input variables, such as the desired coverage probability or the cumulative probability of interest. Last, the time until convergence for a parameter can differ substantially for different cumulative probabilities.

Autocorrelations

The sample autocorrelation of lag h for a parameter $\theta$ is defined in terms of the sample autocovariance function:

$\hat{\rho }_ h \left( \theta \right) =\frac{\hat{\gamma }_ h \left( \theta \right) }{\hat{\gamma }_0 \left( \theta \right) },\ \ |h|<n$

The sample autocovariance function of lag h of $\theta$ is defined by

$\hat{\gamma }_ h \left( \theta \right) =\frac{1}{n-h}\sum _{t=1}^{n-h}\left(\theta ^{t+h}-\bar{\theta }\right) \left( \theta ^{t} -\bar{\theta }\right) ,\ \ 0\leq h<n$

Effective Sample Size

You can use autocorrelation and trace plots to examine the mixing of a Markov chain. A closely related measure of mixing is the effective sample size (ESS) (Kass et al., 1998).

ESS is defined as follows:

$\mr {ESS} = \frac{n}{\tau } = \frac{n}{1 + 2 \sum _{k=1}^\infty \rho _ k(\theta )}$

where n is the total sample size and $\rho _ k(\theta )$ is the autocorrelation of lag k for $\theta$ . The quantity $\tau$ is referred to as the autocorrelation time. To estimate $\tau$ , the Bayesian procedures first find a cutoff point k after which the autocorrelations are very close to zero, and then sum all the $\rho _ k$ up to that point. The cutoff point k is such that $|\rho _ k| < \mbox{min} \left\{ 0.01, 2 s_ k \right\}$ , where $s_ k$ is the estimated standard deviation:

$s_ k = \sqrt {\left( \frac{1}{n} \left( 1 + 2\sum _{j=1}^{k-1} \hat{\rho }_ j^2(\theta ) \right) \right)}$

ESS and $\tau$ are inversely proportional to each other, and low ESS or high $\tau$ indicates bad mixing of the Markov chain.

^[3]This is equivalent to the fast Fourier transformation of the original time series $\theta _{t}$ .

^[4] The von Mises distribution was first introduced by von Mises (1918). The density function is $p(\theta | \mu \kappa ) \thicksim M(\mu , \kappa ) = [2 \pi I_0(\kappa )]^{-1} \exp (\kappa \cos (\theta - \mu )) \, \, ( 0 \leq \theta \leq 2\pi )$ , where the function $I_0(\kappa )$ is the modified Bessel function of the first kind and order zero, defined by $I_0(\kappa ) = (2 \pi )^{-1} \int _0^{2\pi } \exp (\kappa \cos (\theta - \mu )) d \theta$ .