Calcul de la vraisemblance marginale à partir d'échantillons MCMC

24

C'est une question récurrente (voir cet article , cet article et cet article ), mais j'ai un tour différent.

Supposons que j'ai un tas d'échantillons d'un échantillonneur MCMC générique. Pour chaque échantillon , je connais la valeur du log vraisemblance et du log prior . Si cela aide, je connais également la valeur de la vraisemblance du journal par point de données, (ces informations aident avec certaines méthodes, telles que WAIC et PSIS-LOO).log f ( x | θ ) log f ( θ ) log f ( x i | θ )θlogf(x|θ)logf(θ)logf(xi|θ)

Je veux obtenir une estimation (brute) de la probabilité marginale, juste avec les échantillons que j'ai, et éventuellement quelques autres évaluations de fonctions (mais sans relancer un MCMC ad hoc ).

Tout d'abord, effaçons le tableau. Nous savons tous que l'estimateur harmonique est le pire estimateur de tous les temps . Allons-nous en. Si vous faites un échantillonnage de Gibbs avec des prieurs et des postérieurs sous forme fermée, vous pouvez utiliser la méthode de Chib ; mais je ne sais pas comment généraliser en dehors de ces cas. Il existe également des méthodes qui vous obligent à modifier la procédure d'échantillonnage (par exemple via des postérieurs trempés ), mais cela ne m'intéresse pas ici.

L'approche à laquelle je pense consiste à approximer la distribution sous-jacente avec une forme paramétrique (ou non paramétrique) , puis à déterminer la constante de normalisation comme un problème d'optimisation 1-D (c'est-à-dire le qui minimise certaines erreurs) entre et , évalué sur les échantillons). Dans le cas le plus simple, supposons que le postérieur soit à peu près normal à plusieurs variables, je peux adapter comme une normale à plusieurs variables et obtenir quelque chose de similaire à une approximation de Laplace (je pourrais vouloir utiliser quelques évaluations de fonctions supplémentaires pour affiner la position de la mode). Cependant, je pourrais utiliser commeZ Z Z g ( θ ) f ( x | θ ) f ( θ ) g ( θ ) g ( θ )g(θ)ZZZg(θ)f(x|θ)f(θ)g(θ)g(θ)une famille plus flexible telle qu'un mélange variationnel de distributions multivariées .t

J'apprécie que cette méthode ne fonctionne que si est une approximation raisonnable de , mais toute raison ou mise en garde expliquant pourquoi il serait très imprudent de fais le? Une lecture que vous recommanderiez?f ( x | θ ) f ( θ )Zg(θ)f(x|θ)f(θ)

L'approche entièrement non paramétrique utilise une famille non paramétrique, comme un processus gaussien (GP), pour approximer (ou une autre transformation non linéaire de celui-ci, telle que comme racine carrée) et la quadrature bayésienne à intégrer implicitement sur la cible sous-jacente (voir ici et ici ). Cela semble être une approche alternative intéressante, mais analogue dans l'esprit (notez également que les généralistes seraient difficiles à manier dans mon cas).logf(x|θ)+logf(θ)

lacerbi
la source
6
Je pense que Chib, S. et Jeliazkov, I. 2001 "La probabilité marginale de la sortie de Metropolis - Hastings" se généralise aux sorties normales de MCMC - serait intéressé à entendre des expériences avec cette approche. Quant au GP - en gros, cela se résume à l'émulation du postérieur, que vous pourriez également envisager pour d'autres problèmes. Je suppose que le problème est que vous n'êtes jamais sûr de la qualité de l'approximation. Ce que je me demande aussi, c'est si un échantillon MCMC est idéal pour un modèle GP, ou si vous devriez investir davantage dans les queues.
Florian Hartig
2
(+1) Merci pour la référence, semble parfait - je vais le vérifier. Je suis d'accord que toutes les approches basées sur un modèle peuvent être problématiques (la bonne chose avec la quadrature bayésienne est que vous obtenez une estimation de l'incertitude, bien que vous ne soyez pas sûr de son calibrage). Pour le moment, mon objectif modeste est de faire quelque chose de "meilleur qu'une approximation de Laplace".
lacerbi

Réponses:

26

L'extension de Chib et Jeliazkov (2001) devient malheureusement rapidement coûteuse ou très variable, ce qui explique pourquoi elle n'est pas très utilisée en dehors des cas d'échantillonnage de Gibbs.

Bien qu'il existe de nombreuses façons et approches pour résoudre le problème d'estimation de la constante de normalisation (comme l'illustrent les discussions assez diverses de l' atelier d'estimation des constantes que nous avons organisé la semaine dernière à l'Université de Warwick, des diapositives y sont disponibles ), certaines solutions exploitent directement la sortie MCMC .Z

  1. Comme vous l'avez mentionné, l'estimateur de la moyenne harmonique de Newton et Raftery (1994) est presque toujours médiocre pour avoir une variance infinie. Cependant, il existe des moyens d'éviter la malédiction de variance infinie en utilisant à la place une cible de support fini dans l'identité moyenne harmonique en choisissantαcomme indicateur d'une région HPD pour la partie postérieure. Cela garantit une variance finie en supprimant les queues dans la moyenne harmonique. (Les détails se trouvent dansun article que j'ai écrit avec Darren Wraithet dans unchapitre sur la normalisation des constantesécrit avec Jean-Michel Marin.) En bref, la méthode recycle la sortie MCMCθ1,,θMen identifiant leβ( 20% disent) les plus grandes valeurs de la cibleπ(θ)f(x|θ)et créantα

    α(θ)π(θ)f(x|θ)dπ(θ|x)=1Z
    αθ1,,θMβπ(θ)f(x|θ)αcomme simulations uniforme sur l'union des billes centré à celles plus grande densité (HPD) et de rayon ρ , ce qui signifie l'estimation de la constante de normalisation Z est donnée par Z - 1 = 1θi0ρZ sidest la dimension deθ(les corrections s'appliquent pour les boules qui se croisent) et siρest suffisamment petit pour que les boules ne se coupent jamais (ce qui signifie qu'au mieux un seul indicateur sur les boules est différent de zéro). L'explication dudénominateurαM2est qu'il s'agit d'une double somme determesβM2: 1
    Z^1=1βM2m=1Mdouble sum overβM ball centres θi0and M simulations θmI(0,ρ)(mini||θmθi0||){π(θm)f(x|θm)}1/πd/2ρdΓ(d/2+1)1volume of ball with radius ρβMα(θm)π(θm)f(x|θm)
    dθραM2βM2 avec chaque terme dansθmintégrant àZ-1.
    1βMi=1βM1Mm=1MU(θi0,ρ)(θm)same as with min×1π(θm)f(x|θm)
    θmZ1
  2. Z

    i=1nf(xi|θ)nlogexpf(x|θ)dx
    i=1n[f(xi|θ)+ν]nexp[f(x|θ)+ν]dx
    which is the log-likelihood of a Poisson point process with intensity function
    exp{f(x|θ)+ν+logn}
    This is an alternative model in that the original likelihood does not appear as a marginal of the above. Only the modes coincide, with the conditional mode in ν providing the normalising constant. In practice, the above Poisson process likelihood is unavailable and Guttmann and Hyvärinen (2012) offer an approximation by means of a logistic regression. To connect even better with your question, Geyer's estimate is a MLE, hence solution to a maximisation problem.
  3. A connected approach is Charlie Geyer's logistic regression approach. The fundamental notion is to add to the MCMC sample from π(θ|x) another sample from a known target, e.g., your best guess at π(θ|x), g(θ), and to run logistic regression on the index of the distribution behind the data (1 for π(θ|x) and 0 for g(θ)). With the regressors being the values of both densities, normalised or not. This happens to be directly linked with Gelman and Meng (1997) bridge sampling, which also recycles samples from different targets. And later versions, like Meng's MLE.
  4. A different approach that forces one to run a specific MCMC sampler is Skilling's nested sampling. While I [and others] have some reservations on the efficiency of the method, it is quite popular in astrostatistics and cosmology, with software available like multinest.
  5. A last [potential if not always possible] solution is to exploit the Savage-Dickey representation of the Bayes factor in the case of an embedded null hypothesis. If the null writes as H0:θ=θ0 about a parameter of interest and if ξ is the remaining [nuisance] part of the parameter of the model, assuming a prior of the form π1(θ)π2(ξ), the Bayes factor of H0 against the alternative writes as
    B01(x)=πθ(θ0|x)π1(θ0)
    where πθ(θ0|x) denotes the marginal posterior density of θ at the specific value θ0. In case the marginal density under the null H0:θ=θ0
    m0(x)=Ξf(x|θ0,ξ)π2(ξ)dξ
    is available in closed form, one can derive the marginal density for the unconstrained model
    ma(x)=Θ×Ξf(x|θ,ξ)π1(θ)π2(ξ)dθdξ
    from the Bayes factor. (This Savage-Dickey representation relies on specific versions of three different densities and so is fraught with danger, not even mentioning the computational challenge of producing the marginal posterior.)

[Here is a set of slides I wrote about estimating normalising constants for a NIPS workshop last December.]

Xi'an
la source
2
(+1) Incredibly rich answer, thank you. This will be useful to me and, I suppose, many other people. It will take me some time to have a look at the various approaches, and then I might come back with specific questions.
lacerbi
2
Starting from point (1)... I read the relevant articles. The "corrected" harmonic mean estimator seems exactly what I was looking for. It's neat and easy to compute given a MCMC output. So... what's the catch? It doesn't look like the method is being widely used, judging from a quick search on Google Scholar. What are its limitations? (besides the need to identify the HPD regions, which I imagine might become an issue for very complicated posteriors in high dimension). I am definitely going to give it a try -- but I wonder if there is something I need to be wary of.
lacerbi
2
I added a few more details: the issue in implementing the HPD uniform is to figure out a proper compact approximation for the HPD region. The convex hull of points with high posterior values is (NP?) hard to determine while balls centred at those points may intersect, which creates a secondary normalising constant problem.
Xi'an
2
@Xi'an : very helpful, thanks! Can I ask: of all the mentioned approaches, what would currently be your recommendation if one looks for a general approach that tends to work out of the box (i.e. no tuning / checking required from the user)? I would be especially interested in the case of models with a low (< 50) number of parameters, non-normal posteriors, and strong correlations between parameters.
Florian Hartig
1
@FlorianHartig: the fact that a generic software like BUGS does not return a generic estimate of Z is sort of revealing the extent of the problem. The many solutions that one can find in the specialised literature have not produced a consensus estimate. Hence, my recommendation would be to opt for Geyer's logistic regression solution, which is somewhat insensitive to dimension.
Xi'an