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 π, jeśli dla dowolnych A,B⊆X mamy

∫Aπ⁢d⁢x⁢P⁢x,B=∫Bπ⁢d⁢y⁢P⁢y,A.

W skrócie,

⁢π⁢d⁢x⁢P⁢x,d⁢y=π⁢d⁢y⁢P⁢y,d⁢x.

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

Twierdzenie 11.1

Jeśli π⁢d⁢x⁢P⁢x,d⁢y=π⁢d⁢y⁢P⁢y,d⁢x to π⁢P=P.

Dowód

π⁢P⁢B=∫Xπ⁢d⁢x⁢P⁢x,B=∫Bπ⁢d⁢y⁢P⁢y,X=∫Bπ⁢d⁢y=π⁢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 π. 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: π⁢d⁢x=π⁢x⁢d⁢x.

  • Rozkład ,,propozycji”: q⁢x,d⁢y=q⁢x,y⁢d⁢y.

  • Reguła akceptacji:

    a⁢x,y=π⁢y⁢q⁢y,xπ⁢x⁢q⁢x,y∧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 π. Pojedynczy krok algorytmu jest następujący.

function K⁢r⁢o⁢k⁢M⁢H⁢x
          Gen y∼q(x,⋅);     { propozycja }
          Gen U∼U⁢0,1
      if U>a⁢x,y}  then
        begin
        y:=x;     { ruch odrzucony z pr-stwem 1-a⁢x,y }
        end
   K⁢r⁢o⁢k⁢M⁢H:=y
Graficznie to można przedstawić w takiej postaci.

Xn=x[d][dl]a⁢x,yy∼q(x,⋅)[dr]1-a⁢x,yXn+1=yXn+1=x

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

Gen X0∼π0;  { start }
    for n:=1 to ∞
        begin
        Xn:=K⁢r⁢o⁢k⁢M⁢H⁢Xn-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)=∫Bdyq(x,y)a(x,y)+I(x∈B)∫Xdyq(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≠y,

P⁢x,y=q⁢x,y⁢a⁢x,y.
Twierdzenie 11.2

Jądro przejścia MH jest odwracalne względem π.

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=π⁢y⁢q⁢y,xπ⁢x⁢q⁢x,y⩽1,⁢a⁢y,x=1.

Wtedy

⁢π⁢x⁢P⁢x,y⁢=π⁢x⁢q⁢x,y⁢a⁢x,y⁢⁢=π⁢x⁢q⁢x,y⁢π⁢y⁢q⁢y,xπ⁢x⁢q⁢x,y⁢⁢=π⁢y⁢q⁢y,x⁢⁢=π⁢y⁢P⁢y,x⁢⁢bo ⁢⁢a⁢y,x=1.⁢⁢
∎

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=π⁢yπ⁢x∧1.
  • Hastings w 1970 uogólnił rozważania na przypadek niesymetrycznego q.

Uwaga ważna:

  • Algorytm MH wymaga znajomości gęstość π 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 π ma strukturę produktową: X=∏i=1dXi. Przyjmijmy następujące oznaczenia:

  • Jeśli X∋x=xii=1d to x-i=xjj≠i: wektor z pominiętą i-tą współrzędną.

  • Rozkład docelowy (gęstość): π⁢d⁢x=π⁢x⁢d⁢x.

  • Pełne rozkłady warunkowe (full conditionals):

    π(xi|x-i)=π⁢xπ⁢x-i

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

x=(x1,…,xi,…,xd)↓Gen yi∼π(⋅|x-i)↓Y=(x1,…,yi,…,xd).

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

Pi(x,y)=π(yi|x-i)I(x-i=y-i).
Twierdzenie 11.3

Mały krok PG jest π-odwracalny.

Dowód

Niech x-i=y-i. Wtedy

π(x)Pi(x,y)=π(x)π(yi|x-i)=π(x-i)π(xi|x-i)π(yi|x-i)=π(y-i)π(xi|x-i)π(yi|y-i)=π(y)Pi(y,x).

(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 L⁢o⁢s⁢P⁢G⁢x
      Gen i∼c(⋅);
      Gen yi:=π(⋅|x-i); { zmieniamy i-tą współrzędną }
      y-i:=x-i; { wszystkie inne współrzędne pozostawiamy bez zmian }
    L⁢o⁢s⁢P⁢G:=y

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

⁢P=∑i=1dc⁢i⁢Pi,

LosPG jest odwracalny.

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

function S⁢y⁢s⁢t⁢e⁢m⁢P⁢G⁢x
     begin
     Gen y1∼π(⋅|x2,…,xd);
     Gen y2∼π(⋅|y1,x3,…,xd);
     ⋯
     Gen yd∼π(⋅|y1,…,yd-1);
     S⁢y⁢s⁢t⁢e⁢m⁢P⁢G:=y
     end

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

P=P1⁢P2⁢⋯⁢Pd,

Systematyczny PG nie jest odwracalny. Ale jest π-stacjonarny, bo π⁢P1⁢P2⁢⋯⁢Pd=π.

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(⋅) 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.