Czy mogę częściowo zautomatyzować diagnostykę konwergencji MCMC, aby ustawić długość wygrzewania?

Chciałbym zautomatyzować wybór wypalania dla łańcucha MCMC, np. Usuwając pierwsze n wierszy na podstawie diagnostyki zbieżności.

W jakim stopniu można bezpiecznie zautomatyzować ten krok? Nawet jeśli nadal dokładnie sprawdzę autokorelację, ślad MCcm i pliki PDF, byłoby miło mieć automatyczny wybór długości wypalenia.

Moje pytanie jest ogólne, ale byłoby wspaniale, gdybyś mógł podać szczegóły postępowania z obiektem R mcmc.ob; Korzystam z pakietów rjags i coda w R.

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.

1. obliczyć wektor maksymalnej mediany współczynnika skurczu diagnostycznego Gelmana-Rubina (grsf) dla wszystkich zmiennych w
2. 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
3. podpróbkować łańcuchy od tego punktu do końca łańcucha
4. rozcieńczyć łańcuch za pomocą autokorelacji najbardziej autokorelowanego łańcucha
5. 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:

1. Nie można zdiagnozować konwergencji, a jedynie zdiagnozować brak konwergencji (Brooks, Giudici i Philippe, 2003)

2. 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.

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