Introduction to Bayesian Analysis Procedures

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. If the chain has converged, early and later parts of the chain should look similar. The statistic is constructed as follows. Two subsequences of the Markov chain are taken out, with and , where . The first subsequence is called the “first window,” and the second sequence is called the “second window.” Let denote the length of the second window, and define the means of the two windows as

Let and denote consistent spectral density estimates at zero frequency (see the section Spectral Density Estimate at Zero Frequency for estimation details) for each window. If and , where the fractions and are fixed, , and the chain is stationary, then the following statistic converges to a standard normal distribution as :

This is a two-sided test, and the p-value is computed as , where X is a standard normal distribution. A p-value less than indicates a rejection of the null hypothesis that the first window of the chain and the second window of the chain come from the same distribution. If the null hypothesis is rejected, the Markov chain has not yet converged.

For one sequence of the Markov chain , the relationship between the h-lag covariance sequence of a time series and the spectral density, f, is

where i indicates that is the complex argument. Inverting this Fourier integral,

It follows that

which gives an autocorrelation adjusted estimate of the variance. In this equation, is the naive variance estimate of the sequence and is the lag h autocorrelation. Due to obvious computational difficulties, such as calculation of autocorrelation at infinity, you cannot effectively estimate by using the preceding formula. The usual route is to first obtain the periodogram of the sequence, and then estimate by smoothing the estimated periodogram. The periodogram is defined to be

The procedures use the following way to estimate from p (Heidelberger and Welch 1981). In , let and .^[2] A smooth spectral density in the domain of is obtained by fitting a gamma model with the log link function, using as response and as the only regressor. The predicted value is given by

where and are the estimates of the intercept and slope parameters, respectively.

The Heidelberger-Welch test (Heidelberger and Welch 1981, 1983) consists of two parts: a stationary 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 , set , , and . You can construct the following sequence with s coordinates on values from :

where is the rounding operator and 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, converges in distribution to a Brownian bridge (Billingsley 1986). So you can construct a test statistic by using . The statistic that these procedures use is the Cramér–von Mises statistic;^[3] that is, . As , the statistic converges in distribution to a standard Cramér–von Mises distribution. The integral is numerically approximated using Simpson’s rule.

Let , where and . If n is even, let ; otherwise, let . The Simpson’s approximation to the integral is

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

This is a one-sided test, and the p-value is computed as , where X has a standard Cramér–von Mises distribution. A p-value less than indicates rejection of the null hypothesis that the chain is stationary. 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, , is first used to construct the Cramér–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 of the chain and redo the test by using the remaining . 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 ).

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, whereas running it longer would not provide any additional information that is needed. The statistic that is calculated here is the relative half-width (RHW) of the confidence interval. The RHW for a confidence interval of level is

where is the z-score of the th percentile (for example, if ), is the variance of the chain estimated using the spectral density method (see the explanation in the section Spectral Density Estimate at Zero Frequency), n is the length, and is the estimated mean. The RHW quantifies accuracy of the level confidence interval of the mean estimate by measuring the ratio between the sample standard error of the mean and the mean itself. The test assumes that you can stop the Markov chain if the variability of the mean stabilizes with respect to the mean—in other words, if the RHW is small enough. 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.

To perform the stationarity test, you need to select an level (the default is 0.05). To perform the half-width test, you need to select another level (again, the default is 0.05) and a predetermined tolerance value (the default is 0.1). If the calculated RHW is greater than , you conclude that there are not enough data to accurately estimate the mean with confidence under a tolerance of .

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. The Raftery-Lewis test evaluates the accuracy of a quantile estimate of order . In particular, it determines the number of MCMC iterations, M, that need to be discarded (burn-ins) and the number of iterations that are needed, N, in order to achieve a desired level of precision in the estimate . Details follow. Notation and deductions here closely resemble those in Raftery and Lewis (1995).

The quantile is defined such that , where q can be an arbitrary cumulative probability, such as 0.025. This can be empirically estimated by finding the 100nqth number of the sorted . Let denote the estimand, which corresponds to an estimated probability . Because the simulated posterior distribution converges to the true distribution as the simulation sample size grows, can achieve any degree of accuracy if the simulator is run for a very long time. However, running a simulation for too long 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 of the true cumulative probability q, with probability s, such as . 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 percentage points with probability 0.95. In order to perform the diagnostic, a value must be chosen for the tolerance level of a stationary test, described in the following paragraphs. For example, .

The Raftery-Lewis diagnostic test finds the number of iterations, M, that need to be discarded (burn-ins) and the number of iterations that are needed, N, in order to achieve a desired precision. Given a predefined cumulative probability q, these procedures first find , and then they construct a binary process by setting if and 0 otherwise for all t. The sequence is itself not a Markov chain, but you can construct a subsequence of that is approximately Markovian if it is sufficiently k-thinned. When k becomes reasonably large, 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 depends on the previous j values. For example, in a second-order Markov model,

where . Given , you can construct two transition count matrices for a second-order Markov model:

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

where is the expected cell count of under the null model, the first-order Markov model, where the assumption holds. The formula for the expected cell count is

The BIC is , where is the k-thinned sample size (every kth sample, starting with the first), with the last two data points discarded because of 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 :

Because is a Markov chain, its equilibrium distribution exists and is estimated by

where and . The goal is to find an iteration number m such that after m steps, the estimated transition probability is within of equilibrium for . Let and . The estimated transition probability after step m is

which holds when

assuming .

Therefore, by time m, is sufficiently close to its equilibrium distribution, and you know that a total size of should be discarded as the burn-in.

Next, the procedures estimate N, the number of simulations that are needed in order to achieve the desired accuracy in percentile estimation. The estimate of is . For large n, is normally distributed with mean q, the true cumulative probability, and variance

is satisfied if

Therefore, .

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

If does not have that required sample size, the Raftery-Lewis test is not performed. If you still want to perform the test, increase the number of Markov chain iterations.

The ratio is sometimes referred to as the dependence factor. It measures deviation from posterior sample independence: the closer it is to 1, the less correlated the samples are. 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 are inversely proportional to , 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. Finally, the time until convergence for a parameter can differ substantially for different cumulative probabilities.

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

The sample autocovariance function of lag h of is defined by

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 can be interpreted as the equivalent number of independent draws from the posterior distribution that the Markov chain represents.

ESS is defined as

where n is the total sample size and is the autocorrelation of lag k for . The quantity is referred to as the autocorrelation time. To estimate , the Bayesian procedures first find a cutoff point k after which the autocorrelations are very close to zero, and then sum all the up to that point. The cutoff point k is such that , where is the estimated standard deviation:

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

ESS is costly to compute when the Markov chain is highly autocorrelated because it requires so many autocorrelation lags to be computed first. To remedy this problem, you can specify the cutoff point to be no larger than some maximum lag. In the case where more than the maximum number of lags would be required, this reduces computational cost but results in poor-quality estimates of ESS and . However, when the maximum lag is set high enough, this problem occurs only when the Markov chain is mixing very poorly, which is precisely the issue that ESS and are meant to diagnose.