« 14. Markowowskie Monte Carlo V. Elementarna teoria łańcuchów Markowa  Spis treści  Literatura »

Zagadnienia

15. Markowowskie Monte Carlo VI. Oszacowania dokładności

W tym rozdziale zajmiemy się problemem oszacowania błędu markowowskich algorytmów Monte Carlo. W odróżnieniu od wstępnych rozważań w Rozdziale 10, będą nas interesowały ścisłe nierówności, a nie oceny oparte na twierdzeniach granicznych. Skupimy się na rozkładzie spektralnym macierzy przejścia. Jest to najbardziej znana i zapewne najskuteczniejsza metoda otrzymywania dobrych oszacowań, przynajmniej dla łańcuchów odwracalnych na przestrzeni skończonej.

15.1. Reprezentacja spektralna macierzy odwracalnej

Rozważmy łańcuch nieprzywiedlny i odwracalny z macierzą przejścia P i rozkładem stacjonarnym π⊤. Zakładamy więc, że dla dowolnych x,y∈X spełniona jest zależność

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

Niech

Π=diag⁢π=π⁢10⋯00π⁢2⋯0⋮⋮⋱⋮00⋯π⁢d.

Warunek odwracalności możemy zapisać w postaci Π⁢P=P⊤⁢Π. Wyposażmy przestrzeń Rd w iloczyn skalarny

f,gπ=f⊤⁢Π⁢g=∑xf⁢x⁢g⁢x⁢π⁢x,

gdzie f,g∈R traktujemy jak wektory kolumnowe. Oczywiście, norma funkcji f jest zdefiniowana wzorem fπ2=f,fπ. Macierz P traktujemy jako operator działający z lewej strony na funkcje, z prawej strony na rozkłady prawdopodobieństwa. Zgodnie z regułami mnożenia macierzy i wektorów, wzory na P⁢f i ξ⊤⁢P są następujące:

P⁢f⁢x=∑yP⁢x,y⁢f⁢y,⁢ξ⊤⁢P⁢y=∑xξ⁢x⁢P⁢x,y.

Odwracalność P implikuje

f,P⁢gπ=f⊤⁢Π⁢P⁢g=f⊤⁢P⊤⁢Π⁢g=g⊤⁢Π⁢P⁢f=P⁢f,gπ.

Znaczy to, że P jest macierzą operatora samosprzężonego względem iloczynu skalarnego 〈⋅,⋅〉π. W skrócie powiemy, że P jest π-samosprzężona. Wartości własne P są rzeczywiste i zawarte w przedziale -1,1. Uporządkujmy je w kolejności malejącej:

1=λ0>λ1≥⋯⁢λd-1≥-1.

Wiadomo, że λ0=1 jest pojedynczą wartością własną. Jeśli łańcuch jest nieokresowy, to λd-1>-1. Niech

Λ=diag⁢1,λ1,…,λd-1.

Reprezentacja spektralna macierzy P jest następująca:

P=V⁢Λ⁢V⊤⁢Π,

gdzie

V⊤⁢Π⁢V=I

lub, równoważnie, V⁢V⊤⁢Π=I. Zauważmy, że kolumny macierzy

V=1,v1,…,vd-1

są prawostronnymi wektorami własnymi macierzy P tworzącymi bazę π-ortonormalną. Mamy więc P⁢vi=λ⁢vi (oczywiście, v0=1 jest tu wektorem jedynek) oraz

vi,vjπ=vi⊤⁢Π⁢vj=Ii=j.

Lewostronne wektory własne P (zapisane wierszowo) są postaci vi⊤⁢Π: mamy vi⊤⁢Π⁢P=λi⁢vi⊤⁢Π. W szczególności, v0⊤⁢Π=π⊤. Zapiszmy reprezentację spektralną P w bardziej jawnej formie:

P=∑i≥0λi⁢vi⁢vi⊤⁢Π,

przy tym

∑i≥0vi⁢vi⊤⁢Π=I.

Zauważmy jeszcze, że pierwszy (lub raczej - zerowy) składnik jest macierzą stabilną (o jednakowych wierszach):

v0⁢v0⊤⁢Π=1⁢π⊤.

Reprezentacja spektralna prowadzi do zgrabnych wyrażeń na potęgi macierzy. Łatwo zauważyć, że Pn=V⁢Λn⁢V⊤⁢Π, czyli

Pn=∑i≥0λin⁢vi⁢vi⊤⁢Π=1⁢π⊤+∑i≥1λin⁢vi⁢vi⊤⁢Π. (15.1)

Wiemy, że wszystkie wartości własne λi z wyjątkiem zerowej (λ0=1) oraz, być może, ostatniej (λd-1≥-1) są co do modułu mniejsze niż 1. Jeżeli więc λd-1>-1, to wszystkie składniki sumy we wzorze (15.1) z wyjątkiem początkowego zmierzają do zera i w rezultacie

Pn→1⁢π⊤⁢n→∞.

Jest to nic innego jak teza Słabego Twierdzenia Ergodycznego (Twierdzenie 14.5), otrzymana zupełnie inną metodą, przy założeniu odwracalności. Można pokazać, że warunek λd-1>-1 jest równoważny nieokresowości, i jest konieczny (jeśli λd-1=-1 to łańcuch ma okres 2).

Następujący lemat będzie podstawą dalszych rozważań i umożliwi ,,przerobienie” STE na jawne wyniki. Zdefiniujmy

λ=max⁡λ1,⁢λd-1⁢.

Załóżmy, że łańcuch jest nieokresowy, więc λ<1. Pokażemy, że operator P ograniczony do podprzestrzeni f:f⟂1⊂Rd ortogonalnej do funkcji stałych jest zwężający, ze stałą λ<1.

Lemat 15.1

Jeżeli f,1π=0 to P⁢fπ≤λ⁢fπ.

Dowód

Wystarczy zauważyć, że

⁢P⁢f,P⁢fπ⁢=f,P2⁢fπ⁢⁢=f⊤⁢Π⁢1⁢π⊤⁢f+f⊤⁢∑i≥1λi2⁢n⁢vi⁢vi⊤⁢Π⁢f⁢⁢=∑i≥1λi2⁢n⁢vi,fπ2⁢⁢≤λ2⁢n⁢∑i≥1vi,fπ2⁢⁢=λ2⁢n⁢f,f.⁢

Korzystamy tu z faktu, że P jest samosprzężony, ze wzoru (15.1) i z tego, że f⊤⁢Π⁢1=0.

∎

15.1.1. Oszacowanie szybkości zbieżności

Przejdźmy teraz do jawnych oszacowań szybkości zbieżności w STE. Wyniki zawarte z tym podrozdziale pochodzą z pracy Diaconisa i Strooka [5]. Dla rozkładu prawdopodobieństwa ξ⊤ definiujemy ,,odległość” χ2 od rozkładu stacjonarnego wzorem

χ2⁢π,ξ=∑x⁢ξ⁢x-π⁢x⁢2π⁢x=ξ⊤-π⊤⁢Π-1⁢ξ-π.

Ta ,,odległość” nie ma własności symetrii, więc nie jest metryką, ale to nie przeszkadza. Istotna jest interpretacja χ2 jako ,,odstępstwa od stacjonarności”. Wykorzystamy podejście spektralne, w szczególności Lemat 15.1. Niech χ=χ2.

Stwierdzenie 15.1 (Diaconis i Strook)

Jeśli łańcuch odwracalny, nieprzywiedlny i nieokresowy ma rozkład początkowy ξ, to

χ⁢π,Pn⁢ξ≤λn⁢χ⁢π,ξ.
Dowód

Zastosujmy Lemat 15.1 do wektora Π-1⁢ξ-π który, jak łatwo zauważyć, jest prostopadły do 1. Otrzymujemy

χ⁢π,Pn⁢ξ=Pn⁢Π-1⁢ξ-ππ≤λn⁢Π-1⁢ξ-ππ=λn⁢χ⁢π,ξ.
∎
Wniosek 15.1

Dla łańcucha o rozkładzie początkowym skupionym w punkcie x,

χ(π,Pn(x,⋅))≤λn1-π⁢xπ⁢x≤λnπ⁢x.

Istotnie, mamy jeszcze jedno wyrażenie na ,,odległość” χ2: dla dowolnego rozkładu ξ,

χ2⁢π,ξ=ξ⊤⁢Π-1⁢ξ-1.

Wystarczy teraz podstawić ξ⁢y=Iy=x, aby otrzymać χ2=1/π⁢x-1.

15.1.2. Oszacowanie normy pełnego wahania

Istnieje prosta nierówność pomiędzy normą pełnego wahania i ,,odległością” χ2:

⁢ξ-π⁢tv=12⁢∑x⁢⁢ξ⁢x-π⁢x⁢≤12⁢⁢χ⁢π,ξ.

Wynika to z następującego rachunku:

⁢4⁢⁢⁢ξ-π⁢tv2⁢=⁢∑x⁢⁢ξ⁢x-π⁢x⁢π⁢x⁢⁢π⁢x2⁢⁢≤⁢∑x⁢⁢ξ⁢x-π⁢x⁢2π⁢x⁢⁢∑xπ⁢x⁢⁢[ Cauchy-Schwarz ]⁢⁢=⁢χ2⁢π,ξ.

Stąd natychmiast otrzymujemy wniosek

∥Pn(x,⋅)-π∥tv≤12λn1-π⁢xπ⁢x≤λn2⁢π⁢x.

15.1.3. Oszacowanie obciążenia estymatora

Rozważmy zadanie obliczania wartości oczekiwanej θ=Eπ⁢f=π⊤⁢f dla pewnej funkcji f. Niech f¯=f-1⁢π⊤⁢f oznacza ,,scentrowaną” funkcję f. Natychmiast widać, że f¯⟂1 i możemy zastosować Lemat 15.1. Stąd już tylko mały krok do oszacowania różnicy między wartością oczekiwaną Eξ⁢f⁢Xn=ξ⊤⁢Pn⁢f i wartością stacjonarną θ.

Stwierdzenie 15.2

Jeśli łańcuch jest odwracalny, nieprzywiedlny, nieokresowy i ma rozkład początkowy ξ, to

⁢Eξ⁢f⁢Xn-θ⁢≤λn⁢χ⁢π,ξ⁢σ⁢f,

gdzie σ2⁢f=f¯π2 jest wariancją stacjonarną funkcji f.

Dowód

Mamy

⁢⁢Eξ⁢f⁢Xn-θ⁢⁢=⁢ξ⊤-π⊤⁢Pn⁢f⁢=⁢ξ⊤-π⊤⁢Pn⁢f¯⁢⁢⁢=⁢Π-1⁢ξ-π,Pn⁢f¯π⁢⁢⁢≤⁢Π-1⁢ξ-ππ⁢⁢Pn⁢f¯π⁢⁢[ Cauchy-Schwarz ]⁢⁢≤χ⁢π,ξ⁢λn⁢f¯π⁢⁢⁢[ Lemat 15.1 ].
∎

Rozważmy teraz naturalny estymator

θ⌃t,n=1n⁢∑i=tt+n-1f⁢Xi.

Jest to średnia wzdłuż trajektorii łańcucha, dlugości n i ,,opóźniona” o t. Idea jest jasna: ignorujemy początkowy odcinek trajektorii długości t (tak zwany okres burn-in) aby dać łańcuchowi czas na zbliżenie od rozkładu stacjonarnego. Później obliczmy średnią. W ten sposób redukujemy obciążenie. Precyzuje to następujący wniosek.

Wniosek 15.2

Dla dowolnego n mamy

⁢Eξ⁢θ⌃t,n-θ⁢≤λt1-λ⁢χ⁢π,ξ⁢σ⁢f.

Wynika to z nierówności trójkąta i wzoru na sumę szeregu geometrycznego:

⁢Eξ⁢θ⌃t,n-θ⁢≤∑i=tt+n-1⁢Eξ⁢f⁢Xi-θ⁢≤∑i=t∞λi⁢χ⁢π,ξ⁢σ⁢f.

Zwróćmy uwagę, że obciążenie maleje w tempie geometrycznym przy t→∞ ale zachowuje sią zaledwie jak O⁢1/n przy ustalonym t i n→∞ (ponieważ początkowe wyrazy sumy mają na obciążenie wpływ dominujący).

15.2. Oszacowanie błędu średniokwadratowego estymatora

Oszacowanie błędu średniokwadratowego (BŚK) estymatora MCMC jest znaczne trudniejsze i subtelniejsze, niż obciążenia.

15.2.1. Asymptotyczna wariancja

Zacznijmy od wyprowadzenia kolejnego wzoru na asymptotyczną wariancję. Przypomnijmy, że zgodnie ze Stwierdzeniem 14.1,

σas2⁢f=f⊤⁢C⁢f,

gdzie

⁢C⁢=Π⁢2⁢Z-I-1⁢π⊤=2⁢Π⁢Z-Π-π⁢π⊤⁢⁢=Π⁢2⁢A-I+1⁢π⊤=2⁢Π⁢A-Π+π⁢π⊤.⁢

Skorzystajmy z reprezentacji spektralnej macierzy P:

P-1⁢π⊤=∑i≥1λi⁢vi⁢vi⊤⁢Π=V⁢0⋯0λ1⋯0⋮⋮⋱⋮00⋯λd-1.V⊤⁢Π

Z definicji macierzy A i wzoru (15.1), ponieważ ∑n=0∞λin=1/1-λi, więc

A=∑n=0∞Pn-1⁢π⊤=V⁢0⋯0⋱⋯⋮⋮1/1-λi⋮0⋯⋱⁢V⊤⁢Π.

Wreszcie, ponieważ 2/1-λi-1=1+λi/1-λi, więc

Π⁢2⁢A-1+1⁢π⊤=Π⁢V⁢0⋯0⋱⋯⋮⋮1+λi/1-λi⋮0⋯⋱⁢V⊤⁢Π.

Zauważmy teraz, że wektor α=V⊤⁢Π⁢f zawiera współrzędne wektora f w bazie ON złożonej z prawych wektorów własnych: αi=vi⊤⁢Π⁢f=⁢vi,⁢f⁢π. Udowodniliśmy w ten sposób następujący fakt.

Stwierdzenie 15.3

Dla łańcucha odwracalnego, wzór na asymptotyczną wariancję przybiera postać

σas2⁢f=∑i≥1⁢αi2⁢⁢1+λi1-λi,

gdzie αi=vi⊤⁢Π⁢f=⁢vi,⁢f⁢π.

Wynika stąd ważna nierówność.

Wniosek 15.3

Dla łańcucha odwracalnego mamy następujące oszacowanie asymptotycznej wariancji:

σas2⁢f≤1+λ11-λ1⁢σ2⁢f

gdzie λ1 jest największą wartością własną mniejszą od 1, a σ2⁢f jest wariancją stacjonarną.

Istotnie,

σas2⁢f≤∑i≥1⁢αi2⁢⁢1+λi1-λi≤1+λ11-λ1⁢∑i≥1⁢αi2,

a łatwo widzieć, że ∑i≥1⁢αi2=f¯π2=σ2⁢f. Zwróćmy uwagę, że we Wniosku 15.3 występuje λ1, a nie λ=max⁡λ1,λd-1 (największa co do modułu wartość własna mniejsza od 1).

Na zakończenie przytoczę jeszcze kilka sugestywnych wzorów.

Uwaga 15.1

Macierz L=I-P jest nazywana laplasjanem. Zauważmy, że

L=∑i≥11-λi⁢⁢vi⁢vi⊤⁢Π.

Ponieważ

A=∑i≥111-λi⁢⁢vi⁢vi⊤⁢Π,

uzasadnia to interpretację macierzy A jako ,,uogólnionej odwrotności” laplasjanu.

Dorzućmy przy okazji jeszcze jedno wyrażenie na asymptotyczną wariancję:

⁢σas2⁢f⁢=⁢f,⁢2⁢A-I+1⁢πT⁢f⁢π=⁢f,⁢A⁢f⁢π+⁢A⁢f,⁢f⁢π-σ2⁢f.⁢

Jeśli π⊤⁢f=0, to σ2⁢f=Varπ⁢f=⁢f,f⁢π i ostatni wzór możemy przepisać w postaci

σas2⁢f=σ2⁢f⁢2⁢⁢f,⁢A⁢f⁢π⁢f,⁢f⁢π-1

15.2.2. Oszacowanie BŚK

Wyniki w tym podrozdziale zostały otrzymane przez Aldousa [1]. Pomysł polega na tym, żeby najpierw otrzymać nierówność dla łańcucha stacjonarnego, a potem postarać się o uogólnienie dla łańcucha o dowolnym rozkładzie początkowym. Przypomnijmy oznaczenia θ=Eπ⁢f, θ⌃n=∑i=0n-1f⁢Xi.

Stwierdzenie 15.4 (Aldous, 1987)

Dla dla łańcucha nieprzywiedlnego, odwracalnego i stacjonarnego, MSEπ można oszacować w następujący sposób:

Eπ⁢θ⌃n-θ2≤1+λ~1-λ~⁢σ2⁢f

gdzie λ~=max⁡λ1,0, zaś λ1 jest największą wartością własną mniejszą od 1.

Dowód

Korzystamy ze stacjonarności i z rozkładu spektralnego macierzy P.

Eπ(θ⌃n-θ)2=Eπ(1n∑i=0n-1f¯(Xi))2=1n2∑i=0n-1Eπf¯(Xi)2+2n2∑i=0n-1∑j=i+1n-1Eπf¯(Xi)f¯(Xj)=1nEπf¯(X0)2+2n2∑k=1n-1(n-k)Eπf¯(X0)f¯(Xk)[k=j-i]=1n〈f¯,f¯〉π+2n2∑k=1n-1(n-k)〈f¯,Pkf¯〉π=1n〈f¯,f¯〉π+2n2∑k=1n-1(n-k)∑i≥1λ~ik〈vi,f〉π2=1n〈f¯,f¯〉π+2n2∑i≥1〈vi,f〉π2∑k=1n-1(n-k)λik≤1nσ2(f)+2n2∑i≥1λi>0〈vi,f〉π2∑k=1n-1(n-k)λik[pomijamy składniki dla λi<0]≤1nσ2(f)+2n2σ2(f)∑k=1n-1(n-k)λ~k≤1nσ2(f)+2nσ2(f)∑k=1n-1λ~k≤2nσ2(f)(1+2⁢λ~1-λ~)=1nσ2(f)1+λ~1-λ~.

Zwróćmy uwagę na miejsce, w którym pomijamy składniki odpowiadające ujemnym wartościom własnym. Uzasadnienie jest takie, że dla , λi<0, składniki sumy ∑k=1n-1n-k⁢λik są naprzemian ujemne i dodatnie, o malejących wartościach bezwzględnych. Stąd wynika, że ∑k=1n-1n-k⁢λik<0 i można tę sumę w nierówności opuścić. Jest to ciekawe zjawisko: ujemne wartości własne pomagają, zmniejszają błąd! Działają podobnie jak zmienne antytetyczne.

∎

Niech teraz E⁢x oznacza błąd średniokwadratowy dla łańcucha startującego z punktu x∈X,

en⁢x=Ex⁢θ⌃n-θ2.

Rozpatrzmy łańcuch o rozkładzie początkowym ξ. Niech

θ⌃t,n=1n⁢∑i=tt+n-1f⁢Xi.

będzie średnią długości n obliczany po odrzuceniu t początkowych zmiennych.

Stwierdzenie 15.5

Dla dla łańcucha nieprzywiedlnego, odwracalnego startującego z dowolnego rozkładu ξ, MSEξ można oszacować w następujący sposób:

Eξ⁢θ⌃t,n-θ2≤1+λ~1-λ~⁢σ2⁢f+λt⁢χ⁢π,ξ⁢maxx⁡f¯⁢x2

gdzie λ~=max⁡λ1,0 i λ=max⁡λ1,λd-1.

Dowód

Na mocy Stwierdzeń 15.4 i 15.2 mamy

⁢Eξ⁢θ⌃t,n-θ2⁢=Eξ⁢en⁢Xt≤Eπ⁢en+Eξ⁢en⁢Xt-Eπ⁢en⁢⁢≤1n⁢σ2⁢f⁢1+λ~1-λ~+λt⁢χ⁢π,ξ⁢σ⁢en⁢⁢≤1n⁢σ2⁢f⁢1+λ~1-λ~+λt⁢χ⁢π,ξ⁢maxx⁡f¯⁢x2,

bo en⁢y≤maxx⁡f¯⁢x2, a więc σ⁢en≤maxx⁡f¯⁢x2.

∎

Nierówność podane przez Aldousa w cytowanej pracy była nieco inna.

Stwierdzenie 15.6 (Aldous, 1987)

Dla dla łańcucha nieprzywiedlnego, odwracalnego startującego z dowolnego rozkładu ξ mamy następujące oszacowanie błędu M⁢S⁢Eξ:

Eξ⁢θ⌃t,n-θ2≤1+λtπ*⁢1+λ~1-λ~⁢σ2⁢f,

gdzie π*=minx⁡π⁢x.

Dowód

Istotnie, załóżmy, że łańcuch startuje z deterministycznie wybranego punktu x∈X, czyli ξ=P(x,⋅). Jeśli otrzymamy oszacowanie niezależne od x, dowód będzie zakończony.

⁢Ex⁢θ⌃t,n-θ2⁢=∑yPt⁢x,y⁢en⁢y⁢⁢=∑yPt⁢x,y⁢en⁢y=∑yPt⁢x,yπ⁢y⁢π⁢y⁢en⁢y⁢⁢≤1+λtπ*⁢∑yπ⁢y⁢en⁢y=1+λtπ*⁢MSEπ.

i zastosujmy Stwierdzenie 15.4. Wykorzystaliśmy tu nierówność Pt⁢x,y/π⁢y≤1+λt/π*, która wynika z rozkładu spektralnego:

⁢Pt⁢x,yπ⁢y⁢=Pt⁢Π-1⁢x,y=1x⊤⁢Pt⁢Π-1⁢1y⁢⁢=1x⊤⁢∑i=0d-1λit⁢vi⁢vi⊤⁢1y=1+1x⊤⁢∑i=1d-1λit⁢vi⁢vi⊤⁢1y⁢⁢≤1+λt⁢∑i=1d-11x⊤⁢vi⁢vi⊤⁢1y⁢⁢≤1+λt⁢∑i=1d-11x⊤⁢vi2⁢∑i=1d-1vi⁢1y⊤2⁢⁢≤1+λtπ⁢x⁢π⁢y≤1+λtπ*.

W powyższym wzorze symbol 1x⊤ oznacza wektor 0,…,0,1,0,…,0, gdzie jedynka stoi na x-tym miejscu. Skorzystaliśmy z nierówności Schwarza.

∎
« 14. Markowowskie Monte Carlo V. Elementarna teoria łańcuchów Markowa  Spis treści  Literatura »

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

strona główna | webmaster | o portalu | pomoc

© Wydział Matematyki, Informatyki i Mechaniki UW, 2009-2010. Niniejsze materiały są udostępnione bezpłatnie na licencji Creative Commons Uznanie autorstwa-Użycie niekomercyjne-Bez utworów zależnych 3.0 Polska.

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.