Obliczanie marginalnego prawdopodobieństwa na podstawie próbek MCMC

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).$\theta$$\log f(\textbf{x} | \theta)$$\log f(\theta)$$\log f(x_i | \theta)$

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ć jako$g(\theta)$$Z$$Z$$Z g(\theta)$$f(\textbf{x}|\theta) f(\theta)$$g(\theta)$$g(\theta)$bardziej elastyczna rodzina, taka jak wariacyjna mieszanka wielowymiarowych rozkładów $t$ .

Rozumiem, że ta metoda działa tylko wtedy, gdy $Z g(\theta)$ jest rozsądnym przybliżeniem do $f(\textbf{x}|\theta) f(\theta)$ , 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 $\log f(\textbf{x}|\theta) + \log f(\theta)$ (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).

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.$\mathfrak{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

   $$\int \dfrac{\alpha(\theta)}{\pi(\theta)f(x|\theta)}\text{d}\pi(\theta|x)=\frac{1}{\mathfrak{Z}}$$ $\alpha$$\theta_1,\ldots,\theta_M$przez 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 $\beta$$\pi(\theta)f(x|\theta)$$\alpha$$\theta^0_i$$\rho$$\mathfrak{Z}$ $$\hat{\mathfrak{Z}}^{-1}=\underbrace{\frac{1}{\beta M^2}\sum_{m=1}^M}_{\text{double sum over}\\\beta M\text{ ball centres }\theta_i^0\\\text{and $M$ simulations } \theta_m} \underbrace{\mathbb{I}_{(0,\rho)}(\min_i||\theta_m-\theta^0_i||)\{\pi(\theta_m)f(x|\theta_m)\}^{-1}\big/\overbrace{\pi^{d/2}\rho^d\Gamma(d/2+1)^{-1}}^{\text{volume of ball with radius $\rho$}}}_{\dfrac{\beta M\alpha(\theta_m)}{\pi(\theta_m)f(x|\theta_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 suma$d$$\theta$$\rho$$\alpha M^2$$\beta M^2$ warunki: z każdym terminem w integrującym się z . $$\frac{1}{\beta M}\sum_{i=1}^{\beta M} \underbrace{\frac{1}{M}\sum_{m=1}^M {\cal U}(\theta_i^0,\rho)(\theta_m)}_{\text{same as with $\min$}} \times \frac{1}{\pi(\theta_m)f(x|\theta_m)}$$ $\theta_m$${\mathfrak{Z}}^{-1}$

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 $\mathfrak{Z}$

   $$\sum_{i=1}^n f(x_i|\theta) - n \log \int \exp f(x|\theta) \text{d}x$$ $$\sum_{i=1}^n[f(x_i|\theta)+\nu]-n\int\exp[f(x|\theta)+\nu]\text{d}x$$ $$\exp\{ f(x|\theta) + \nu +\log n\}$$ 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.$\pi(\theta|x)$$\pi(\theta|x)$$g(\theta)$$\pi(\theta|x)$$g(\theta)$
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ści$H_0: \theta=\theta_0$$\xi$$\pi_1(\theta)\pi_2(\xi)$$H_0$ $$\mathfrak{B}_{01}(x)=\dfrac{\pi^\theta(\theta_0|x)}{\pi_1(\theta_0)}$$ $\pi^\theta(\theta_0|x)$$\theta$$\theta_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).$H_0: \theta=\theta_0$ $$m_0(x)=\int_\Xi f(x|\theta_0,\xi)\pi_2(\xi)\text{d}\xi$$ $$m_a(x)=\int_{\Theta\times\Xi} f(x|\theta,\xi)\pi_1(\theta)\pi_2(\xi)\text{d}\theta\text{d}\xi$$

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