MCMC Geweke Diagnostik

Ich verwende einen Metropolis-Sampler (C ++) und möchte die Konvergenzrate anhand der vorherigen Samples schätzen.

Eine einfach zu implementierende Diagnose, die ich gefunden habe, ist die Geweke-Diagnose , die die Differenz zwischen den beiden Stichprobenmitteln dividiert durch ihren geschätzten Standardfehler berechnet. Der Standardfehler wird aus der spektralen Dichte bei Null geschätzt.

$$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)}},$$

wobei , B zwei Fenster innerhalb der Markov-Kette sind. Ich habe einige Nachforschungen angestellt, was ^ S A θ ( 0 ) und ^ S B θ ( 0 ) sind , bin aber in ein Durcheinander von Literatur über spektrale Energiedichte und -leistung geraten$A$$B$$\hat{S_{\theta}^A}(0)$$\hat{S_{\theta}^B}(0)$Leistungsdichtespektrum , aber ich bin kein Experte zu diesen Themen; Ich brauche nur eine kurze Antwort: Entsprechen diese Größen der Stichprobenvarianz? Wenn nicht, wie lautet die Formel für die Berechnung?

Ein weiterer Zweifel an dieser Geweke-Diagnose ist, wie man auswählt . Die obige Literatur gesagt , dass es einige funktionelle ist θ ( X ) und sollte eine Existenz einer Spektraldichte implizieren ^ S A θ ( 0 $\theta$$\theta(X)$$\hat{S_{\theta}^A}(0)$ halber besteht der einfachste Weg jedoch darin, die Identitätsfunktion zu verwenden (verwenden Sie die Proben selbst). Ist das richtig?

Das R- Coda-Paket enthält eine Beschreibung, gibt jedoch auch nicht an, wie die Werte berechnet werden sollen.$S$

Sie können den Code für die geweke.diagFunktion im codaPaket durchsehen, um zu sehen, wie die Varianz über den Aufruf der spectrum.ar0Funktion berechnet wird .

---

Hier ist eine kurze Begründung für die Berechnung der spektralen Dichte eines AR ($p$) bei Null verarbeiten.

Die spektrale Dichte eines AR ($p$) mit der Frequenz verarbeiten $\lambda$ ist gegeben durch den Ausdruck:

$$f(\lambda) = \dfrac{\sigma^2}{(1-\sum_{j=1}^p\alpha_j\exp(-2\pi\iota j\lambda))^2}$$ wo $\alpha_j$ sind die autoregressiven Parameter.

Dieser Ausdruck vereinfacht die Berechnung der spektralen Dichte eines AR ($p$) verarbeiten um $0$:

$$f(0) = \dfrac{\sigma^2}{(1-\sum_{j=1}^p\alpha_j)^2}$$

Die Berechnung würde dann ungefähr so ​​aussehen (anstelle der üblichen Schätzer für Parameter):

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)

Überprüfen Sie die Wikipedia-Seite . Du wirst sehen$S_{xx}(\omega)$, das ist die spektrale Dichte. In Ihrem Fall sollten Sie verwenden$S_{xx}(0)$.