Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

I'm running a Metropolis sampler (C++) and want to use the previous samples to estimate the convergence rate.

One easy to implement diagnostic I found is the Geweke diagnostic, which computes the difference between the two sample means divided by itsestimated standard error. The standard error is estimated from the spectral density at zero.

$$Z_n=\frac{\bar{\theta}_A-\bar{\theta}_B}{\sqrt{\frac{1}{n_A}\hat{S_{\theta}^A}(0)+\frac{1}{n_B}\hat{S_{\theta}^B}(0)}},$$

where $A$, $B$ are two windows within the Markov chain. I did some research on what are $\hat{S_{\theta}^A}(0)$ and $\hat{S_{\theta}^B}(0)$ but get into a mess of literature on energy spectral density and power spectral density. (http://en.wikipedia.org/wiki/Spectral_density) but I'm not an expert on these topics; I just need a quick answer: are these quantities the same as sample variance? If not, what's the formula for computing them?

Another doubt on this Geweke diagnostic is how to pick $\theta$? The above literature said that it is some functional $\theta(X)$ and should imply an existence of a spectral density $\hat{S_{\theta}^A}(0)$, but for convenience I guess the simplest way is to use the identity function (use samples themselves). Is this correct?

The R coda package has a description but neither does it specify how to compute the $S$ values.

share|improve this question
    
you could look in the guts of the coda function geweke.diag to see what it does ... –  Ben Bolker Aug 21 '12 at 16:54

2 Answers 2

You can look through the code for the geweke.diag function in the coda package to see how the variance is computed, via the call to the spectrum.ar0 function.


Here is a short motivation of the computation of the spectral density of an AR($p$) process at zero.

The spectral density of an AR($p$) process at frequency $\lambda$ is given by the expression: $$ f(\lambda) = \dfrac{\sigma^2}{(1-\sum_{j=1}^p\alpha_j\exp(-2\pi\iota j\lambda))^2} $$ where $\alpha_j$ are the autoregressive parameters.

This expression simplifies considerably when computing the spectral density of an AR($p$) process at $0$: $$ f(0) = \dfrac{\sigma^2}{(1-\sum_{j=1}^p\alpha_j)^2} $$

The computation then would look something like this (substituting the usual estimators for parameters):

tsAR2 = arima.sim(list(ar = c(0.01, 0.03)), n = 1000)  # simulate an AR(2) process
ar2 = ar(tsAR2, aic = TRUE)  # estimate it with AIC complexity selection

# manual estimate of spectral density at zero
sdMan = ar2$var.pred/(1-sum(ar2$ar))^2

# coda computation of spectral density at zer0
sdCoda = coda::spectrum0.ar(tsAr2)$spec

# assert equality
all.equal(sdCoda, sdMan)
share|improve this answer

Check the wikipedia page. You'll see $S_{xx}(\omega)$, which is the spectral density. In your case, you should use $S_{xx}(0)$.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.