机器算法验证 - MCMC格维克诊断 - 吾爱随笔录

MCMC格维克诊断

机器算法验证马尔可夫链蒙特卡罗诊断

2022-02-28 23:28:53

我正在运行 Metropolis 采样器 (C++)，并希望使用之前的样本来估计收敛速度。

我发现一个易于实施的诊断是Geweke 诊断，它计算两个样本均值之间的差异除以其估计的标准误差。标准误差是从零处的光谱密度估计的。

Z_{n} = \frac{{\bar{θ}}_{A} - {\bar{θ}}_{B}}{\sqrt{\frac{1}{n_{A}} \hat{S_{θ}^{A}} (0) + \frac{1}{n_{B}} \hat{S_{θ}^{B}} (0)}},

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

在哪里 $A$ , $B$ 是马尔可夫链中的两个窗口。我做了一些研究是什么 $\hat{S_{\theta}^A}(0)$ 和 $\hat{S_{\theta}^B}(0)$ 但是进入一堆关于能量谱密度和功率谱密度的文献，但我不是这些主题的专家；我只需要一个快速的答案：这些数量与样本方差相同吗？如果不是，计算它们的公式是什么？

对这个 Geweke 诊断的另一个疑问是如何选择 $\theta$ ? 上面的文献说是一些功能性的 $\theta(X)$ 并且应该暗示存在谱密度 $\hat{S_{\theta}^A}(0)$ ，但为方便起见，我想最简单的方法是使用标识函数（使用样本本身）。它是否正确？

R coda 包有一个描述，但它也没有指定如何计算 $S$ 价值观。

2个回答

您可以通过对geweke.diag函数的调用查看包中函数的代码，coda以了解如何计算方差spectrum.ar0。

这是计算 AR 的谱密度的简短动机（ $p$ ) 过程为零。

AR的光谱密度（ $p$ ) 频率处理 $\lambda$ 由表达式给出：

f (λ) = \frac{σ^{2}}{(1 - \sum_{j = 1}^{p} α_{j} \exp (- 2 π ι j λ))^{2}}

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

α_{j}

$\alpha_j$ 是自回归参数。

在计算 AR 的谱密度时，这个表达式大大简化了（ $p$ ) 处理在 $0$ ：

f (0) = \frac{σ^{2}}{(1 - \sum_{j = 1}^{p} α_{j})^{2}}

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

然后计算看起来像这样（用通常的估计器代替参数）：

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)$ .

其它你可能感兴趣的问题