Diagnostik MCMC Geweke

Saya menjalankan sampler Metropolis (C ++) dan ingin menggunakan sampel sebelumnya untuk memperkirakan tingkat konvergensi.

Satu yang mudah untuk mengimplementasikan diagnostik yang saya temukan adalah diagnostik Geweke , yang menghitung perbedaan antara dua rata-rata sampel dibagi dengan kesalahan standar yang diestimasi. Kesalahan standar diperkirakan dari kepadatan spektral di nol.

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

di mana $A$ , adalah dua jendela dalam rantai Markov. Saya melakukan penelitian tentang apa yang dan tetapi masuk ke dalam literatur literatur tentang kepadatan spektral energi dan spektral daya kepadatan tapi aku bukan ahli tentang topik ini; Saya hanya perlu jawaban cepat: apakah jumlah ini sama dengan varians sampel? Jika tidak, apa rumus untuk menghitungnya?$B$$\hat{S_{\theta}^A}(0)$$\hat{S_{\theta}^B}(0)$

Keraguan lain pada diagnostik Geweke ini adalah bagaimana cara memilih ? Literatur di atas mengatakan bahwa itu adalah beberapa fungsional dan harus menyiratkan adanya kepadatan spektral , tetapi untuk kenyamanan saya kira cara paling sederhana adalah dengan menggunakan fungsi identitas (gunakan sampel itu sendiri). Apakah ini benar?$\theta$$\theta(X)$$\hat{S_{\theta}^A}(0)$

R paket coda memiliki deskripsi tapi juga tidak menentukan bagaimana untuk menghitung nilai-nilai.$S$

Anda dapat melihat kode untuk geweke.diagfungsi dalam codapaket untuk melihat bagaimana varians dihitung, melalui panggilan ke spectrum.ar0fungsi.

---

Berikut ini adalah motivasi singkat dari perhitungan kepadatan spektral dari proses AR ( ) pada nol.$p$

$p$$\lambda$

$$f(\lambda) = \dfrac{\sigma^2}{(1-\sum_{j=1}^p\alpha_j\exp(-2\pi\iota j\lambda))^2}$$ $\alpha_j$

$p$$0$

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

Perhitungan kemudian akan terlihat seperti ini (menggantikan penduga yang biasa dengan 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)

$S_{xx}(\omega)$$S_{xx}(0)$