MCMC Geweke Diagnostik

14

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.

Zn=θ¯Aθ¯B1nASθA^(0)+1nBSθ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 geratenABSθA^(0)Sθ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 )θθ(X)Sθ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

Yang
quelle
Sie könnten in die Eingeweide der codaFunktion schauen, um geweke.diagzu sehen, was es tut ...
Ben Bolker

Antworten:

4

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 λ ist gegeben durch den Ausdruck:

f(λ)=σ2(1-j=1pαjexp(-2πιjλ))2
wo αj sind die autoregressiven Parameter.

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

f(0)=σ2(1-j=1pα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)
tchakravarty
quelle
0

Überprüfen Sie die Wikipedia-Seite . Du wirst sehenSxx(ω), das ist die spektrale Dichte. In Ihrem Fall sollten Sie verwendenSxx(0).

xuhdev
quelle