raftery.diag {coda}R Documentation

Raftery and Lewis's diagnostic

Description

raftery.diag is a run length control diagnostic based on a criterion of accuracy of estimation of the quantile q. It is intended for use on a short pilot run of a Markov chain. The number of iterations required to estimate the quantile qq to within an accuracy of +/- rr with probability pp is calculated. Separate calculations are performed for each variable within each chain.

If the number of iterations in data is too small, an error message is printed indicating the minimum length of pilot run. The minimum length is the required sample size for a chain with no correlation between consecutive samples. Positive autocorrelation will increase the required sample size above this minimum value. An estimate I (the ‘dependence factor’) of the extent to which autocorrelation inflates the required sample size is also provided. Values of I larger than 5 indicate strong autocorrelation which may be due to a poor choice of starting value, high posterior correlations or ‘stickiness’ of the MCMC algorithm.

The number of ‘burn in’ iterations to be discarded at the beginning of the chain is also calculated.

Usage

raftery.diag(data, q=0.025, r=0.005, s=0.95, converge.eps=0.001)

Arguments

data

an mcmc object

q

the quantile to be estimated.

r

the desired margin of error of the estimate.

s

the probability of obtaining an estimate in the interval (q-r,q+r).

converge.eps

Precision required for estimate of time to convergence.

Value

A list with class raftery.diag. A print method is available for objects of this class. the contents of the list are

tspar

The time series parameters of data

params

A vector containing the parameters r, s and q

Niters

The number of iterations in data

resmatrix

A 3-d array containing the results: MM the length of "burn in", NN the required sample size, NminNmin the minimum sample size based on zero autocorrelation and I=(M+N)/NminI = (M+N)/Nmin the "dependence factor"

Theory

The estimated sample size for variable U is based on the process Zt=d(Ut<=u)Z_t = d(U_t <= u) where dd is the indicator function and u is the qth quantile of U. The process ZtZ_t is derived from the Markov chain data by marginalization and truncation, but is not itself a Markov chain. However, ZtZ_t may behave as a Markov chain if it is sufficiently thinned out. raftery.diag calculates the smallest value of thinning interval kk which makes the thinned chain ZtkZ^k_t behave as a Markov chain. The required sample size is calculated from this thinned sequence. Since some data is ‘thrown away’ the sample size estimates are conservative.

The criterion for the number of ‘burn in’ iterations mm to be discarded, is that the conditional distribution of ZmkZ^k_m given Z0Z_0 should be within converge.eps of the equilibrium distribution of the chain ZtkZ^k_t.

Note

raftery.diag is based on the FORTRAN program ‘gibbsit’ written by Steven Lewis, and available from the Statlib archive.

References

Raftery, A.E. and Lewis, S.M. (1992). One long run with diagnostics: Implementation strategies for Markov chain Monte Carlo. Statistical Science, 7, 493-497.

Raftery, A.E. and Lewis, S.M. (1995). The number of iterations, convergence diagnostics and generic Metropolis algorithms. In Practical Markov Chain Monte Carlo (W.R. Gilks, D.J. Spiegelhalter and S. Richardson, eds.). London, U.K.: Chapman and Hall.


[Package coda version 0.19-4.1 Index]