Les bookmarks de Clem
2369 shaares
43 results
tagged
MCMC
Un package qui commence à pas mal faire parler de lui. À creuser.
Une approche pour bien sélectionner la proposal en MCMC
Dans SAS, la phase de tuning de la proposal dans l'algorithme de Metropolis consiste dans le cas multivarié, à estimer la matrice de covariance d'après les réalisations déjà obtenues. Le principe est le suivant: on a des tuning loops de 500 itérations par défaut chacunes. Lors de chaque loop, on lance une chaîne avec une proposal multinormale caractérisée par une matrice de covariance donnée. À la loop suivante, même principe, mais la matrice de covariance est mise à jour comme une moyenne pondérée de l'ancienne matrice de covariance et de la nouvelle matrice obtenue d'après les réalisations. En plus, un side product, c'est qu'on trouve le MAP de la posterior, dont on peut se servir comme valeurs initiales.
Très intéressant: en MCMC lorsque la proposal (e.g. une gaussienne) n'a pas le même support que la cible (e.g. une variable strictement positive), on est tenté de répéter l'échantillonnage de la proposal jusqu'à ce que la valeur proposée soit dans le support de la cible. Il est bien connu que cette approche n'est pas valide quand elle est implémentée de façon bourrine. Mais l'auteur montre ici qu'elle peut être valide, à condition d'utiliser la bonne distribution dans le calcul du rapport pour déterminer la proba d'acceptation.
M'a l'air super intéressant ce blog...
M'a l'air super intéressant ce blog...
Article de blog intéressant sur l'ajustement de spline 2D avec JAGS. La fonction jagam du package mgcv semble très intéressante. Décortiquer cette fonction, programmée par S.Wood lui-même (!), risque d'apporter pas mal d'infos super intéressantes!
Ben décidément, yen a des choses sur Nimble aujourd'hui
Plein de sources pour la mise en œuvre du RJMCMC...
Je stocke l'expérience ici car je viens de perdre une semaine à cause de cette connerie.
Quand on ajuste un modèle par MCMC avec un sampler de Gibbs (JAGS, Winbugs, etc.), toutes les variables pour lesquelles on va générer des valeurs jouent sur le comportement du sampler.
Ça a l'air évident à dire, mais je viens de me faire avoir comme un bleu en oubliant ça.
Je cherchais à ajuster un modèle estimant des limites de populations, avec estimation de probabilités de migration d'une pop à l'autre. Et je m'étais dit, tiens, ben comme j'ai tous les éléments nécessaires, pour chaque individu de mon échantillon, je vais échantillonner une appartenance de l'individu à une population dans la posterior. Pour chaque individu, je considère cette appartenance comme une variable cachée n'ayant aucune influence sur mon estimation. Cette appartenance dépend des paramètres estimés, du coup, à chaque itération, je génère une valeur et comme ça à la fin, j'ai déjà les dépendances dont je peux me servir directement à titre exploratoire.
Et j'avais pas pensé que ça marche aussi dans l'autre sens: le sampler de Gibbs, quand il va sampler dans la postérior conditionnelle pour les limites de population, il va le faire conditionnellement aux appartenances des individus à une population à l'itération t. Donc pour modifier une limite de population, il faut simultanément modifier les paramètres qui contrôlent cette limite ET les appartenance de tous les individus "frontaliers" concerné (changer leur "nationalité"). Or, la question du simultané avec un sampler de Gibbs, c'est très compliqué...
Résultat: mélange franchement dégueulasse, et résultat inutilisable.
Bon ben je le saurai pour la prochaine fois. Conclusion évidente toujours bonne à rappeler: On ne met dans le modèle à ajuster que les paramètres que l'on souhaite ajuster. Et rien d'autre.
Quand on ajuste un modèle par MCMC avec un sampler de Gibbs (JAGS, Winbugs, etc.), toutes les variables pour lesquelles on va générer des valeurs jouent sur le comportement du sampler.
Ça a l'air évident à dire, mais je viens de me faire avoir comme un bleu en oubliant ça.
Je cherchais à ajuster un modèle estimant des limites de populations, avec estimation de probabilités de migration d'une pop à l'autre. Et je m'étais dit, tiens, ben comme j'ai tous les éléments nécessaires, pour chaque individu de mon échantillon, je vais échantillonner une appartenance de l'individu à une population dans la posterior. Pour chaque individu, je considère cette appartenance comme une variable cachée n'ayant aucune influence sur mon estimation. Cette appartenance dépend des paramètres estimés, du coup, à chaque itération, je génère une valeur et comme ça à la fin, j'ai déjà les dépendances dont je peux me servir directement à titre exploratoire.
Et j'avais pas pensé que ça marche aussi dans l'autre sens: le sampler de Gibbs, quand il va sampler dans la postérior conditionnelle pour les limites de population, il va le faire conditionnellement aux appartenances des individus à une population à l'itération t. Donc pour modifier une limite de population, il faut simultanément modifier les paramètres qui contrôlent cette limite ET les appartenance de tous les individus "frontaliers" concerné (changer leur "nationalité"). Or, la question du simultané avec un sampler de Gibbs, c'est très compliqué...
Résultat: mélange franchement dégueulasse, et résultat inutilisable.
Bon ben je le saurai pour la prochaine fois. Conclusion évidente toujours bonne à rappeler: On ne met dans le modèle à ajuster que les paramètres que l'on souhaite ajuster. Et rien d'autre.
Première tentative d'ajustement de modèle bayésien sous STAN... échec: STAN ne permet pas l'ajustement de modèles avec des variables cachées entières (e.g. un effectif détecté comme réponse, un effectif réel comme paramètre). Ce qui rétrospectivement semble assez logique, quand on connaît le principe du Monte Carlo Hamiltonien. La seule solution est de marginaliser le paramètre entier quand on en a un... Ce qui n'est pas toujours simple à réaliser. Enfin, dans le cas présent, je n'ai pas le choix, mon modèle est caractérisé par un mélange moisi avec un sampler de Gibbs, et je pense que le fait de m'appuyer sur un paramètre latent entier n'y est pas pour rien.
Edit: le manuel de Stan, section 11.3, évoque les modèles de CMR. Il présente deux cas de figure:
* l'estimateur de Lincoln-Petersen de la taille de population N à partir d'animaux marqués à une première occasion de capture, puis recapturés à une deuxième occasion. Dans ce cas, on peut traiter le modèle en prenant N comme paramètre continu.
* le modèle de Cormack-Jolly-Seber, dans lequel on a une variable latente z_i(t) pour chaque animal i qui prend la valeur 1 si l'animal est vivant au temps t et 0 sinon. Pour pouvoir ajuster ce modèle, il faut marginaliser pour se débarrasser du paramètre.
Oui, donc Stan, uniquement quand on a des paramètres continus (éviter les variables latentes discrètes genre effectif réel non observé).
Edit: le manuel de Stan, section 11.3, évoque les modèles de CMR. Il présente deux cas de figure:
* l'estimateur de Lincoln-Petersen de la taille de population N à partir d'animaux marqués à une première occasion de capture, puis recapturés à une deuxième occasion. Dans ce cas, on peut traiter le modèle en prenant N comme paramètre continu.
* le modèle de Cormack-Jolly-Seber, dans lequel on a une variable latente z_i(t) pour chaque animal i qui prend la valeur 1 si l'animal est vivant au temps t et 0 sinon. Pour pouvoir ajuster ce modèle, il faut marginaliser pour se débarrasser du paramètre.
Oui, donc Stan, uniquement quand on a des paramètres continus (éviter les variables latentes discrètes genre effectif réel non observé).
Il y a même un mode emacs pour faire du STAN! Je sens que je vais vraiment tester l'outil très prochainement...
Tiens? une alternative à BUGS et JAGS? À suivre...
Très intéressant!
À suivre de près.
À suivre de près.
Une autre solution au slicer coincé à l'infini. Et une explication au problème...
grrrrrr
Edit: si le problème se produit, il ne faut pas oublier que la distribution gamma inclue une *fonction* gamma, et que Gamma(x) = (x-1)!
Et la factorielle augmente assez fortement avec x. Donc pour une distribution gamma G(a,b), si a est trop important, Gamma(a) ne pourra être calculé numériquement, et ça se traduira par une valeur infinie. Un bidouillage consiste à tronquer la prior de a si a est une variable à estimer.
grrrrrr
Edit: si le problème se produit, il ne faut pas oublier que la distribution gamma inclue une *fonction* gamma, et que Gamma(x) = (x-1)!
Et la factorielle augmente assez fortement avec x. Donc pour une distribution gamma G(a,b), si a est trop important, Gamma(a) ne pourra être calculé numériquement, et ça se traduira par une valeur infinie. Un bidouillage consiste à tronquer la prior de a si a est une variable à estimer.
La loi gamma est une loi de merde. Mais malheureusement, parfois inévitable.
Je me stocke ça ici pour si jamais un jour je retombe sur le même problème. Quand on ajuste un modèle de tendance temporelle avec des résidus autocorrélés selon un modèle de corrélation exponentielle (Diggle et al. 2002, Analysis of longitudinal data, p. 56) avec JAGS, il ne faut pas générer les résidus en tirant dans une loi multinormale de moyenne nulle que l'on va ajouter à une tendance. Il faut tirer dans une loi multinormale dont le vecteur moyenne EST la tendance. Je me suis arraché les cheveux là-dessus, mais le mélange est bien meilleur si l'on procède comme ça. Autrement dit, il ne faut pas faire ça:
residus~dmnorm(vecteurDeZeros, Omega)
for (i in 1:10) {
esperance[i]<- mu[i] + residus[i]
}
Mais plutôt faire ça:
esperance~dmnorm(mu, Omega)
Les chaînes se mélangent mieux. Fallait le savoir...
Edit: Bon apparemment, la stratégie est connue et est valable pour tout modèle linéaire réclamant un résidu quelconque: on ne doit jamais tirer au sort un résidu dans une loi normale de moyenne nulle. Mais on indique que la réponse suit une loi normale de moyenne correspondant à l'espérance modélisée. Ça s'appelle du hierarchical centring (voir Browne et al. 2009, dans Journal of the Royal Statistical Society). Par contre, il paraît que ça marche pas top quand la variance est faible.
residus~dmnorm(vecteurDeZeros, Omega)
for (i in 1:10) {
esperance[i]<- mu[i] + residus[i]
}
Mais plutôt faire ça:
esperance~dmnorm(mu, Omega)
Les chaînes se mélangent mieux. Fallait le savoir...
Edit: Bon apparemment, la stratégie est connue et est valable pour tout modèle linéaire réclamant un résidu quelconque: on ne doit jamais tirer au sort un résidu dans une loi normale de moyenne nulle. Mais on indique que la réponse suit une loi normale de moyenne correspondant à l'espérance modélisée. Ça s'appelle du hierarchical centring (voir Browne et al. 2009, dans Journal of the Royal Statistical Society). Par contre, il paraît que ça marche pas top quand la variance est faible.
Ah ben je suis pas le seul à pas aimer WinBUGS et OpenBUGS.
Utiliser JAGS pour ajuster du modèle avec autocorrélation spatiale. À garder sous le coude pour le jour où...
Qui pourrait bien arriver plus rapidement que prévu (je crains cependant qu'on ne soit limité en termes de nombre d'unités spatiales que l'on peut intégrer le modèle).
À lire un jour
Qui pourrait bien arriver plus rapidement que prévu (je crains cependant qu'on ne soit limité en termes de nombre d'unités spatiales que l'on peut intégrer le modèle).
À lire un jour
Faudrait que je prenne le temps de creuser c't'affaire. Pas mal de références fournies.
Vinzou, faut vraiment que je prenne le temps de creuser ce logiciel. Ça m'a l'air d'être une alternative intéressante à JAGS pour le MCMC. Plus rapide apparemment... Et ya même une interface pour R