Oto jedno podejście do automatyzacji. Informacje zwrotne bardzo mile widziane. Jest to próba zastąpienia wstępnej kontroli wzrokowej obliczeniami, a następnie kolejnej kontroli wzrokowej, zgodnie ze standardową praktyką.
To rozwiązanie zawiera dwa potencjalne rozwiązania, po pierwsze, oblicz wypalenie, aby usunąć długość łańcucha przed osiągnięciem pewnego progu, a następnie użyj matrycy autokorelacyjnej, aby obliczyć przedział przerzedzania.
- obliczyć wektor maksymalnej mediany współczynnika skurczu diagnostycznego Gelmana-Rubina (grsf) dla wszystkich zmiennych w
- znaleźć minimalną liczbę próbek, przy których grsf we wszystkich zmiennych spada poniżej pewnego progu, np. 1,1 w przykładzie, być może mniej w praktyce
- podpróbkować łańcuchy od tego punktu do końca łańcucha
- rozcieńczyć łańcuch za pomocą autokorelacji najbardziej autokorelowanego łańcucha
- wizualnie potwierdź zbieżność za pomocą wykresów śledzenia, autokorelacji i gęstości
Obiekt mcmc można pobrać tutaj: jags.out.Rdata
# jags.out is the mcmc.object with m variables
library(coda)
load('jags.out.Rdata')
# 1. calculate max.gd.vec,
# max.gd.vec is a vector of the maximum shrink factor
max.gd.vec <- apply(gelman.plot(jags.out)$shrink[, ,'median'], 1, max)
# 2. will use window() to subsample the jags.out mcmc.object
# 3. start window at min(where max.gd.vec < 1.1, 100)
window.start <- max(100, min(as.numeric(names(which(max.gd.vec - 1.1 < 0)))))
jags.out.trunc <- window(jags.out, start = window.start)
# 4. calculate thinning interval
# thin.int is the chain thin interval
# step is very slow
# 4.1 find n most autocorrelated variables
n = min(3, ncol(acm))
acm <- autocorr.diag(jags.out.trunc)
acm.subset <- colnames(acm)[rank(-colSums(acm))][1:n]
jags.out.subset <- jags.out.trunc[,acm.subset]
# 4.2 calculate the thinning interval
# ac.int is the time step interval for autocorrelation matrix
ac.int <- 500 #set high to reduce computation time
thin.int <- max(apply(acm2 < 0, 2, function(x) match(T,x)) * ac.int, 50)
# 4.3 thin the chain
jags.out.thin <- window(jags.out.trunc, thin = thin.int)
# 5. plots for visual diagnostics
plot(jags.out.thin)
autocorr.plot(jags.win.out.thin)
--aktualizacja--
Jak zaimplementowano w R, obliczanie macierzy autokorelacji jest wolniejsze niż byłoby to pożądane (> 15 minut w niektórych przypadkach), w mniejszym stopniu, podobnie jak obliczanie współczynnika kurczenia GR. Jest to kwestia o tym, jak przyspieszyć kroku 4 na stackoverflow tutaj
- aktualizacja część 2--
dodatkowe odpowiedzi:
Nie można zdiagnozować konwergencji, a jedynie zdiagnozować brak konwergencji (Brooks, Giudici i Philippe, 2003)
Funkcja autorun.jags z pakietu runjags automatyzuje obliczanie diagnostyki długości przebiegu i zbieżności. Nie rozpoczyna monitorowania łańcucha, dopóki diagnostyka rubinowa Gelmana nie spadnie poniżej 1,05; oblicza długość łańcucha za pomocą diagnostyki Raftery i Lewis.
Gelman i wsp. (Gelman 2004 Bayesian Data Analysis, s. 295, Gelman i Shirley, 2010 ) twierdzą, że stosują konserwatywne podejście do odrzucania pierwszej połowy łańcucha. Chociaż jest to stosunkowo proste rozwiązanie, w praktyce wystarczy, aby rozwiązać problem z moim konkretnym zestawem modeli i danych.
#code for answer 3
chain.length <- summary(jags.out)$end
jags.out.trunc <- window(jags.out, start = chain.length / 2)
# thin based on autocorrelation if < 50, otherwise ignore
acm <- autocorr.diag(jags.out.trunc, lags = c(1, 5, 10, 15, 25))
# require visual inspection, check acceptance rate
if (acm == 50) stop('check acceptance rate, inspect diagnostic figures')
thin.int <- min(apply(acm2 < 0, 2, function(x) match(TRUE, x)), 50)
jags.out.thin <- window(jags.out.trunc, thin = thin.int)