Obliczanie marginalnego prawdopodobieństwa na podstawie próbek MCMC

24

To jest powtarzające się pytanie (patrz ten post , ten post i ten post ), ale mam inny obrót.

Załóżmy, że mam kilka próbek z ogólnego próbnika MCMC. Dla każdej próbki znam wartość prawdopodobieństwa dziennika i dziennika przed . Jeśli to pomaga, znam również wartość prawdopodobieństwa dziennika na punkt danych, (ta informacja pomaga w przypadku niektórych metod, takich jak WAIC i PSIS-LOO).θlogf(x|θ)logf(θ)logf(xi|θ)

Chcę uzyskać (przybliżony) szacunek krańcowego prawdopodobieństwa, tylko na podstawie próbek, które mam, i ewentualnie kilku innych ocen funkcji (ale bez ponownego uruchamiania MCMC ad hoc ).

Przede wszystkim wyczyśćmy stół. Wszyscy wiemy, że estymator harmonicznych jest najgorszym estymatorem w historii . Przejdźmy dalej. Jeśli próbujesz Gibbsa z priory i posteriorami w formie zamkniętej, możesz użyć metody Chiba ; ale nie jestem pewien, jak uogólniać poza tymi przypadkami. Istnieją również metody, które wymagają zmodyfikowania procedury pobierania próbek (na przykład za pomocą temperowanych tylnych ), ale nie interesuje mnie to tutaj.

Podejście, o którym myślę, polega na aproksymacji rozkładu podstawowego parametrycznym (lub nieparametrycznym) kształcie , a następnie ustaleniu stałej normalizacyjnej jako problemu optymalizacji 1-D (tj. który minimalizuje pewien błąd między a , obliczone na próbkach). W najprostszym przypadku, załóżmy, że tył jest w przybliżeniu wielowymiarowy normalny, mogę dopasować jako wielowymiarową normalną i uzyskać coś podobnego do aproksymacji Laplace'a (mógłbym chcieć użyć kilku dodatkowych ocen funkcji w celu uściślenia pozycji tryb). Mógłbym jednak użyć jakog(θ)ZZZg(θ)f(x|θ)f(θ)g(θ)g(θ)bardziej elastyczna rodzina, taka jak wariacyjna mieszanka wielowymiarowych rozkładów t .

Rozumiem, że ta metoda działa tylko wtedy, gdy Zg(θ) jest rozsądnym przybliżeniem do f(x|θ)f(θ) , ale z jakiegokolwiek powodu lub przestrogi, dlaczego byłoby to bardzo nierozsądne Zrób to? Jakieś lektury, które poleciłbyś?

Podejście w pełni nieparametryczne wykorzystuje pewną rodzinę nieparametryczną, taką jak proces Gaussa (GP), do przybliżenia logf(x|θ)+logf(θ) (lub kilka innych jego nieliniowych transformacji, takich jak jako pierwiastek kwadratowy) i kwadraturę bayesowską, aby pośrednio zintegrować się nad podstawowym celem (patrz tutaj i tutaj ). To wydaje się być ciekawym alternatywnym podejściem, ale analogicznym w duchu (zauważ też, że lekarze ogólni byliby niewygodni w moim przypadku).

Lacerbi
źródło
6
Myślę, że Chib, S. i Jeliazkov, I. 2001 „Marginalna wiarygodność z Metropolis - wyjście Hastingsa uogólnia na normalne wyjścia MCMC - byłaby zainteresowana usłyszeniem doświadczeń z tym podejściem. Jeśli chodzi o lekarza rodzinnego - w zasadzie sprowadza się to do emulacji tylnej części ciała, którą można również rozważyć w przypadku innych problemów. Myślę, że problem polega na tym, że nigdy nie masz pewności co do jakości przybliżenia. Zastanawiam się również, czy próbka MCMC jest idealna do modelu GP, czy też należy zainwestować więcej w reszkę.
Florian Hartig,
2
(+1) Dzięki za referencję, wygląda na miejscu - sprawdzę to. Zgadzam się, że wszystkie podejścia oparte na modelach mogą być problematyczne (dobrą rzeczą w kwadraturze bayesowskiej jest to, że otrzymujesz oszacowanie niepewności, chociaż nie jestem pewien, jak jest skalibrowane). Na razie moim skromnym celem jest zrobienie czegoś, co jest „lepsze niż przybliżenie Laplace'a”.
lacerbi

Odpowiedzi:

26

Rozszerzenie autorstwa Chiba i Jeliazkowa (2001) niestety szybko staje się kosztowne lub bardzo zmienne, co powoduje, że nie jest często używane poza przypadkami próbkowania Gibbsa.

Chociaż istnieje wiele sposobów i podejść do problemu szacowania stałej normalizacji (jak ilustrują to dość różnorodne rozmowy podczas warsztatów Estimating Constant, które przeprowadziliśmy w zeszłym tygodniu na University of Warwick, dostępne tam slajdy ), niektóre rozwiązania wykorzystują bezpośrednio wyjście MCMC.Z

  1. Jak wspomniałeś, estymator średniej harmonicznej Newtona i Raftery'ego (1994) jest prawie niezmiennie słaby z powodu nieskończonej wariancji. Istnieją jednak sposoby na uniknięcie nieskończonej klątwy wariancji poprzez użycie zamiast tego skończonego celu wsparcia w harmonicznej średniej tożsamości , wybierając jako wskaźnik regionu HPD tylnej. Zapewnia to skończoną wariancję poprzez usunięcie ogonów ze średniej harmonicznej. (Szczegóły można znaleźć w artykule, który napisałem z Darrenem Wraithem oraz w rozdziale o normalizowaniu stałych napisanym przez Jean-Michela Marina.) W skrócie, metoda przetwarza dane wyjściowe MCMC aθ1,...,θMbetagatunku(θ)F(x|θ)aθ 0 I ρZ Z -1=

    α(θ)π(θ)f(x|θ)dπ(θ|x)=1Z
    αθ1,,θMprzez identyfikację (powiedzmy 20%) największych wartości celu i utworzenie jako jednolitego elementu nad połączeniem kulek wyśrodkowanych w tych symulacjach największej gęstości (HPD) oraz z promieniem , co oznacza oszacowanie stałej normalizującej jest podane przez βπ(θ)f(x|θ)αθi0ρZ
    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)
    jeśli jest wymiarem (poprawki dotyczą przecinających się kulek) i jeśli jest wystarczająco mały, aby kule nigdy się nie przecinały (co oznacza, że ​​w najlepszym razie tylko jeden wskaźnik na kule różnią się od zera). Wyjaśnienie mianownika jest takie, że jest to podwójna sumadθραM2βM2 warunki: z każdym terminem w integrującym się z .
    1βMi=1βM1Mm=1MU(θi0,ρ)(θm)same as with min×1π(θm)f(x|θm)
    θmZ1
  2. Innym podejściem jest przekształcenie stałej normalizującej w parametr. Brzmi to jak herezja statystyczna, ale artykuł Guttmanna i Hyvärinena (2012) przekonał mnie do czegoś przeciwnego. Bez wnikania w szczegóły, dobrym pomysłem jest obrócenie obserwowanego prawdopodobieństwa logarytmu do wspólnego prawdopodobieństwa dziennika który jest logarytmicznym prawdopodobieństwem procesu punktu Poissona z funkcją intensywności Z

    i=1nf(xi|θ)nlogexpf(x|θ)dx
    i=1n[f(xi|θ)+ν]nexp[f(x|θ)+ν]dx
    exp{f(x|θ)+ν+logn}
    Jest to model alternatywny, ponieważ pierwotne prawdopodobieństwo nie pojawia się jako margines powyższego. Tylko tryby się pokrywają, a tryb warunkowy w ν zapewnia stałą normalizującą. W praktyce powyższe prawdopodobieństwo procesu Poissona jest niedostępne, a Guttmann i Hyvärinen (2012) oferują przybliżenie za pomocą regresji logistycznej. Aby jeszcze lepiej połączyć się z pytaniem, szacunek Geyera jest MLE, stąd rozwiązanie problemu maksymalizacji.
  3. Podejście powiązane to podejście regresji logistycznej Charliego Geyera . Podstawowym pojęciem jest dodanie do próbki MCMC z innej próbki od znanego celu, np. Najlepsze zgadywanie w , i uruchomienie regresja logistyczna indeksu rozkładu za danymi (1 dla i 0 dla ). Regresory są wartościami obu gęstości, znormalizowanymi lub nie. Zdarza się to bezpośrednio związane z próbkowaniem pomostowym Gelmana i Menga (1997), które również przetwarza próbki z różnych celów. I późniejsze wersje, jak MLE Menga.π(θ|x)π(θ|x)g(θ)π(θ|x)g(θ)
  4. Innym podejściem, które zmusza do uruchomienia określonego próbnika MCMC, jest próbkowanie zagnieżdżone przez Skilling . Chociaż ja [i inni] mam pewne zastrzeżenia co do wydajności metody, jest ona dość popularna w astrostatyce i kosmologii, z oprogramowaniem dostępnym jak multinest .
  5. Ostatnim [potencjalnym, jeśli nie zawsze możliwym] rozwiązaniem jest wykorzystanie przedstawionej przez Savage'a-Dickeya reprezentacji czynnika Bayesa w przypadku osadzonej hipotezy zerowej. Jeśli null zapisuje jako o interesującym parametrze, a jeśli jest pozostałą [uciążliwą] częścią parametru modelu, przyjmując wcześniejszą postać , współczynnik Bayesa porównaniu z alternatywnymi zapisuje jako gdzie oznacza marginalną tylną gęstość o określonej wartościH0:θ=θ0ξπ1(θ)π2(ξ)H0
    B01(x)=πθ(θ0|x)π1(θ0)
    πθ(θ0|x)θθ0. W przypadku gdy marginalna gęstość poniżej wartości zerowej jest dostępny w zamknięta forma, można uzyskać gęstość krańcową dla nieograniczonego modelu od współczynnika Bayesa. (Ta reprezentacja Savage'a-Dickeya opiera się na konkretnych wersjach trzech różnych gęstości, a więc jest obarczona niebezpieczeństwem, nie wspominając nawet o obliczeniowym wyzwaniu polegającym na wytworzeniu marginalnego tylnego).H0:θ=θ0
    m0(x)=Ξf(x|θ0,ξ)π2(ξ)dξ
    ma(x)=Θ×Ξf(x|θ,ξ)π1(θ)π2(ξ)dθdξ

[Oto zestaw slajdów, które napisałem o szacowaniu stałych normalizujących na warsztatach NIPS w grudniu.]

Xi'an
źródło
2
(+1) Niezwykle bogata odpowiedź, dziękuję. Będzie to dla mnie przydatne i, jak sądzę, wielu innych ludzi. Zajmie mi trochę czasu, aby przyjrzeć się różnym podejściom, a potem mogę wrócić z konkretnymi pytaniami.
lacerbi 30.04.16
2
Począwszy od punktu (1) ... Czytam odpowiednie artykuły. „Skorygowany” estymator średniej harmonicznej wydaje się dokładnie tym , czego szukałem. Jest czysty i łatwy do obliczenia, biorąc pod uwagę wyjście MCMC. Więc ... jaki jest haczyk? Nie wygląda na to, że metoda jest szeroko stosowana, sądząc po szybkim wyszukiwaniu w Google Scholar. Jakie są jego ograniczenia? (oprócz potrzeby zidentyfikowania regionów HPD, które, jak sądzę, mogą stać się problemem dla bardzo skomplikowanych bocznych w dużym wymiarze). Na pewno spróbuję - ale zastanawiam się, czy jest coś, o co muszę się uważać.
lacerbi
2
Dodałem jeszcze kilka szczegółów: problem z implementacją munduru HPD polega na znalezieniu odpowiedniego kompaktowego przybliżenia dla regionu HPD. Wypukły kadłub punktów o wysokich wartościach tylnych jest (NP?) Trudny do ustalenia, podczas gdy kulki wycentrowane w tych punktach mogą się przecinać, co stwarza wtórny normalizujący stały problem.
Xi'an
2
@ Xi'an: bardzo pomocny, dzięki! Czy mogę zapytać: ze wszystkich wymienionych podejść, jakie byłoby obecnie twoje zalecenie, jeśli ktoś szuka ogólnego podejścia, które zwykle działa od razu po wyjęciu z pudełka (tj. Nie wymaga strojenia / sprawdzania od użytkownika)? Byłbym szczególnie zainteresowany w przypadku modeli o niskiej (<50) liczbie parametrów, nietypowych tylnych ścianach i silnych korelacjach między parametrami.
Florian Hartig
1
@FlorianHartig: fakt, że ogólne oprogramowanie, takie jak BŁĘDY, nie zwraca ogólnego oszacowania pewnym sensie ujawnia zakres problemu. Wiele rozwiązań, które można znaleźć w literaturze specjalistycznej, nie przyniosło konsensusu. Dlatego moim zaleceniem byłoby wybranie rozwiązania regresji logistycznej Geyera, które jest nieco niewrażliwe na wymiary. Z
Xi'an