Zagadnienia

11.1 Odwracalność
11.2 Algorytm Metropolisa-Hastingsa
11.3 Próbnik Gibbsa

11. Markowowskie Monte Carlo II. Podstawowe algorytmy

11.1. Odwracalność

Najważniejsze algorytmy MCMC są oparte na idei odwracalności łańcucha Markowa.

Definicja 11.1

Łańcuch o jądrze $P$ jest odwracalny względem rozkładu prawdopodobieństwa $\pi$ , jeśli dla dowolnych $A,B\subseteq{\cal X}$ mamy

$\int _{A}\pi({\rm d}x)P(x,B)=\int _{B}\pi({\rm d}y)P(y,A).$

W skrócie,

$\qquad\pi({\rm d}x)P(x,{\rm d}y)=\pi({\rm d}y)P(y,{\rm d}x).$

Odwracalność implikuje, że rozkład $\pi$ jest stacjonarny. jest to dlatego ważne, że sprawdzanie odwracalności jest stosunkowo łatwe.

Twierdzenie 11.1

Jeśli $\pi({\rm d}x)P(x,{\rm d}y)=\pi({\rm d}y)P(y,{\rm d}x)$ to $\pi P=P$ .

Dowód

$\pi P(B)=\int _{{\cal X}}\pi({\rm d}x)P(x,B)=\int _{B}\pi({\rm d}y)P(y,{\cal X})=\int _{B}\pi({\rm d}y)=\pi(B)$ .

∎

11.2. Algorytm Metropolisa-Hastingsa

To jest pierwszy historycznie i wciąż najważniejszy algorytm MCMC. Zakładamy, że umiemy generować łańcuch Markowa z pewnym jądrem $q$ . Pomysł Metropolisa polega na tym, żeby zmodyfikować ten łańcuch wprowadzając specjalnie dobraną regułę akceptacji w taki sposób, żeby wymusić zbieżność do zadanego rozkładu $\pi$ . W dalszym ciągu systematycznie utożsamiamy rozkłady prawdopodobieństwa z ich gęstościami, aby nie mnożyć oznaczeń. Mamy zatem:

Rozkład docelowy: $\pi({\rm d}x)=\pi(x){\rm d}x$ .
Rozkład ,,propozycji”: $q(x,{\rm d}y)=q(x,y){\rm d}y$ .
Reguła akceptacji:

$a(x,y)=\frac{\pi(y)q(y,x)}{\pi(x)q(x,y)}\wedge 1.$

Algorytm Metropolisa-Hastingsa (MH) interpretujemy jako ,,błądzenie losowe” zgodnie z jądrem przejścia $q$ , zmodyfikowane poprzez odrzucanie niektórych ruchów, przy czym reguła akceptacji/odrzucania zależy w specjalny sposób od $\pi$ . Pojedynczy krok algorytmu jest następujący.

function $KrokMH(x)$

Gen $y\sim q(x,\cdot)$ ; { propozycja }

Gen $U\sim{\rm U}(0,1)$

if $U>a(x,y)$ } then

begin

$y:=x$ ; { ruch odrzucony z pr-stwem $1-a(x,y)$ }

end

$KrokMH:=y$

Graficznie to można przedstawić w takiej postaci.

$\xymatrix{&X_{n}=x\ar[d]&\\ &\ar[dl]_{{a(x,y)}}y\sim q(x,\cdot)\ar[dr]^{{1-a(x,y)}}&\\ X_{{n+1}}=y&&X_{{n+1}}=x\\ }$

Łańcuch Markowa powstaje zgodnie z następującym schematem:

Gen $X_{{0}}\sim\pi _{{0}}$ ; { start }

for $n:=1$ to $\infty$

begin

$X_{{n}}:=KrokMH(X_{{n-1}})$ { krok }

end

Oczywista jest analogia z podstawową metodą eliminacji. Zasadnicza różnica polega na tym, że w algorytmie MH nie ,,odrzucamy” zmiennej losowej, tylko ,,odrzucamy propozycje ruchu” – i stoimy w miejscu.

Jądro przejścia M-H jest następujące:

$P(x,B)=\int _{B}{\rm d}y\, q(x,y)a(x,y)+\mathbb{I}(x\in B)\int _{{\cal X}}{\rm d}y\, q(x,y)[1-a(x,y)].$

Dla przestrzeni skończonej, jądro łańcucha MH redukuje się do macierzy prawdopodobieństwprzejścia. Wzór jest w tym przypadku bardzo prosty: dla $x\not=y$ ,

$P(x,y)=q(x,y)a(x,y).$

Twierdzenie 11.2

Jądro przejścia MH jest odwracalne względem $\pi$ .

Dowód

Ograniczmy się do przestrzeni skończonej, żeby nie komplikować oznaczeń. W ogólnym przypadku dowód jest w zasadzie taki sam, tylko napisy stają się mniej czytelne. Niech (bez straty ogólności)

$a(x,y)=\frac{\pi(y)q(y,x)}{\pi(x)q(x,y)}\leqslant 1,\qquad a(y,x)=1.$

Wtedy

$\begin{split}{\pi(x)P(x,y)}&=\pi(x)q(x,y)a(x,y)\\ &=\pi(x)q(x,y)\frac{\pi(y)q(y,x)}{\pi(x)q(x,y)}\\ &=\pi(y)q(y,x)\\ &={\pi(y)P(y,x)}\quad\textrm{ bo }\quad a(y,x)=1.\\ \end{split}$

∎

Uwagi historyczne:

Metropolis w 1953 zaproponował algorytm, w którym zakłada się symetrię rozkładu propozycji, $q(x,y)=q(y,x)$ . Warto zauważyć, że wtedy łańcuch odpowiadający $q$ (błądzenie bez eliminacji ruchów) ma rozkład stacjonarny jednostajny. Reguła akceptacji przybiera postać

$a(x,y)=\frac{\pi(y)}{\pi(x)}\wedge 1.$
Hastings w 1970 uogólnił rozważania na przypadek niesymetrycznego $q$ .

Uwaga ważna:

Algorytm MH wymaga znajomości gęstość $\pi$ tylko z dokładnością do proporcjonalności, bez stałej normującej.

11.3. Próbnik Gibbsa

Drugim podstawowym algorytmem MCMC jest próbnik Gibbsa (PG) (Gibbs Sampler, GS). Załóżmy, że przestrzeń na której żyje docelowy rozkład $\pi$ ma strukturę produktową: ${\cal X}=\prod _{{i=1}}^{{d}}{\cal X}_{{i}}$ . Przyjmijmy następujące oznaczenia:

Jeśli ${\cal X}\ni x=(x_{{i}})_{{i=1}}^{d}$ to $x_{{-i}}=(x_{{j}})_{{j\not=i}}$ : wektor z pominiętą $i$ -tą współrzędną.
Rozkład docelowy (gęstość): $\pi({\rm d}x)=\pi(x){\rm d}x$ .
Pełne rozkłady warunkowe (full conditionals):

$\pi(x_{{i}}|x_{{-i}})=\frac{\pi(x)}{\pi(x_{{-i}})}$

Mały krok PG jest zmianą $i$ -tej współrzędnej (wylosowaniem nowej wartości z rozkładu warunkowego):

$\begin{split} x&=(x_{1},\ldots,x_{i},\ldots,x_{d})\\ &\downarrow\\ \textrm{Gen }y_{i}&\sim\pi(\cdot|x_{{-i}})\\ &\downarrow\\ Y&=(x_{1},\ldots,y_{i},\ldots,x_{d}).\\ \end{split}$

Prawdopodobieństwo przejścia małego kroku PG (w przypadku przestrzeni skończonej) jest takie:

$\qquad P_{{i}}(x,y)=\pi(y_{i}|x_{{-i}})\mathbb{I}(x_{{-i}}=y_{{-i}}).$

Twierdzenie 11.3

Mały krok PG jest $\pi$ -odwracalny.

Dowód

Niech $x_{{-i}}=y_{{-i}}$ . Wtedy

$\qquad\begin{split}{\pi(x)P_{i}(x,y)}&=\pi(x)\pi(y_{i}|x_{{-i}})\\ &=\pi(x_{{-i}})\pi(x_{i}|x_{{-i}})\pi(y_{i}|x_{{-i}})\\ &=\pi(y_{{-i}})\pi(x_{i}|x_{{-i}})\pi(y_{i}|y_{{-i}})\\ &={\pi(y)P_{i}(y,x)}.\end{split}$

(skorzystaliśmy z symetrii).

∎

Oczywiście, trzeba jeszcze zadbać o to, żeby łańcuch generowany przez PG był nieprzywiedlny. Trzeba zmieniać wszystkie współrzędne, nie tylko jedną. Istnieją dwie zasadnicze odmiany próbnika Gibbsa, różniące się spoobem wyboru współrzędnych do zmiany.

Losowy wybór współrzędnych, ,,LosPG”.
Systematyczny wybór współrzędnych, ,,SystemPG”.

Losowy PG. Wybieramy współrzędną $i$ -tą z prawdopodobieństwem $c(i)$ . .

function $LosPG(x)$

Gen $i\sim c(\cdot)$ ;

Gen $y_{{i}}:=\pi(\cdot|x_{{-i}})$ ; { zmieniamy $i$ -tą współrzędną }

$y_{{-i}}:=x_{{-i}}$ ; { wszystkie inne współrzędne pozostawiamy bez zmian }

$LosPG:=y$

Jądro przejścia w ,,dużym” kroku losowego PG jest takie:

$\qquad P=\sum _{{i=1}}^{{d}}c(i)P_{{i}},$

LosPG jest odwracalny.

Systematyczny PG. Współrzędne są zmieniane w porządku cyklicznym.

function $SystemPG(x)$

begin

Gen $y_{{1}}\sim\pi(\cdot|x_{{2}},\ldots,x_{{d}})$ ;

Gen $y_{{2}}\sim\pi(\cdot|y_{{1}},x_{{3}},\ldots,x_{{d}})$ ;

$\cdots$

Gen $y_{{d}}\sim\pi(\cdot|y_{{1}},\ldots,y_{{d-1}})$ ;

$SystemPG:=y$

end

Jądro przejścia w ,,dużym” kroku systematycznego PGjest następujące.

$P=P_{{1}}P_{{2}}\cdots P_{{d}},$

Systematyczny PG nie jest odwracalny. Ale jest $\pi$ -stacjonarny, bo $\pi P_{{1}}P_{{2}}\cdots P_{{d}}=\pi$ .

Przy projektowaniu konkretnych realizacji PG pojawia się szereg problemów, ważnych zarówno z praktycznego jak i teoretycznego punktu widzenia. Jak wybrać rozkład $c(\cdot)$ w losowym PG? Jest raczej jasne, że niektóre współrzędne powinny być zmieniane częściej, a inne rzadziej. Jak dobrać kolejność współrzędnych w systematycznym PG? Czy ta kolejność ma wpływ na tempo zbieżności łańcucha. Wreszczie, w wielu przykładach można zmieniać całe ,,bloki” współrzędnych na raz.

Systematyczny PG jest uważany za bardziej efektywny i częściej stosowany w praktyce. Z drugiej strony jest trudniejszy do analizy teoretycznej, niż losowy PG.

Rys. 11.1. Trajektoria próbnika Gibbsa w przestrzeni dwuwymiarowej.

Chmurka punktów wygenerowanych przez próbnika Gibbsa

Rys. 11.2. Chmurka punktów wygenerowanych przez próbnik Gibbsa.

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

Symulacje stochastyczne i metody Monte Carlo

Zagadnienia

11. Markowowskie Monte Carlo II. Podstawowe algorytmy

11.1. Odwracalność

Definicja 11.1

Twierdzenie 11.1

Dowód

11.2. Algorytm Metropolisa-Hastingsa

Twierdzenie 11.2

Dowód

11.3. Próbnik Gibbsa

Twierdzenie 11.3

Dowód