Berechnung der Grenzwahrscheinlichkeit aus MCMC-Proben

24

Dies ist eine wiederkehrende Frage (siehe diesen Beitrag , diesen Beitrag und diesen Beitrag ), aber ich habe einen anderen Dreh.

Angenommen, ich habe ein paar Samples von einem generischen MCMC-Sampler. Für jede Probe , weiß ich den Wert der Log - Likelihood und des log vor . Wenn es hilft, kenne ich auch den Wert der Log-Wahrscheinlichkeit pro Datenpunkt, (diese Information hilft bei bestimmten Methoden, wie WAIC und PSIS-LOO).θlogf(x|θ)logf(θ)logf(xi|θ)

Ich möchte eine (grobe) Schätzung der Grenzwahrscheinlichkeit erhalten, nur mit den Stichproben, die ich habe, und möglicherweise ein paar anderen Funktionsbewertungen (aber ohne ein Ad-hoc- MCMC erneut auszuführen ).

Lassen Sie uns zuerst den Tisch abräumen. Wir alle wissen, dass der harmonische Schätzer der schlechteste Schätzer aller Zeiten ist . Lass uns weitermachen. Wenn Sie Gibbs-Sampling mit Priors und Posteriors in geschlossener Form durchführen, können Sie die Chib-Methode verwenden . Ich bin mir jedoch nicht sicher, wie ich außerhalb dieser Fälle verallgemeinern soll. Es gibt auch Methoden, bei denen Sie das Probenahmeverfahren ändern müssen (z. B. über temperierte Seitenzähne ), aber das interessiert mich hier nicht.

Der Ansatz, an den ich denke, besteht darin, die zugrunde liegende Verteilung mit einer parametrischen (oder nichtparametrischen) Form ) zu approximieren und dann die Normalisierungskonstante als ein 1-D-Optimierungsproblem (dh das , das einen Fehler minimiert herauszufinden zwischen und , bewertet an den Abtastwerten). Nehmen wir im einfachsten Fall an, der posterior ist ungefähr multivariate Normalen, ich kann als multivariate Normalen anpassen und etwas Ähnliches wie eine Laplace-Näherung erhalten (ich möchte möglicherweise ein paar zusätzliche Funktionsauswertungen verwenden, um die Position von zu verfeinern der Modus). Ich könnte jedoch alsg(θ)ZZZg(θ)f(x|θ)f(θ)g(θ)g(θ)eine flexiblere Familie wie eine Variationsmischung multivariater Verteilungen.t

Ich schätze, dass diese Methode nur funktioniert, wenn eine vernünftige Annäherung an , aber jeder Grund oder jede warnende Geschichte, warum dies sehr unklug wäre TU es? Eine Lektüre, die Sie empfehlen würden?Zg(θ)f(x|θ)f(θ)

Der vollständig nichtparametrische Ansatz verwendet eine nichtparametrische Familie, wie einen Gaußschen Prozess (GP), um (oder eine andere nichtlineare Transformation davon, wie z als Quadratwurzel) und Bayes'sche Quadratur zur impliziten Integration über das zugrunde liegende Ziel (siehe hier und hier ). Dies scheint ein interessanter alternativer Ansatz zu sein, der jedoch sinngemäß analog ist (beachten Sie auch, dass Allgemeinmediziner in meinem Fall unhandlich wären).logf(x|θ)+logf(θ)

lacerbi
quelle
6
Ich denke, Chib, S. und Jeliazkov, I. 2001 "Geringe Wahrscheinlichkeit aus der Metropole - Hastings-Ausgabe" verallgemeinert auf normale MCMC-Ausgaben - wären interessiert, Erfahrungen mit diesem Ansatz zu hören. Was den GP betrifft, so läuft dies im Grunde auf die Emulation des Seitenzahns hinaus, die Sie auch für andere Probleme in Betracht ziehen könnten. Ich denke, das Problem ist, dass Sie sich über die Qualität der Approximation nie sicher sind. Was ich mich auch frage, ist, ob ein MCMC-Beispiel ideal für ein GP-Modell ist oder ob Sie mehr in die Schwänze investieren sollten.
Florian Hartig
2
(+1) Vielen Dank für den Hinweis, sieht genau richtig aus - ich werde es überprüfen. Ich stimme zu, dass alle modellbasierten Ansätze problematisch sein können (das Gute an der Bayes'schen Quadratur ist, dass Sie eine Schätzung der Unsicherheit erhalten, obwohl Sie nicht sicher sind, wie kalibriert sie ist). Im Moment ist es mein bescheidenes Ziel, etwas zu tun, das "besser als eine Laplace-Annäherung" ist.
Lacerbi

Antworten:

26

Die Erweiterung von Chib und Jeliazkov (2001) wird leider schnell kostspielig oder sehr variabel, was ein Grund dafür ist, dass sie außerhalb von Gibbs-Sampling-Fällen nicht häufig verwendet wird.

Während es viele Möglichkeiten und Ansätze für das Problem der Schätzung der Normalisierungskonstante (wie die recht unterschiedlichen Vorträge im Workshop zur Schätzung der Konstanten, den wir letzte Woche an der Universität von Warwick durchgeführt haben und die dort zur Verfügung stehen ), nutzen einige Lösungen die MCMC-Ausgabe direkt .Z

  1. Wie Sie bereits erwähnt haben, ist der harmonische Mittelwertschätzer von Newton and Raftery (1994) fast immer schlecht für eine unendliche Varianz. Es gibt jedoch Möglichkeiten, den Fluch der unendlichen Varianz zu vermeiden, indem stattdessen ein endliches Unterstützungsziel in der harmonischen mittleren Identität durch Auswählen vonαals Indikator für eine HPD-Region für den posterioren Bereich. Dies stellt eine endliche Varianz sicher, indem die Schwänze im harmonischen Mittelwert entfernt werden. (Details finden Sie ineinem Artikel, den ich mit Darren Wraith geschrieben habe,und in einemKapitel über das Normalisieren von Konstanten,die mit Jean-Michel Marin geschrieben wurden.) Kurz gesagt, die Methode recycelt die MCMC-Ausgabeθ1,,θMdurch Identifizieren desβ( 20% sagen) größte Werte des Zielsπ(θ)f(x|θ)und Erzeugen vonα

    α(θ)π(θ)f(x|θ)dπ(θ|x)=1Z
    αθ1,,θMβπ(θ)f(x|θ)αθich0ρZ , wennddie Dimension istθ(Korrekturen gelten für schneidende Kugeln) und wennρfür die Kugeln nie schneiden klein genug ist (was bedeutetdass am besten nureinem Indikator für die Kugeln von Null verschieden). Die Erklärung für denαM2-Nenner ist, dass dies eine doppelte Summe vonβM2-Termen ist: 1
    Z^-1=1βM2m=1Mdoppelte Summe überβM Ballzentren θich0und M Simulationen θmich(0,ρ)(Mindestich||θm-θich0||){π(θm)f(x|θm)}-1/πd/2ρdΓ(d/2+1)-1Volumen der Kugel mit Radius ρβMα(θm)π(θm)f(x|θm)
    dθραM2βM2 mit jedem Term in& theta;mzu integrierendenZ-1.
    1βMi=1βM1Mm=1MU(θi0,ρ)(θm)same as with min×1π(θm)f(x|θm)
    θmZ-1
  2. Ein anderer Ansatz besteht darin, die Normalisierungskonstante in einen Parameter umzuwandeln. Das klingt nach einer statistischen Ketzerei, aber die Arbeit von Guttmann und Hyvärinen (2012) hat mich vom Gegenteil überzeugt. Zu viel in die Details , ohne sich darin die nette Idee ist die beobachtete Log-Likelihood drehen n Σ i = 1 f ( x i | & thgr; ) - n log exp f ( x | & thgr; ) d x in einer gemeinsamen Log-Likelihood n i = 1 [ fZ

    ich=1nf(xich|θ)-nLogexpf(x|θ)dx
    was die logarithmische Wahrscheinlichkeit eines Poisson-Punkt-Prozesses mit der Intensitätsfunktion exp { f ( x | θ ) + ν + ist log n }
    ich=1n[f(xich|θ)+ν]-nexp[f(x|θ)+ν]dx
    exp{f(x|θ)+ν+Logn}
    Dies ist insofern ein alternatives Modell, als die ursprüngliche Wahrscheinlichkeit nicht als Randerscheinung des oben Gesagten erscheint. Nur die Modi stimmen überein, wobei der bedingte Modus in ν die Normalisierungskonstante liefert. In der Praxis ist die obige Wahrscheinlichkeit des Poisson-Prozesses nicht verfügbar, und Guttmann und Hyvärinen (2012) bieten eine Annäherung mittels einer logistischen Regression an. Um Ihre Frage noch besser beantworten zu können, ist Geyers Schätzung ein MLE und damit eine Lösung für ein Maximierungsproblem.
  3. π(θ|x)π(θ|x)G(θ)π(θ|x)G(θ)). Wenn die Regressoren die Werte beider Dichten sind, normalisiert oder nicht. Dies steht in direktem Zusammenhang mit der Brückenentnahme nach Gelman und Meng (1997), bei der auch Proben von verschiedenen Targets recycelt werden. Und spätere Versionen, wie Mengs MLE.
  4. Ein anderer Ansatz, der die Ausführung eines bestimmten MCMC-Samplers erzwingt, ist das verschachtelte Sampling von Skilling . Obwohl ich [und andere] einige Bedenken hinsichtlich der Effizienz der Methode habe, ist sie in der Astrostatistik und Kosmologie sehr beliebt, da Software wie Multinest verfügbar ist .
  5. H0:θ=θ0ξπ1(θ)π2(ξ)H0 wobeiπθ(θ0|x)die marginale hintere Dichte von bezeichnet
    B01(x)=πθ(θ0|x)π1(θ0)
    πθ(θ0|x)θθ0H0:θ=θ0
    m0(x)=Ξf(x|θ0,ξ)π2(ξ)dξ
    mein(x)=Θ×Ξf(x|θ,ξ)π1(θ)π2(ξ)dθdξ

[Hier ist eine Reihe von Folien, die ich über das Schätzen von Normalisierungskonstanten für einen NIPS-Workshop im vergangenen Dezember geschrieben habe.]

Xi'an
quelle
2
(+1) Unglaublich reichhaltige Antwort, danke. Dies wird mir und vermutlich vielen anderen Menschen von Nutzen sein. Ich werde einige Zeit brauchen, um mir die verschiedenen Ansätze anzuschauen, und dann werde ich möglicherweise mit spezifischen Fragen zurückkommen.
Lacerbi
2
Ausgehend von Punkt (1) ... habe ich die relevanten Artikel gelesen. Der "korrigierte" harmonische Mittelwertschätzer scheint genau das zu sein , wonach ich gesucht habe. Mit einer MCMC-Ausgabe ist es ordentlich und einfach zu berechnen. Also ... was ist der Haken? Es sieht nicht so aus, als ob die Methode weit verbreitet ist. Dies geht aus einer schnellen Suche in Google Scholar hervor. Was sind ihre Grenzen? (neben der Notwendigkeit, die HPD-Regionen zu identifizieren, die meiner Meinung nach ein Problem für sehr komplizierte Posterioren in großen Dimensionen sein könnten). Ich werde es auf jeden Fall versuchen - aber ich frage mich, ob ich etwas beachten muss.
Lacerbi
2
Ich habe ein paar weitere Details hinzugefügt: Bei der Implementierung der HPD-Uniform geht es darum, eine richtige kompakte Näherung für die HPD-Region zu finden. Die konvexe Hülle von Punkten mit hohen posterioren Werten ist (NP?) Schwer zu bestimmen, während sich Kugeln, die an diesen Punkten zentriert sind, schneiden können, was ein Problem der sekundären Normalisierungskonstante erzeugt.
Xi'an
2
@ Xi'an: sehr hilfreich, danke! Darf ich fragen: Welche der genannten Ansätze würden Sie derzeit empfehlen, wenn Sie nach einem allgemeinen Ansatz suchen, der in der Regel sofort funktioniert (dh keine Abstimmung / Überprüfung durch den Benutzer erforderlich ist)? Ich würde mich besonders für Modelle mit einer geringen (<50) Anzahl von Parametern, nicht normalen Seitenzähnen und starken Korrelationen zwischen Parametern interessieren.
Florian Hartig
1
Z