[R] CODA vs. BOA discrepancy

[Martyn Plummer]

On Tue, 2005-01-25 at 03:58 -0500, gc4 at duke.edu wrote:
> Dear List:
> 
> the CODA and BOA packages for the analysis of MCMC output yield different
> results on two dignostic test of convergence: 1) Geweke's convergence
> diagnostic; 2) Heidelberger and Welch's convergence diagnostic. Does that
> imply that the CODA and BOA packages implement different ``flavors'' of
> the same test?
> 
> I paste below an example.
> 
> Geweke's test
> cbind(coda.gwk,boa.gwk)
> 
>       z-score p-value Z-Score p-value
> b[1]    1.143   0.253   0.000   1.000
> b[2]    0.470   0.638   0.056   0.956
> b[3]   -1.037   0.300  -0.388   0.698
> b[4]   -0.618   0.536  -0.085   0.933
> tau    -0.206   0.837  -0.008   0.994
> sigma   0.716   0.474   0.437   0.662
> 
> CODA -- Heidelberger and Welch
> 
>       Stationarity start     p-value
>       test         iteration
> b[1]  passed       1         0.2649
> b[2]  passed       1         0.6709
> b[3]  passed       1         0.6376
> b[4]  passed       1         0.3673
> tau   passed       1         0.1944
> sigma passed       1         0.0725
> 
>       Halfwidth Mean    Halfwidth
>       test
> b[1]  passed    -39.800 0.303994
> b[2]  passed      0.714 0.003505
> b[3]  passed      1.297 0.010317
> b[4]  passed     -0.153 0.004025
> tau   passed      0.106 0.000918
> sigma passed      3.193 0.014986
> 
> BOA -- Heidelberger and Welch
> 
>       Stationarity Test Keep Discard      C-von-M
> b[1]             passed 5000       0 2.841419e-13
> b[2]             passed 5000       0 1.512525e-03
> b[3]             passed 5000       0 1.162634e-02
> b[4]             passed 5000       0 2.589256e-03
> tau              passed 5000       0 2.742828e-04
> sigma            passed 5000       0 1.086017e-01
> 
>        Halfwidth Test        Mean    Halfwidth
> b[1]           failed -39.8000008 2.569228e+04
> b[2]           passed   0.7137679 2.777682e-02
> b[3]           passed   1.2968816 2.812990e-02
> b[4]           failed  -0.1527912 2.778146e-02
> tau            failed   0.1063554 2.772179e-02
> sigma          passed   3.1930083 2.872075e-02
> 
> Thank you very much!
> giacomo
> 

I have a long list of bug reports for CODA, with patches, from Russell
Almond. When I have folded these in I suggest you try again, and if the
problem persists, send me your data so I can look into it.

Both of these diagnostics require an estimate of the spectral density at
zero.  I am using a modification of the method suggested by Heidelburger
and Welch (1981) - using GLMs insted of smoothing and transforming the
data to fit a linear regression model - but I haven't validated it, so
there might be circumstances where it performs badly.  In particular, if
you have a short run with high autocorrelation, there is no reliable
method of estimating the spectral density at zero, so if BOA is using
another method, such as batch means, there may be big differences.

Have a look at the help page for the function spectrum0().

Martyn