Puis-je semi-automatiser les diagnostics de convergence MCMC pour définir la longueur de rodage?

13

Je voudrais automatiser le choix du burn-in pour une chaîne MCMC, par exemple en supprimant les n premières lignes sur la base d'un diagnostic de convergence.

Dans quelle mesure cette étape peut-elle être automatisée en toute sécurité? Même si je vérifie toujours l'autocorrélation, la trace mcmc et les fichiers PDF, il serait bien d'avoir le choix de la durée de gravure automatisé.

Ma question est générale, mais ce serait formidable si vous pouviez fournir des détails pour traiter un objet R mcmc.object; J'utilise les packages rjags et coda dans R.

David LeBauer
la source
bien qu'il ne soit pas inclus dans la question d'origine, il serait également utile de définir automatiquement l'intervalle d'amincissement comme proposé dans ma réponse.
David LeBauer
1
Je voudrais juste mentionner qu'en tant que personne intéressée à faire des algorithmes MCMC génériques, facilement applicables à de nombreux problèmes, je suis très intéressé par ce sujet.
John Salvatier

Réponses:

6

Voici une approche à l'automatisation. Commentaires très appréciés. Il s'agit d'une tentative de remplacer l'inspection visuelle initiale par un calcul, suivie d'une inspection visuelle ultérieure, conformément à la pratique standard.

Cette solution intègre en fait deux solutions potentielles, tout d'abord, calculer le rodage pour supprimer la longueur de la chaîne avant qu'un certain seuil ne soit atteint, puis utiliser la matrice d'autocorrélation pour calculer l'intervalle d'amincissement.

  1. calculer un vecteur du facteur de rétrécissement diagnostique de convergence Gelman-Rubin médian maximal (grsf) pour toutes les variables
  2. trouver le nombre minimum d'échantillons auquel le grsf à travers toutes les variables passe en dessous d'un certain seuil, par exemple 1,1 dans l'exemple, peut-être plus bas en pratique
  3. sous-échantillonner les chaînes de ce point à la fin de la chaîne
  4. amincir la chaîne en utilisant l'autocorrélation de la chaîne la plus autocorrélée
  5. confirmer visuellement la convergence avec les tracés de trace, d'autocorrélation et de densité

L'objet mcmc peut être téléchargé ici: 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)

--mise à jour--

Comme implémenté dans R, le calcul de la matrice d'autocorrélation est plus lent que ce qui serait souhaitable (> 15 min dans certains cas), dans une moindre mesure, tout comme le calcul du facteur de rétrécissement GR. Il y a une question sur la façon d'accélérer l'étape 4 sur stackoverflow ici

--update partie 2--

réponses supplémentaires:

  1. Il n'est pas possible de diagnostiquer la convergence, seulement de diagnostiquer le manque de convergence (Brooks, Giudici et Philippe, 2003)

  2. La fonction autorun.jags du package runjags automatise le calcul de la longueur d'exécution et les diagnostics de convergence. Il ne commence à surveiller la chaîne que lorsque le diagnostic de la rubine Gelman est inférieur à 1,05; il calcule la longueur de chaîne en utilisant le diagnostic Raftery et Lewis.

  3. Gelman et al (Gelman 2004 Bayesian Data Analysis, p. 295, Gelman et Shirley, 2010 ) déclarent qu'ils utilisent une approche conservatrice de rejet de la première moitié de la chaîne. Bien qu'il s'agisse d'une solution relativement simple, en pratique, cela suffit pour résoudre le problème pour mon ensemble particulier de modèles et de données.


#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)
David LeBauer
la source
2
Deux principes s'appliquent: vous ne pouvez jamais savoir si votre chaîne a convergé vers sa distribution stationnaire. Et tout test de convergence que vous pouvez faire manuellement, vous pouvez automatiser. Votre approche semble donc suffisamment solide.
Tristan
Dans la documentation runjags, je vois que autorun.jags indique que le modèle est automatiquement évalué pour la convergence et la taille d'échantillon adéquate avant d'être renvoyé. Pourriez-vous m'indiquer où vous avez constaté que autorun.jags ne commence pas à surveiller la chaîne jusqu'à ce que le diagnostic de la rubine Gelman soit inférieur à 1,05? Merci
user1068430
@ user1068430 in autorun.jags, ...permet de transmettre les paramètres à la add.summaryfonction. La add.summaryfonction a un argument psrf.targetavec une valeur par défaut de 1,05
David LeBauer