Kann ich die MCMC-Konvergenzdiagnose halbautomatisch durchführen, um die Einbrennlänge festzulegen?

Ich möchte die Auswahl des Einbrennens für eine MCMC-Kette automatisieren, z. B. durch Entfernen der ersten n Zeilen basierend auf einer Konvergenzdiagnose.

Inwieweit kann dieser Schritt sicher automatisiert werden? Selbst wenn ich die Autokorrelation, die mcmc-Ablaufverfolgung und die pdfs noch zweimal überprüfe, wäre es schön, wenn die Auswahl der Einbrennlänge automatisiert wäre.

Meine Frage ist allgemein, aber es wäre großartig, wenn Sie Einzelheiten zum Umgang mit einem R mcmc.object angeben könnten. Ich benutze die Pakete rjags und coda in R.

Hier ist ein Ansatz bei der Automatisierung. Feedback sehr geschätzt. Dies ist ein Versuch, die anfängliche Sichtprüfung durch eine Berechnung zu ersetzen, gefolgt von einer anschließenden Sichtprüfung gemäß der Standardpraxis.

Diese Lösung enthält tatsächlich zwei mögliche Lösungen: Zuerst das Einbrennen berechnen, um die Kettenlänge zu entfernen, bevor ein bestimmter Schwellenwert erreicht ist, und dann die Autokorrelationsmatrix verwenden, um das Ausdünnungsintervall zu berechnen.

1. Berechnen Sie einen Vektor des maximalen Medians des diagnostischen Gelman - Rubin - Konvergenzschrumpffaktors (grsf) für alle Variablen in der
2. Finden Sie die minimale Anzahl von Samples, bei denen die grsf über alle Variablen unter einen bestimmten Schwellenwert fällt, z. B. 1,1 im Beispiel, möglicherweise in der Praxis niedriger
3. Nehmen Sie eine Teilprobe der Ketten von diesem Punkt bis zum Ende der Kette
4. Die Kette mit der Autokorrelation der am meisten autokorrelierten Kette verdünnen
5. Bestätigen Sie die Konvergenz visuell mit Kurven-, Autokorrelations- und Dichtediagrammen

Das mcmc-Objekt kann hier heruntergeladen werden: 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)

--aktualisieren--

Wie in R implementiert, ist die Berechnung der Autokorrelationsmatrix langsamer als wünschenswert (in einigen Fällen> 15 min), in geringerem Maße auch die Berechnung des GR-Schrumpffaktors. Es gibt eine Frage, wie Sie Schritt 4 beim Stackoverflow beschleunigen können hier

--Update Teil 2--

zusätzliche Antworten:

1. Es ist nicht möglich, eine Konvergenz zu diagnostizieren, sondern nur einen Mangel an Konvergenz (Brooks, Giudici und Philippe, 2003).

2. Die Funktion autorun.jags aus dem Paket runjags automatisiert die Berechnung der Lauflängen- und Konvergenzdiagnose. Die Überwachung der Kette wird erst gestartet, wenn die Gelman-Rubindiagnose unter 1,05 liegt. Er berechnet die Kettenlänge mithilfe der Raftery- und Lewis-Diagnose.

3. Gelman et al. (Gelman 2004 Bayesian Data Analysis, S. 295, Gelman und Shirley, 2010 ) geben an, dass sie einen konservativen Ansatz zum Verwerfen der ersten Hälfte der Kette verwenden. Obwohl dies eine relativ einfache Lösung ist, reicht dies in der Praxis aus, um das Problem für meinen speziellen Satz von Modellen und Daten zu lösen.

#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)