Zagadnienia

12.1 Statystyka bayesowska
- 12.1.1 Hierarchiczny model klasyfikacji
- 12.1.2 Próbnik Gibbsa w modelu hierarchicznym
12.2 Estymatory największej wiarogodności
- 12.2.1 Model auto-logistyczny
  - 12.2.1.1 Próbnik Gibbsa w modelu auto-logistycznym
  - 12.2.1.2 Aproksymacja wiarogodności

12. Markowowskie Monte Carlo III. Przykłady zastosowań

12.1. Statystyka bayesowska

Algorytmy MCMC zrewolucjonizowały statystykę bayesowską. Stworzyły możliwość obliczania (w przybliżeniu) rozkładów a posteriori w sytuacji, gdy dokładne, analityczne wyrażenia są niedostępne. W ten sposób statystycy uwolnili się od koniecznośći używania nadmiernie uproszczonych modeli. Zaczęli śmiało budować modele coraz bardziej realistyczne, zwykle o strukturze hierarchicznej. Przedstawię to na jednym dość typowym przykładzie, opartym na pracy [21]. Inne przykłady i doskonały wstęp do tematyki zastosowań MCMC można znaleźć w pracy Geyera [7].

12.1.1. Hierarchiczny model klasyfikacji

Przykład 12.1 (Statystyka małych obszarów)

Zacznijmy od opisania problemu tak zwanych ,,małych obszarów”, który jest dość ważny w dziedzinie badań reprezentacyjnych, czyli w tak zwanej ,,statystyce oficjalnej”. Wyobraźmy sobie, że w celu zbadania kondycji przedsiębiorstw losuje się próbkę, która liczy (powiedzmy) 3500 przedsiębiorstw z całego kraju. Na podstawie tej losowej próbki można dość wiarygodnie szacować (estymować) pewne parametry opisujące populację przedsiębiorstw w kraju. Czy jednak można z rozsądną dokłdnością oszacować sprzedaż w powiecie garwolińskim? W Polsce mamy ponad 350 powiatów. Na jeden powiat przypada średnio 10 przedsiębiorstw wybranych do próbki. Małe obszary to pod-populacje w których rozmiar próbki nie jest wystarczający, aby zastosować ,,zwykłe” estymatory (średnie z próbki). Podejście bayesowskie pozwala ,,pożyczać informację” z innych obszarów. Zakłada się, że z każdym małym obszarem związany jest nieznany parametr, który staramy się estymować. Obserwacje pochodzące z określonego obszaru mają rozkład prawdopodobieństwa zależny od odpowiadającego temu obszarowi parameru. Parametry, zgodnie z filozofią bayesowską, traktuje się jak zmienne losowe. W najprostszej wersji taki model jest zbudowany w następujący sposób.

12.1.1.1. Model bayesowski

$y_{{ij}}\sim{\rm N}(\theta _{{i}},\sigma^{2})$ – badana cecha dla $j$ -tej wylosowanej jednostki $i$ -tego obszaru, $(j=1,\ldots,n_{i})$ , $(i=1,\ldots,k)$ ,
$\theta _{{i}}\sim{\rm N}(\mu,\upsilon^{2})$ – interesująca nas średnia w $i$ -tym obszarze,
$\mu$ – średnia w całej populacji.

Ciekawe, że ten sam model pojawia się w różnych innych zastosowanianich, na przykład w matematyce ubezpieczeniowej. Przytoczymy klasyczny rezultat dotyczący tego modelu, aby wyjaśnić na czym polega wspomniane ,,pożyczanie informacji”.

12.1.1.2. Estymator bayesowski

W modelu przedstawionym powyżej, łatwo obliczyć estymator bayesowski (przy kwadratowej funkcji straty), czyli wartość oczekiwaną a posteriori. Następujący wzór jest bardzo dobrze znany specjalistom od małych obszarów i aktuariuszom.

$\hat{\theta}_{{i}}=\mathbb{E}(\theta _{i}|y)=z_{{i}}\bar{y}_{{i}}+(1-z_{i})\mu,\qquad z_{i}=\frac{n_{i}\upsilon^{2}}{n_{i}\upsilon^{2}+\sigma^{2}}.$

Estymator bayesowski dla $i$ -tego obszaru jest średnią ważoną $\bar{y}$ (estymatora opartego na danych z tego obszaru) i wielkości $\mu$ , która opisuje całą populację, a nie tylko $i$ -ty obszar. Niestety, proste estymator napisany powyżej zależy od parametrów $\mu$ , $\sigma$ i $v$ , które w praktyce są nieznane i które trzeba estymować. Konsekwentnie bayesowskie podejście polega na traktowaniu również tych parametrów jako zmiennych losowych, czyli nałożeniu na nie rozkłądów a priori. Powstaje w ten sposób model hierarchiczny.

12.1.1.3. Hierarchiczny model bayesowski

Uzupełnijmy rozpatrywany powyżej model, dobudowując ,,wyższe piętra” hierarchii. potraktujemy mianowicie parametry rozkładów a priori: $\mu$ , $\sigma$ i $v$ jako zmienne losowe i wyspecyfikujemy ich rozkłady a priori.

$y_{{ij}}\sim{\rm N}(\theta _{{i}},\sigma^{2})$ ,
$\theta _{{i}}\sim{\rm N}(\mu,\upsilon^{2})$ ,
$\mu\sim{\rm N}(m,\tau^{2})$ ,
$\sigma^{{-2}}\sim{\rm Gamma}(p,\lambda)$ ,
$\upsilon^{{-2}}\sim{\rm Gamma}(q,\kappa)$ .

Zakładamy przy tym, że $\mu$ , $\sigma$ i $v$ są a priori niezależne (niestety, są one zależne a posteriori). Na szczycie hierarchii mamy ,,hiperparametry” $m$ , $\tau$ , $p$ , $\lambda$ , $q$ , $\kappa$ , o których musimy założyć, ze są znanymi liczbami.

Łączny rozkład prawdopodobieństwa wszystkich zmiennych losowych w modelu ma postać

$p(y,\theta,\mu,\sigma^{{-2}},\upsilon^{{-2}})=p(y|\theta,\sigma^{{-2}})p(\theta|\mu,\upsilon^{{-2}})p(\mu)p(\sigma^{{-2}})p(\upsilon^{{-2}}).$

We wzorze powyżej i w dalej traktujemy (trochę nieformalnie) $\sigma^{{-2}}$ i $\upsilon^{{-2}}$ jako pojedyncze symbole nowych zmiennych, żeby nie mnożyć oznaczeń. Rozkład prawdopodobieństwa a posteriori jest więc taki:

$p(\theta,\mu,\sigma^{{-2}},\upsilon^{{-2}}|y)=\frac{p(y,\theta,\mu,\sigma^{{-2}},\upsilon^{{-2}})}{p(y)}.$

To jest rozkład ,,docelowy” $\pi$ , na przestrzeni ${\cal X}=\mathbb{R}^{{k+3}}$ , ze nieznaną stałą normującą $1/p(y)$ . Choć wygląda na papierze dość prosto, ale obliczenie rozkładów brzegowych, wartości oczekiwanych i innych charakterystyk jest, łagodnie mówiąc, trudne.

12.1.2. Próbnik Gibbsa w modelu hierarchicznym

Rozkłady warunkowe poszczególnych współrzędnych są proste i łatwe do generowania. Można te rozkłady ,,odczytać” uważnie patrząc na rozkład łączny:

$\begin{split} p(\theta,\mu,{\upsilon^{{-2}}},\sigma^{{-2}}|y)&\propto(\sigma^{{-2}})^{{n/2}}\exp\left\{-\frac{\sigma^{{-2}}}{2}\sum _{{i=1}}^{k}\sum _{{j=1}}^{{n_{i}}}(y_{{ij}}-\theta _{i})^{2}\right\}\\ &\cdot{(\upsilon^{{-2}})^{{k/2}}\exp\left\{-\frac{\upsilon^{{-2}}}{2}\sum _{{i=1}}^{k}(\theta _{i}-\mu)^{2}\right\}}\\ &\cdot\exp\left\{-\frac{\tau^{{-2}}}{2}(\mu-m)^{2}\right\}\\ &\cdot(\sigma^{{-2}})^{{q-1}}\exp\{-\kappa\sigma^{{-2}}\}\\ &\cdot{(\upsilon^{{-2}})^{{p-1}}\exp\{-\lambda\upsilon^{{-2}}\}}.\end{split}$

Dla ustalenia uwagi zajmijmy się rozkładem warunkowym zmiennej $\upsilon^{{-2}}$ . Kolorem niebieskim oznaczyliśmy te czynniki łącznej gęstości, które zawierają $\upsilon^{{-2}}$ . Pozostałe, czarne czynniki traktujemy jako stałe. Stąd widać, jak wygląda rozkład warunkowy $\upsilon^{{-2}}$ ,, przynajmniej z dokładnością do proporcjonalności:

$\begin{split} p({\upsilon^{{-2}}}|y,\theta,\mu,\sigma^{{-2}})&\propto{(\upsilon^{{-2}})^{{k/2+p-1}}}\\ &{\cdot\exp\left\{-\left(\frac{1}{2}\sum _{{i=1}}^{k}(\theta _{i}-\mu)^{2}+\lambda\right)\upsilon^{{-2}}\right\}}.\end{split}$

Jest to zatem rozkład ${\rm Gamma}(k/2+p,\sum _{{i=1}}^{k}(\theta _{i}-\mu)^{2}/2+\lambda)$ . Zupełnie podobnie rozpoznajemy inne (pełne) rozkłady warunkowe:

$\begin{split}\upsilon^{{-2}}|y,\theta,\mu,\sigma^{{-2}}\frac{}{}\quad&\sim{\rm Gamma}\left(\frac{k}{2}+p,\frac{1}{2}\sum _{{i=1}}^{k}(\theta _{i}-\mu)^{2}+\lambda\right),\\ \sigma^{{-2}}|y,\theta,\mu,\upsilon^{{-2}}\quad&\sim{\rm Gamma}\left(\frac{n}{2}+q,\frac{1}{2}\sum _{{i=1}}^{k}\sum _{{j=1}}^{{n_{i}}}(y_{{ij}}-\theta _{i})^{2}+\kappa\right),\\ \mu|y,\theta,\sigma^{{-2}},\upsilon^{{-2}}\quad&\sim{\rm N}\left(\frac{k\tau^{{2}}}{k\tau^{{2}}+\upsilon^{{2}}}\bar{\theta}+\frac{\upsilon^{{2}}}{k\tau^{{2}}+\upsilon^{{2}}}m,\frac{\tau^{{2}}\upsilon^{{2}}}{k\tau^{{2}}+\upsilon^{{2}}}\right),\\ \theta _{i}|y,\theta _{{-i}},\mu,\sigma^{{-2}},\upsilon^{{-2}}\quad&\sim{\rm N}\left(\frac{n\upsilon^{{2}}}{n\upsilon^{{2}}+\sigma^{{2}}}\bar{y}_{i}+\frac{\sigma^{{2}}}{n\upsilon^{{2}}+\sigma^{{2}}}\mu,\frac{\upsilon^{{2}}\sigma^{{2}}}{n\upsilon^{{2}}+\sigma^{{2}}}\right),\end{split}$

(12.1)

gdzie, rzecz jasna, $n=\sum n_{i}$ , $\bar{\theta}=\sum _{i}\theta _{i}/k$ i $\theta _{{-i}}=(\theta _{k})_{{k\not=i}}$ . Zwróćmy uwagę, że współrzędne wektora $\theta$ są warunkowo niezależne (pełny rozkład warunkowy $\theta _{i}$ nie zależy od $\theta _{{-i}}$ ). Dzięki temy możemy w próbniku Gibbsa potraktować $\theta$ jako cały ,,blok” współrządnych i zmieniać ,,na raz”.

Próbnik Gibbsa ma w tym modelu przestrzeń stanów ${\cal X}$ składającą się z punktów $x=(\theta,\mu,\sigma^{{-2}},\upsilon^{{-2}})\in\mathbb{R}^{{k+3}}$ . Reguła przejścia próbnika w wersji systematycznej (duży krok ,,SystemPG”),

$\underbrace{\Big(\theta,\mu,\sigma^{{-2}},\upsilon^{{-2}}\Big)}_{{X_{t}}}\longmapsto{\underbrace{\Big(\theta,\mu,\sigma^{{-2}},\upsilon^{{-2}}\Big)}_{{X_{{t+1}}}},}$

jest złożona z następujących ,,małych kroków”:

Wylosuj ${\upsilon^{{-2}}}\sim p(\upsilon^{{-2}}|y,\theta,\mu,\sigma^{{-2}}\phantom{,\upsilon^{{-2}}}){\quad=\quad{\rm Gamma}(\ldots)}$ ,
Wylosuj ${\sigma^{{-2}}}\sim p(\sigma^{{-2}}|y,\theta,\mu,\phantom{\sigma^{{-2}},}{\upsilon^{{-2}}}){\quad=\quad{\rm Gamma}(\ldots)}$ ,
Wylosuj ${\mu}\phantom{{}^{{-2}}}\sim p(\mu\phantom{{}^{{-2}}}|y,\theta,\phantom{\mu,}{\sigma^{{-2}},\upsilon^{{-2}}}){\quad=\quad{\rm N}(\ldots)}$ ,
Wylosuj ${\theta}\phantom{{}^{{-2}}}\sim p(\theta\phantom{{}^{{-2}}}|y,\phantom{\theta,}{\mu,\sigma^{{-2}},\upsilon^{{-2}}}){\quad=\quad{\rm N}(\ldots)}$ .

Łańcuch Markowa jest zbieżny do rozkładu a posteriori:

$X_{t}\to\pi(\cdot)=p(\theta,\mu,\sigma^{{-2}},\upsilon^{{-2}}|y).$

Najbardziej interesujące są w tym hierarchicznym modelu zmienne $\theta _{i}$ (pozostałe zmienne można uznać za ,,parametry zkłócające). Dla ustalenia uwagi zajmijmy się zmienną $\theta _{1}$ (powiedzmy, wartością średnią w pierwszym małym obszarze). Estymator bayesowski jest to wartość oczekiwana a posteriori tej zmiennej:

$\mathbb{E}(\theta _{1}|y)=\int\cdots\int\theta _{1}p(\theta,\mu,\sigma^{{-2}},\upsilon^{{-2}}|y){\rm d}\theta _{2}\cdots{\rm d}\theta _{k}{\rm d}\mu{\rm d}\sigma^{{-2}}{\rm d}\upsilon^{{-2}}.$

Aproksymacją MCMC interesującej nas wielkości są średnie wzdłuż trajektorii łańcucha:

$\theta _{1}(X_{0}),\theta _{1}(X_{1}),\ldots,\theta _{1}(X_{t}),\ldots$

gdzie $\theta _{1}(x)=\theta _{1}$ dla $x=(\theta _{1},\ldots,\theta _{k},\mu,\sigma^{{-2}},\upsilon^{{-2}})$ .

$Trajektorie zmiennej $v^{{-2}}$ dla dwóch punktów startowych$

Rys. 12.1. Trajektorie zmiennej $v^{{-2}}$ dla dwóch punktów startowych.

$Skumulowane średnie zmiennej $v^{{-2}}$ dla dwóch punktów startowych$

Rys. 12.2. Skumulowane średnie zmiennej $v^{{-2}}$ dla dwóch punktów startowych.

Na Rysunku 12.1 pokazane są dwie przykładowe trajektorie współrzędnej $\upsilon^{{-2}}$ dla PG poruszającego się po przestrzeni $k+3=1003$ wymiarowej (model uwzględniający 1000 małych obszarów). Dwie trajektorie odpowiadają dwum różnym punktom startowym. Dla innych zmiennych rysunki wyglądają bardzo podobnie. Uderzające jest to, jak szybko trajektoria zdaje się ,,osiągać” rozkład stacjonarny, przynajmniej wizualnie. Na Rysunku 12.2 pokazane są kolejne ,,skumulowane” średnie dla tych samych dwóch trajektorii zmiennej $\upsilon^{{-2}}$ .

12.2. Estymatory największej wiarogodności

Metody MCMC w połączeniu z ideą losowania istotnego (Rozdział 8) znajdują zastosowanie również w statystyce ,,częstościowej”, czyli nie-bayesowskiej. Pokażę przykład, w którym oblicza się meodami Monte carlo estymator największej wiarogodności.

12.2.1. Model auto-logistyczny

Niech $x=(x_{1},\ldots,x_{d})$ będzie wektorem (konfiguracją) binarnych zmiennych losowych na przestrzeni ${\cal X}=\{ 0,1\}^{d}$ . Rozważmy następujący rozkład Gibbsa:

$p_{{\theta}}(x)=\frac{1}{Z(\theta)}\,\exp\left\{\sum _{{i,j=1}}^{{d}}\theta _{{ij}}x_{{i}}x_{{j}}\right\}.$

Rolę parametru gra macierz $(\theta _{{ij}})$ . Zakłąd się, że jest to maciezrz symetryczna . Stała normująca $Z(\theta)=\sum _{{x\in{\cal X}}}\exp\left\{\sum _{{i,j=1}}^{{d}}\theta _{{ij}}x_{{i}}x_{{j}}\right\}$ jest typowo (dla dużych $d$ ) niemożliwa do obliczenia. Stanowi to, jak dalej zobaczymy, poważny problem dla statystyków.

Przykład 12.2 (Statystyka przestrzenna)

W zastosowaniach ,,przestrzennych” indeks $i\in\{ 1,\ldots,d\}$ interpretuje się jako ,,miejsce”. Zbiór miejsc wyposażony jest w strukturę grafu. Krawędzie łączą miejsca ,,sąsiadujące”. Piszemy $i\sim j$ . Tego typu modele mogą opisywać na przykład rozprzestrzenianie się chorób lub występowanie pewnych gatunków. Wartośc $x_{i}=1$ oznacza obecność gatunku lub występowanie choroby w miejscu $i$ . Najprostszy model zakłada, że każda zmienna $x_{i}$ zależy tylko od swoich ,,sąsiadów” i to w podobny sposób w całym rozpatrywanym obszarze. W takim modelu mamy tylko dwa parametry $\theta=(\theta _{0},\theta _{1})$ :

$\theta _{{ij}}=\begin{cases}0&i\not\sim j,i\not=j;\\ \theta _{1}&i\sim j;\\ \theta _{0}&i=j.\\ \end{cases}$

Parametr $th_{0}$ opisuje ,,skłonność” pojedynczej zmiennej do przyjmowania wartości 1, zaś parametr $\theta _{1}$ odpowiada za zależność od zmiennych sąsiadujących (zakaźność choroby, powiedzmy). W typowej dla statystyki przestrzennej sytuacji, rozpatruje się nawet dziesiątki tysięcy ,,miejsc”. Stała $Z(\theta)$ jest wtedy sumą niewyobrażalnie wielu (dokładnie $2^{d}$ ) składników.

12.2.1.1. Próbnik Gibbsa w modelu auto-logistycznym

,,Pełne” rozkłady warunkowe (full conditionals) są w modelu auto-logistycznym bardzo proste:

$p_{{\theta}}(x_{{i}}=1\mid x_{{-i}})=\frac{\exp\left(\theta _{{ii}}+\sum\limits _{{j=1\atop j\neq i}}^{{d}}\theta _{{ij}}x_{{j}}\right)}{1+\exp\left(\theta _{{ii}}+\sum\limits _{{j=1\atop j\neq i}}^{{d}}\theta _{{ij}}x_{{j}}\right)},$

gdzie $x_{{-i}}=(x_{{j}},j\neq i)$ . Zatem:

Symulowanie $x\sim p_{\theta}$ jest łatwe za pomocą próbnika Gibbsa (PG):
$x_{1}\sim p_{\theta}(x_{1}|x_{{-1}})$ ,
$x_{2}\sim p_{\theta}(x_{2}|x_{{-2}})$ , $\ldots$

12.2.1.2. Aproksymacja wiarogodności

Estymator największej wiarogodności obliczany metodą Monte Carlo został zaproponowany w pracy Geyer and Thopmpson (1992, JRSS). Rozważmy bardziej ogólną wykładniczą rodzinę rozkładów prawdopodobieństwa:

$p_{{\theta}}(x)=\frac{1}{Z(\theta)}\exp\left[\theta^{\top}T(x)\right],$

gdzie $T(x)$ jest wektorem statystyk dostatecznych. Rozkłady autologistyczne tworzą rodzinę wykładniczą. Wystarczy $\theta$ ustawić w wektor, statystykami $T(x)$ są $x_{i}$ i $x_{i}x_{j}$ . Logarytm wiarogodności $L(\theta)=\log p_{\theta}(x)$ jest dany pozornie prostym wzorem:

$L(\theta)=\theta^{\top}T(x)-\log Z(\theta).$

Kłopot jest ze stałą normującą $Z(\theta)$ , która w bardziej skomplikowanych modelach jest trudna lub wręcz niemożliwa do obliczenia. Wspomnieliśmy, że tak właśnie jest dla typowych modeli autologistycznych. Istnieje jednak prosty sposób aproksymacji stałej $Z(\theta)$ metodą MC. Istotnie, wybierzmy (w zasadzie dowolny) punkt $\theta _{*}$ i zauważmy, że

$Z(\theta)=\sum _{x}\exp\left[\theta^{\top}T(x)\right]\\ =\sum _{x}\exp\left[(\theta-\theta _{*})^{\top}T(x)\right]\exp\left[\theta _{*}^{\top}T(x)\right].$

Pozwala to wyrazić wielkość $Z(\theta)/{Z(\theta _{*})}$ jako wartość oczekiwaną względem rozkładu $p_{{\theta _{{*}}}}$ :

$\begin{split} Z(\theta)/{Z(\theta _{*})}&=\sum _{x}\exp\left[(\theta-\theta _{*})^{\top}T(x)\right]p_{{\theta _{*}}}(x)\\ &=\mathbb{E}_{{\theta _{*}}}\exp\left[(\theta-\theta _{*})^{\top}T(x)\right].\end{split}$

Powyższe rozważania nie są niczym nowym: to jest po prostu przedstawiony w Rozdziale 8 schemat losowania istotnego, w specjalnym przypadku rodziny wykładniczej. Jeśli celem jest maksymalizacja wiarogodności $Z(\theta)$ , to nieznana stała $Z(\theta _{*})$ zupełnie nie przeszkadza (bo interesującą nas funkcję wiarygodności wystarczy znać z dokładnością do stałej).

Jeśli umiemy generować zmienne losowe o rozkładzie $p_{{\theta _{{*}}}}$ , to sposób przybliżonego obliczania $\hat{Z}(\theta)$ jest w zasadzie oczywisty. Generujemy próbkę MC: $x_{{*}}(k)\sim p_{{\theta _{{*}}}}$ , $k=1,\ldots,n_{{*}}$ gdzie $\theta _{*}$ jest w zasadzie dowolne, zaś $n_{{*}}$ jest możliwie największe. Obliczamy

$\hat{Z}(\theta)/Z(\theta _{*})=\frac{1}{n_{*}}\sum _{{k=1}}^{{n_{*}}}\exp\left[(\theta-\theta _{{*}})^{\top}T(x_{*}(k))\right].$

Aproksymacja logarytmu wiarogodności polega na wstawieniu $\hat{Z}(\theta)$ w miejsce $\hat{Z}(\theta)$ we wzorze na $L(\theta)=\log p_{\theta}(x)$ :

$\hat{L}_{{\rm MC}}(\theta)=\theta^{\top}T(x)-\log\underbrace{\sum _{{k=1}}^{{n_{*}}}{\exp[(\theta-\theta _{{*}})^{\top}T(x_{*}(k))}]}_{{\text{ przybliżenie MC }Z(\theta)}}+{\rm const}.$

Pozostaje wiele szczegółów do dopracowania. Umiemy obliczać wiarogodność jako funkcję $\theta$ , ale trzeba jeszcze znależć maksimum tej funkcji. Dobór $\theta _{*}$ nie wpływa na poprawność rozumowania ale może mieć zasadniczy wpływ na efektywność algorytmu. Nie będziemy się w to zagłębiać. Wspomnimy tylko, jak omawiana metoda jest związana z markowowskimi metodami Monte Carlo, MCMC. Wróćmy do rozkładu auto-logistycznego. Jak losować próbkę z tego rozkładu? Tu przychodzi z pomocą próbnik Gibbsa. W rzeczy samej, mamy do czynienia z kolejną aproksymacją: PG jest tylko przybliżonym algorytmem generowania z rozkładu $p_{{\theta _{{*}}}}$ , które to losowanie pozwala obliczać wiarogodność w przybliżeniu. Momo wszystko, metody MCMC oferują pewne wyjście lepsze niż bezradne rozłożenie rąk.

Treść automatycznie generowana z plików źródłowych LaTeXa za pomocą oprogramowania wykorzystującego LaTeXML.

Symulacje stochastyczne i metody Monte Carlo