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 – 13. Markowowskie Monte Carlo IV. Pola losowe – MIM UW

Zagadnienia

13. Markowowskie Monte Carlo IV. Pola losowe

13.1. Definicje

Niech ({\cal S},{\cal E}) będzie nieskierowanym grafem. Wyobraźmy sobie, że elementy s\in S reprezentują ,,miejsca” w przestrzeni lub na płaszczyźnie, zaś krawędzie grafu łączą miejsca ,,sąsiadujące” ze sobą. Taka interpretacja jest związana z zastosowaniami do statystyki ,,przestrzennej” i przetwarzania obrazów. Model, który przedstawimy ma również zupełnie inne interpretacje, ale pozostaniemy przy sugestywnej terminologii ,,przestrzennej”:

  • {\cal S} — zbiór miejsc,

  • \{ s,t\}\in{\cal E} — miejsca s i t sąsiadują — będziemy wtedy pisać s\sim t,

  • \partial t=\{ s:s\sim t\} — zbiór sąsiadów miejsca t.

Niech \Lambda=\{ 1,\ldots,l\} będzie skończonym zbiorem. Powiedzmy, że elementy a\in\Lambda są ,,kolorami” które mogą być przypisane elementom zbioru {\cal S}. Konfiguracją nazywamy dowolną funkcję x:{\cal S}\to\Lambda. Będziemy mówić że x_{s}=x(s) jest s-tą współrzędną konfiguracji x i stosować oznaczenia podobne jak dla wektorów:

x=(x_{{s}})=(x_{{s}}:s\in{\cal S}). (13.1)

W zadaniach przetwarzania obrazów, miejsca są pikslami na ekranie i konfigurację utożsamiamy z ich pokolorowaniem, a więc z cyfrową reprezentacją obrazu. Zbiór \Lambda gra rolę ,,palety kolorów”. Niekiedy założenie o skończoności zbioru \Lambda staje się niewygodne. Dla czarno-szaro-białych obrazów ,,pomalowanych” różnymi odcieniami szarości, wygodnie przyjąć, że \Lambda=[0,1] lub \lambda=[0,\infty[. Tego typu modyfikacje są dość oczywiste i nie będę się nad tym zatrzymywał. Dla ustalenia uwagi, wzory w tym podrozdziale dotyczą przypadku skończonego zbioru ,,kolorów”. Przestrzenią konfiguracji jest zbiór {\cal X}=\Lambda^{{{\cal S}}}. Dla konfiguracji x i miejsca t, niech

  • x_{{-t}}=(x_{{s}}:s\not=t)=(x_{{s}}:s\in{\cal S}\setminus\{ t\}) — konfiguracja z pominiętą t-tą współrzędną,

  • x_{{\partial t}}=(x_{{s}}:s\in\partial t) — konfiguracja ograniczona do sąsiadów miejsca t.

Jeśli H:{\cal X}\to\mathbb{R} i \beta\geqslant 0 to rozkładem Gibbsa nazywamy rozkład prawdopodobieństwa na przestrzeni konfiguracji dany wzorem

\pi _{{\beta}}(x)=\frac{1}{Z(\beta)}\exp[-\beta H(x)].

Ze względu na inspiracje pochodzące z fizyki statystycznej, funkcję H nazywamy energią, \beta jest (z dokładnością do stałej) odwrotnością temperatury. Stała normująca wyraża się wzorem

Z(\beta)=\sum _{{x\in{\cal X}}}\exp[-\beta H(x)]. (13.2)

i jest typowo niemożliwa do obliczenia.

Oczywiście, każdy rozkład prawdopodobieństwa \pi na {\cal X} dale się zapisać jako rozkład Gibbsa, jeśli położyć H(x)=-\log\pi(x), \beta=1 i umownie przyjąć, że -\log 0=\infty (czyli konfiguracje niemożliwe mają nieskończoną energię). Nie o to jednak chodzi. Ciekawe są rozkłady Gibbsa, dla których fukcja energii ma specjalną postać związaną z topologią grafu ,,sąsiedztw”. Ograniczymy się do ważnej podklasy markowowskich pól losowych (MPL), mianowicie do sytuacji gdy energia jest sumą ,,oddziaływań” lub ,,interakcji” między parami miejsc sąsiadujących i składników zależnych od pojedynczych miejsc. Dokładniej, założymy że

H(x)=\sum _{{s\sim t}}V(x_{{s}},x_{{t}})+\sum _{{s}}U_{{s}}(x_{{s}}), (13.3)

dla pewnych funkcji V:\Lambda\times\Lambda\to\mathbb{R} i U_{{s}}:\Lambda\to\mathbb{R}. Funkcja V(x_{{s}},x_{{t}}) opisuje ,,potencjał interakcji pomiędzy s i t”, zaś U_{{s}}(a) jest wielkością związaną z ,,tendencją miejsca s do przybrania koloru a. Zwróćmy uwagę, że potencjał V jest jednorodny (V(a,b) zależy tylko od ,,kolorów” a,b\in\Lambda ale nie od miejsc), zaś U_{s}(a) może zależeć zarówno od a\in\Lambda jak i od s\in{\cal S}. W modelach fizyki statystycznej zazwyczaj U_{s}(a)=U(a) jest jednorodnym ,,oddziaływaniem zewnętrznym” ale w modelach rekonstrukcji obrazów nie można tego zakładać.

Przykład 13.1 (Model Pottsa)

Niech \Lambda będzie zbiorem skończonym i

H(x)=J\sum _{{s\sim t}}\mathbb{I}(x_{{s}}\not=x_{{t}}).

Ta funkcja opisuje ,,tendencję sąsiednich miejsc do przybierania tego samego koloru”. Jeśli J>0 to preferowane są konfiguracje złożone z dużych, jednobarwnych plam.

13.2. Generowanie markowowskich pól losowych

Użyteczność MPL w różnorodnych zastosowaniach związana jest z istnieniem efektywnych algorytmów symulacyjnych. Wszystko opiera się na próbniku Gibbsa i następującym prostym fakcie.

Twierdzenie 13.1 (Pełne rozkłady warunkowe dla MPL)

Jeżeli \pi _{{\beta}} jest rozkładem Gibbsa z energią daną wzorem (13.3), to

\pi _{{\beta}}(x_{{s}}|x_{{-s}})=\pi _{{\beta}}(x_{{s}}|x_{{\partial s}})=\frac{1}{Z_{{s}}(\beta)}\exp\left[-\beta H_{{s}}(x)\right], (13.4)

gdzie

H_{{s}}(x)=\sum _{{t\in\partial s}}V(x_{s},x_{t})+U(x_{s}), (13.5)

Z_{{s}}(\beta)=\sum _{{a\in\Lambda}}\exp[H_{s}(x_{{a\rightsquigarrow s}}). Symbol x_{{a\rightsquigarrow s}} oznacza konfigurację powstałą z x przez wpisanie koloru a w miejscu s.

Dowód

Skorzystamy z elementarnej definicji prawdopodobieństwa warunkowego (poniżej piszemy \pi _{{\beta}}(\cdot)=\pi(\cdot), bo parametr \beta jest ustalony):

\begin{split}\pi(x_{{s}}|x_{{-s}})&=\frac{\pi(x)}{\pi(x_{{-s}})}=\frac{\pi(x)}{\sum\limits _{a}\pi(x_{{a\rightsquigarrow s}})}\\
&=\frac{\exp-\beta H(x)}{\sum _{a}\exp-\beta H(x_{{a\rightsquigarrow s}})}\\
&=\frac{\exp-\beta\left(\sum\limits _{{t:t\sim s}}V(x_{{s}},x_{{t}})+\sum\limits _{{t\sim w,t\not=s,w\not=s}}V(x_{{t}},x_{{w}})+U_{s}(x_{s})+\sum\limits _{{t\not=s}}U_{{t}}(x_{{t}})\right)}{\sum\limits _{a}\exp-\beta\left(\sum\limits _{{t:t\sim s}}V(a,x_{{t}})+\sum\limits _{{t\sim w,t\not=s,w\not=s}}V(x_{{t}},x_{{w}})+U_{s}(a)+\sum\limits _{{t\not=s}}U_{{t}}(x_{{t}})\right)}\\
&=\frac{\exp-\beta\left(\sum\limits _{{t:t\sim s}}V(x_{{s}},x_{{t}})+U_{s}(x_{s})\right)}{\sum\limits _{a}\exp-\beta\left(\sum\limits _{{t:t\sim s}}V(a,x_{{t}})+U_{s}(a)\right)}\\
&=\frac{\exp-\beta H_{s}(x)}{Z_{s}(\beta)}.\end{split} (13.6)

Ponieważ otrzymany wynik zależy tylko od x_{s} i x_{{\partial s}}, więc \pi(x_{{s}}|x_{{-s}})=\pi(x_{{s}}|x_{{\partial s}}). Ten wniosek jest pewną formą własności Markowa.

Zauważmy, że obliczenie H_{s}(x) jest łatwe, bo suma \sum _{{t\in\partial s}}\cdots zawiera tylko tyle składników, ile jest sąsiadów miejsca s. Obliczenie Z_{s}(\beta) też jest łatwe, bo suma \sum _{{a\in\Lambda}}\cdots zawiera tylko l=|\Lambda| składników. Ale nawet nie musimy obliczać stałej normującej Z_{s}(\beta) żeby generować z rozkładu

\pi(x_{s}=a|x_{{\partial s}})\propto exp-\beta\left(\sum\limits _{{t:t\sim s}}V(a,x_{{t}})+U_{s}(a)\right). (13.7)

Na tym opiera się implementacja próbnika Gibbsa. Wersję PG z ,,systemstycznym przeglądem miejsc” można zapisać tak:

for  s\in{\cal S} do
    begin
    Gen a\sim\pi(x_{s}=\cdot|x_{{\partial s}});
    x:=x_{{a\rightsquigarrow s}}
    end
Faktycznie już ten algorytm spotkaliśmy w Podrozdziale 12.2, dla szczególnego przypadku modelu auto-logistycznego.

13.3. Rekonstrukcja obrazów

Bayesowski model rekonstrukcji obrazów został zaproponowany w pracy Gemana i Gemana w 1987 roku. Potem zdobył dużą popularność i odniósł wiele sukcesów. Model łączy idee zaczerpnięte ze statystyki bayesowskiej i fizyki statystycznej. Cyfrową reprezentację obrazu utożsamiamy z konfiguracją kolorów na wierzchołkach grafu, czyli z elementem przestrzeni {\cal X}=\Lambda^{{\cal S}}, zdefiniowanej w Podrozdziale 13.1. Przyjmijmy, że ,,idealny obraz”, czyli to co chcielibyśmy zrekonstruować jest konfiguracją x\in{\cal X}. Niestety, obraz jest ,,zakłócony” lub ,,zaszumiony”. Możemy tylko obserwować konfigurację y reprezentującą zakłócony obraz. Zbiór kolorów w obrazie y nie musi być identyczny jak w obrazie x. Powiedzmy, że y\in{\cal Y}=\Gamma^{{\cal S}}. Ważne jest to, że zniekształcenie modelujemy probabilistycznie przy pomocy rodziny rozkładów warunkowych f(y|x). Dodatkowo zakładamy, że obraz x pojawia się losowo, zgodnie z rozkładem prawdopodobieństwa \pi(x). Innymi słowy, ,,idealny” obraz x oraz ,,zniekształcony” obraz y traktujemy jako realizacje zmiennych losowych X:\Omega\to{\cal X} i Y:\Omega\to{\cal Y},

\pi(x)=\mathbb{P}(X=x),\qquad f(y|x)=\mathbb{P}(Y=y|X=x). (13.8)

W ten sposób buduje się statystyczny model bayesowski, w którym

  • Y jest obserwowaną zmienną losową,

  • x jest nieznanym parametrem traktowanym jako zmienna losowa X.

Oczywiście, \pi gra rolę rozkładu a priori, zaś f jest wiarogodnością. Być może użycie literki x na oznaczenie parametru jest niezgodne z tradycyjnymi oznaczeniami statystycznymi, ale z drugiej strony jest wygodne. Wzór Bayesa mówi, że rozkład a posteriori jest następujący.

\pi _{y}(x)=\mathbb{P}(X=x|Y=y)\propto f(y|x)\pi(x). (13.9)

Pomysł Gemana i Gemana polegał na tym, żeby modelować rozkład a priori \pi jako MPL. Załóżmy, że \pi jest rozkładem Gibbsa,

\pi(x)\propto\exp\left(-H(x)\right), (13.10)

gdzie

H(x)=J\sum _{{s\sim t}}V(x_{{s}},x_{{t}}). (13.11)

Energia ,,a priori” zawiera tu tylko składniki reprezentujące oddziaływania między parami miejsc sąsiednich. Funkcja V(a,b) zazwyczaj ma najmniejszą wartość dla a=b i rośnie wraz z ,,odległością” między a i b (jakkolwiek tę odległość zdefiniujemy). W ten sposób ,,nagradza” konfiguracje w których sąsiedznie miejsca są podobnie pokolorowane. Im większy parametr J>0, tym bardziej prawdopodobne są obrazy zawierające jednolite plamy kolorów.

Trzeba jeszcze założyć coś o ,,wiarogodności” f. Dla uproszczenia opiszę tylko najprostszy model, w którym kolor y_{s} na obserwowanym obrazie zależy tylko od koloru x_{s} na obrazie idealnym. Intuicyjnie znaczy to, że ,,zaszumienie” ma ściśle lokalny charakter. Matematycznie znaczy to, że

f(y|x)=\prod _{{s}}f(y_{{s}}|x_{{s}}) (13.12)

(pozwolę sobie na odrobinę nieścisłości aby uniknąć nowego symbolu na oznaczenie f(y_{{s}}|x_{{s}})). Zapiszemy teraz ,,wiarogodność” f w postaci zlogarytmowanej. Jeśli położymy -\log f(y_{s}|x_{s})=U_{s}(x_{s}) pamiętając, że y jest w świecie bayesowskim ustalone, to otrzymujemy następujący wzór:

\pi _{y}(x)\propto\exp\left(-H_{y}(x)\right), (13.13)

gdzie

H_{y}(x)=J\sum _{{s\sim t}}V(x_{{s}},x_{{t}})-\sum _{{s}}\log f(y_{s}|x_{s}). (13.14)

Okazuje się zatem, że rozkład a posteriori ma podobną postać do rozkładu a priori. Też jest rozkładem Gibbsa, a różnica polega tylko na dodaniu składników reprezentujących oddziaływania zewnętrzne U_{s}(x_{s})=-\log f(y_{s}|x_{s}). Pamiętajmy przy tym, że y jest w świecie bayesowskim ustalone. W modelu rekonstrukcji obrazów ,,oddziaływania zewnętrzne” zależą od y i ,,wymuszają podobieństwo” rekonstruowanego obrazu do obserwacji. Z kolei ,,oddziaływania między parami” są odpowiedzialne za wygładzenie obrazu. Lepiej to wyjaśnimy na przykładzie.

Przykład 13.2 (Losowe ,,przekłamanie koloru” i wygładzanie Pottsa)

Załóżmy, że \Lambda=\{ 1,\ldots,l\} jest naprawdę paletą kolorów, na przykład \Lambda=\{{Czerwony},{Niebieski},{Pomarańczowy},{Zielony}\}.

Przypuśćmy, że mechanizm losowego ,,przekłamania” polega na tym, że w każdym pikslu, kolor obecny w idealnym obrazie x jest z prawdopodobieństwem 1-\varepsilon niezmieniony, a z prawdopodobieństwem \varepsilon zmienia się na losowo wybrany inny kolor. Tak więc zarówno x jak i y należą do tej samej przestrzeni \Lambda^{{\cal S}},

f(y_{s}|x_{s})=\begin{cases}1-\varepsilon&\text{ dla }y_{s}=x_{s};\\
\varepsilon/(l-1)&\text{ dla }y_{s}\not=x_{s}.\end{cases} (13.15)

Można za rozkład a priori przyjąć rozkład Pottsa z Przykładu 13.1. Rozkład a posteriori ma funkcję energii daną następującym wzorem (z J>0):

H_{y}(x)=J\sum _{{s\sim t}}\mathbb{I}(x_{{s}}\not=x_{{t}})-\sum _{{s}}\left[\log(1-\varepsilon)\mathbb{I}(x_{s}=y_{s})+\log(\varepsilon/(l-1))\mathbb{I}(x_{s}\not=y_{s})\right]. (13.16)

Pierwszy składnik w tym wzorze pochodzi od rozkładu a priori (z modelu Pottsa) i ,,nagradza” konfiguracje w których dużo sąsiednich punktów jest pomalowanych na ten sam kolor. Powoduje to, że obrazy x składające się z jednolotych dużych ,,plam” są preferowane. Drugi składnik pochodzi od obserwowanegj konfiguracji y i jest najmniejszy dla x=y. Powoduje to, ze obrazy x mało się różniące od y są bardziej prawdopodobne. Rozkład a posteriori jest pewnym kompromisem pomiędzy tymi dwoma konkurującymi składnikami. Parametr J jest ,,wagą” pierwszego składnika i dlatego odgrywa rolę ,,parametru wygładzającego”. Im większe J tym odtwarzany obraz będzie bardziej regularny (a tym mnie będzie starał się upodobnić do y). I odwrotnie, małe J powoduje ściślejsze dopasowanie x do y ale mniejszą ,,regularność” x.

Jeszcze lepiej to samo widać na przykładzie tak zwanego ,,szumu gaussowskiego”.

Przykład 13.3 (Addytywny szum gaussowski)

Załóżmy, że x jest konfiguracją ,,poziomów szarości” czyli, powiedzmy, \Lambda\subseteq[0,\infty[. Mechanizm losowego ,,zaszumienia” polega na tym, że zamiast poziomu szarości x_{s} obserwujemy y_{s}\sim{\rm N}(x_{s},\sigma^{2}). Innymi słowy,

f(y_{s}|x_{s})\propto\exp\left[-\frac{1}{2\sigma^{2}}(y_{s}-x_{s})^{2}\right]. (13.17)

Przestrzenią obserwaowanych konfiguracji y jest tutaj (formalnie) \mathbb{R}^{{\cal S}} (faktycznie, raczej [0,\infty[^{{\cal S}}). Rozkład a posteriori ma funkcję energii daną następującym wzorem:

H_{y}(x)=J\sum _{{s\sim t}}V(x_{{s}}\not=x_{{t}})+\frac{1}{2\sigma^{2}}\sum _{{s}}(y_{s}-x_{s})^{2}. (13.18)

Jeśli rozpatrujemy model ze skończoną liczbą poziomów szarości dla konfiguracji x to można pierwszy składnik określić tak jak w poprzednim przykładzie, czyli zapożyczyć z modelu Pottsa. Bardziej naturalne jest określenie V(a,b) w taki sposób, aby większe różnice pomiędzy poziomami a i b były silniej karane. parametr J jest, jak poprzednio, odpowiedzialny zaa stopień wygładzenia.

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.