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,yX spełniona jest zależność

πxPx,y=πyPy,x.

Niech

Π=diagπ=π1000π2000πd.

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

f,gπ=fΠg=xfxgxπx,

gdzie f,gR 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 Pf i ξP są następujące:

Pfx=yPx,yfy,ξPy=xξxPx,y.

Odwracalność P implikuje

f,Pgπ=fΠPg=fPΠg=gΠPf=Pf,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

Λ=diag1,λ1,,λd-1.

Reprezentacja spektralna macierzy P jest następująca:

P=VΛVΠ,

gdzie

VΠV=I

lub, równoważnie, VVΠ=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 Pvi=λ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=λiviΠ. W szczególności, v0Π=π. Zapiszmy reprezentację spektralną P w bardziej jawnej formie:

P=i0λiviviΠ,

przy tym

i0viviΠ=I.

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

v0v0Π=1π.

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

Pn=i0λinviviΠ=1π+i1λinviviΠ. (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

Pn1π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:f1Rd ortogonalnej do funkcji stałych jest zwężający, ze stałą λ<1.

Lemat 15.1

Jeżeli f,1π=0 to Pfπλfπ.

Dowód

Wystarczy zauważyć, że

Pf,Pfπ=f,P2fπ=fΠ1πf+fi1λi2nviviΠf=i1λi2nvi,fπ2λ2ni1vi,fπ2=λ2nf,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-πx2π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=12xξx-πx12χπ,ξ.

Wynika to z następującego rachunku:

4ξ-πtv2=xξx-πxπxπx2xξx-πx2πxxπx[ Cauchy-Schwarz ]=χ2π,ξ.

Stąd natychmiast otrzymujemy wniosek

Pn(x,)-πtv12λ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ξfXn=ξPnf i wartością stacjonarną θ.

Stwierdzenie 15.2

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

EξfXn-θλnχπ,ξσf,

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

Dowód

Mamy

EξfXn-θ=ξ-πPnf=ξ-πPnf¯=Π-1ξ-π,Pnf¯πΠ-1ξ-ππPnf¯π[ Cauchy-Schwarz ]χπ,ξλnf¯π[ Lemat 15.1 ].

Rozważmy teraz naturalny estymator

θt,n=1ni=tt+n-1fXi.

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-1EξfXi-θi=tλiχπ,ξσf.

Zwróćmy uwagę, że obciążenie maleje w tempie geometrycznym przy t ale zachowuje sią zaledwie jak O1/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,

σas2f=fCf,

gdzie

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

Skorzystajmy z reprezentacji spektralnej macierzy P:

P-1π=i1λiviviΠ=V00λ1000λd-1.VΠ

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

A=n=0Pn-1π=V001/1-λi0VΠ.

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

Π2A-1+1π=ΠV001+λi/1-λi0VΠ.

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ć

σas2f=i1αi21+λ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:

σas2f1+λ11-λ1σ2f

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

Istotnie,

σas2fi1αi21+λi1-λi1+λ11-λ1i1αi2,

a łatwo widzieć, że i1αi2=f¯π2=σ2f. 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=i11-λiviviΠ.

Ponieważ

A=i111-λiviviΠ,

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

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

σas2f=f,2A-I+1πTfπ=f,Afπ+Af,fπ-σ2f.

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

σas2f=σ2f2f,Afπ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-1fXi.

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-θ21+λ~1-λ~σ2f

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π(1ni=0n-1f¯(Xi))2=1n2i=0n-1Eπf¯(Xi)2+2n2i=0n-1j=i+1n-1Eπf¯(Xi)f¯(Xj)=1nEπf¯(X0)2+2n2k=1n-1(n-k)Eπf¯(X0)f¯(Xk)[k=j-i]=1nf¯,f¯π+2n2k=1n-1(n-k)f¯,Pkf¯π=1nf¯,f¯π+2n2k=1n-1(n-k)i1λ~ikvi,fπ2=1nf¯,f¯π+2n2i1vi,fπ2k=1n-1(n-k)λik1nσ2(f)+2n2i1λi>0vi,fπ2k=1n-1(n-k)λik[pomijamy składniki dla λi<0]1nσ2(f)+2n2σ2(f)k=1n-1(n-k)λ~k1nσ2(f)+2nσ2(f)k=1n-1λ~k2nσ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 Ex oznacza błąd średniokwadratowy dla łańcucha startującego z punktu xX,

enx=Exθn-θ2.

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

θt,n=1ni=tt+n-1fXi.

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-θ21+λ~1-λ~σ2f+λtχπ,ξmaxxf¯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ξenXtEπen+EξenXt-Eπen1nσ2f1+λ~1-λ~+λtχπ,ξσen1nσ2f1+λ~1-λ~+λtχπ,ξmaxxf¯x2,

bo enymaxxf¯x2, a więc σenmaxxf¯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 MSEξ:

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

gdzie π*=minxπx.

Dowód

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

Exθt,n-θ2=yPtx,yeny=yPtx,yeny=yPtx,yπyπyeny1+λtπ*yπyeny=1+λtπ*MSEπ.

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

Ptx,yπy=PtΠ-1x,y=1xPtΠ-11y=1xi=0d-1λitvivi1y=1+1xi=1d-1λitvivi1y1+λti=1d-11xvivi1y1+λti=1d-11xvi2i=1d-1vi1y21+λtπxπy1+λ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.

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.