Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 63 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 65 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 67 Notice: Undefined index: mode in /home/misc/mst/public_html/common.php on line 69 Notice: Undefined variable: base in /home/misc/mst/public_html/lecture.php on line 36 Symulacje stochastyczne – 11. Markowowskie Monte Carlo II. Podstawowe algorytmy – MIM UW

Zagadnienia

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.

Trajektoria próbnika Gibbsa w przestrzeni dwuwymiarowej
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.

Projekt współfinansowany przez Unię Europejską w ramach Europejskiego Funduszu Społecznego.

Projekt współfinansowany przez Ministerstwo Nauki i Szkolnictwa Wyższego i przez Uniwersytet Warszawski.