Я пытаюсь повторить расчет ESS (эффективный размер выборки), используя метод Vehtari et al. в: Нормализация рангов, свертывание и локализация: улучшенная Rhat для оценки сходимости MCM C
Я работаю с кодом здесь: https://github.com/avehtari/rhat_ess/blob/master/code/monitornew.R
# Geyer's initial positive sequence
rho_hat_t <- rep.int(0, n_samples)
t <- 0
rho_hat_even <- 1
rho_hat_t[t + 1] <- rho_hat_even
rho_hat_odd <- 1 - (mean_var - mean(acov[t + 2, ])) / var_plus # 251
rho_hat_t[t + 2] <- rho_hat_odd
while (t < nrow(acov) - 5 && !is.nan(rho_hat_even + rho_hat_odd) &&
(rho_hat_even + rho_hat_odd > 0)) {
t <- t + 2
rho_hat_even = 1 - (mean_var - mean(acov[t + 1, ])) / var_plus # 256
rho_hat_odd = 1 - (mean_var - mean(acov[t + 2, ])) / var_plus # 257
if ((rho_hat_even + rho_hat_odd) >= 0) {
rho_hat_t[t + 1] <- rho_hat_even
rho_hat_t[t + 2] <- rho_hat_odd
}
}
Я могу следовать коду из статьи, кроме случаев, когда мы дойдем до уравнения 10 в статье (вычисление перекрестной автокорреляции). Код (строки 251, 256 и 257) отображается в виде:
1 - (mean_var - mean(acov[t + 1, ])) / var_plus
, который близок к уравнению 10, за исключением пропущенных в уравнении 10 членов 's':
Я не вижу нигде в коде, что это каким-то образом учитывается в другом месте в способе вычисления. Я попытался вставить термины «s» обратно в эти строки кода, и это сильно влияет на окончательное значение ESS.
Может ли кто-нибудь помочь мне понять несоответствие между бумагой и кодом?
Спасибо.