Zagadnienia

13.1 Definicje
13.2 Generowanie markowowskich pól losowych
13.3 Rekonstrukcja obrazów

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.

Symulacje stochastyczne i metody Monte Carlo

Zagadnienia

13. Markowowskie Monte Carlo IV. Pola losowe

13.1. Definicje

Przykład 13.1 (Model Pottsa)

13.2. Generowanie markowowskich pól losowych

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

Dowód

13.3. Rekonstrukcja obrazów

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

Przykład 13.3 (Addytywny szum gaussowski)