Markovketten Monte-Carlo Methoden - Friedrich-Schiller

Markovketten Monte-Carlo Methoden
Dr. Daniel Rudolf
24. September 2015
Dieses Vorlesungsskript basiert auf der Lehrveranstaltung Markovketten Monte-Carlo Methoden,
gehalten im Sommersemester 2015 an der Friedrich-Schiller-Universität Jena. Ich bedanke mich
herzlich bei Stephan Wolf, der die Bilder erstellte und den Hauptteil des Textsatzes übernahm.
Inhaltsverzeichnis
1 Motivation
2
2 Bedingter Erwartungswert und bedingte Verteilung
2.1 Bedingter Erwartungswert . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Bedingte Verteilung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Bedingte Verteilungsdichten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
9
13
15
3 Markovketten
3.1 Definitionen und erste Eigenschaften . . . . . . . . . . .
3.2 Grundlegende Algorithmen . . . . . . . . . . . . . . . .
3.2.1 Metropolis-Hastings Algorithmus . . . . . .
3.2.2 Hit-and-run-Algorithmus . . . . . . . . . . . . .
3.2.3 Gibbs sampler . . . . . . . . . . . . . . . . . . .
3.2.4 Slice sampling . . . . . . . . . . . . . . . . . . .
3.3 Markovoperator . . . . . . . . . . . . . . . . . . . . . . .
3.4 Klassische Konvergenzeigenschaften von Markovketten
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
18
18
31
33
37
40
42
45
49
4 Markovketten Monte-Carlo Methoden
4.1 Fehlerabschätzungen . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Hinreichende Bedingungen für explizite Konvergenzabschätzungen
4.2.1 Doeblin-Bedingung . . . . . . . . . . . . . . . . . . . . . . .
4.2.2 Leitfähigkeit und Cheeger-Ungleichung . . . . . . . . . . .
4.3 Anwendungen und Beispiele . . . . . . . . . . . . . . . . . . . . . .
4.3.1 Integration über konvexe Mengen . . . . . . . . . . . . . . .
4.3.2 Integration bezüglich nichtnormierter Dichtefunktion . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
55
56
68
68
70
80
80
86
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Literatur
103
Index
104
1
1
Motivation
• (Ω, F, P) Wraum
∗
• Sei G ⊂ Rd beschränkt, d ∈ N, λd d-dimensionales Lebesguemaß und 0 < λ(G) < ∞
(∗, da G beschränkt)
• B(G) sei die σ-Algebra der Borelmengen (bzgl. Euklidischer Metrik)
• (G, B(G)) messbarer Raum
R
• f : G → R messbar mit G |f (x)|dx < ∞
Ziel: Berechne
S(f ) =
1
λd (G)
Z
f (x)dx
G
(S = Solution-Operator) für f ∈ F(G), wobei
F(G) = {g : G → R : ∫G |g(x)|dx < ∞ + weitere Eigenschaften}.
Bemerkung 1.1
S(f ) ist der Erwartungswert (EW) von f bezüglich der Gleichverteilung in G.
Gegebene Informationen / Annahmen:
(i) Wir besitzen eine „black box“ (auch Orakel genannt), die uns Funktionswerte von f liefert, d.h.
für jedes x ∈ G können wir nach f (x) fragen.
(ii) Wir nehmen an, dass wir einen Zufallszahlengenerator bezüglich der Gleichverteilung in G besitzen, d.h. wir haben eine Folge von unabhängig identisch verteilten (iid, bzw. uiv) Zufallsvariablen
(ZVen) (Xi )i∈N , wobei
Xi : (Ω, F, P) → (G, B(G))
mit
PXi (A) =
λd (A)
,
λd (G)
A ∈ B(G).
Und wir können für n ∈ N nach einer Realisierung xi = Xi (ω), i = 1, . . . , n fragen (ω ∈ Ω fest)
und erhalten eine Antwort.
Bemerkung 1.2
Die zweite Annahme wird meist durch die Transformation von Zufallsvariablen mit „einfachen“ Verteilungen sichergestellt.
Beispiel
(i) G = [0, 1]d . Sei x1 , . . . , xd ∈ [0, 1] eine Folge von Realisierungen von iid ZVen, die in [0, 1]
gleichverteilt sind. Dann ist x = (x1 , . . . , xd ) ∈ [0, 1]d eine Realisierung der Gleichverteilung in
G.
Pd
(ii) G = Bd := {x ∈ Rd : |x| := ( i=1 x2i )1/2 ≤ 1}. Wir benötigen folgende Notation:
• Sd−1 = {x ∈ Rd : |x| = 1} = ∂Bd die Einheitssphäre (Rand der Einheitskugel)
• σ Oberflächenmaß auf Sd−1
2
• d-dimensionale Polarkoordinaten-Transformationsformel: f : Rd → R sei integrierbar.
Dann gilt
Z
Z∞ Z
f (x)dx =
f (θr)rd−1 σ(dθ)dr.
x=θr
0 Sd−1
Rd
(siehe [5, Satz 8 in Kapitel 14])
Lemma 1.3
Sei X : (Ω, F, P) → (Rd , B(Rd )) eine ZV mit rotationsinvarianter Lebesguedichte h : Rd → R, d.h.
∀x, y ∈ Rd : |x| = |y| =⇒ h(x) = h(y).
Desweiteren sei U : (Ω, F, P) → ([0, 1], B([0, 1])) eine gleichverteilte ZV in [0, 1], unabhängig von X.
Dann gilt P({|X| = 0}) = 0 und
X
(i) P( |X|
∈ A) =
σ(A)
σ(Sd−1 ) ,
X
∈ A) =
(ii) P(U 1/d |X|
für A ∈ B(Sd−1 ) (Gleichverteilung auf der Einheitssphäre).
λd (A)
λd (Bd ) ,
für A ∈ B(Bd ) (Gleichverteilung in der Einheitskugel).
Notation: {X ∈ A} = {ω ∈ Ω : X(ω) ∈ A}
Beweis: Es gilt P({|X| = 0}) = P(X ∈ {0}) = 0, da {0} eine Nullmenge bezüglich des Lebesguemaßes
ist. Im Folgenden beweisen wir die Aussagen (i) und (ii).
(i)
P(X/|X| ∈ A) = E[1{X/|X|∈A} ]
Z
= 1A (x/|x|)h(x)dx
Rd
Z∞
Z
1A (θ)h(θr)rd−1 σ(dθ)dr
=
0 Sd−1
Z∞
Z
1A (θ)σ(dθ) ·
=
h(θr)rd−1 dr
(θ ∈ Sd−1 bel. wg. Rotat.inv.)
0
Sd−1
Z∞
= σ(A)
h(θr)rd−1 dr.
0
Sei für den Moment A = Sd−1 . Dann gilt
Z∞
1 = σ(Sd−1 )
h(θr)rd−1 dr.
0
Dies führt für beliebiges A ∈ B(Sd−1 ) zu
P(X/|X| ∈ A) =
σ(A)
.
σ(Sd−1 )
3
(ii) Es gilt
P(U 1/d X/|X| ∈ A) = E[1{U 1/d X/|X|∈A} ]
Z1 Z
1A (u1/d θ)
=
σ(dθ)
du
σ(Sd−1 )
0 Sd−1
{z
|
Z1
1/d
g(u
=
Z1
)du =
0
=
}
=g(u1/d )
d
σ(Sd−1 )
dtd−1 g(t)dt
0
Z1 Z
1A (tθ)σ(dθ)td−1 dt
0 Sd−1
=
d
σ(Sd−1 )
Z∞ Z
1A (tθ)σ(dθ)td−1 dt
0 Sd−1
=
d
σ(Sd−1 )
Z
1A (x)dx = λd (A) ·
d
σ(Sd−1 )
Rd
λd (A)
,
=
λd (Bd )
da λd (Bd ) =
σ(Sd−1 )
.
d
(Zur Erinnerung: λd (Bd ) =
π d/2
Γ(d/2+1) ,
d/2
π
σ(Sd−1 ) = d Γ(d/2+1)
, siehe [5, Kapitel 7, unter (7.4)].)
Algorithmus zur Realisierung gleichverteilter ZV in Bd :
(i) Sei x = (x1 , . . . , xd ) ∈ Rd , wobei die Folge der xi eine Realisierung von iid Standardnormal2
verteilten ZVen ist (ZV hat Dichte g(s) = √12π e−s /2 ). Da die Normalverteilung im Rd ein
Produktmaß ist, ist x eine Realisierung der d-dimensionalen Standardnormalverteilung. Diese
2
ist rotationsinvariant (ZV hat die Dichte g(x) = (2π)1d/2 e−|x| /2 ).
(ii) y = x/|x|, d.h. y ist Realisierung der Gleichverteilung in Sd−1 .
(iii) Sei u ∈ [0, 1] eine Realisierung der Gleichverteilung in [0, 1], dann ist z = u1/d y eine Realisierung
der Gleichverteilung in Bd .
Bemerkung 1.4
Sei a ∈ Rd , r > 0 und B(a, r) = {x ∈ Rd : |x − a| ≤ r}. B(a, r) ist die Kugel mit Radius r um a. Dann
x
ist a + ru1/d |x|
(mit u aus (iii) und x aus (i)) eine Realisierung einer gleichverteilten ZV in B(a, r).
Zusammenfassend können wir sagen, dass die zweite Annahme für „gute“ Mengen G ⊂ Rd erfüllt ist.
Klassische Monte-Carlo Methode Sei (Ω, F, P) Wahrscheinlichkeitsraum und X1 , . . . , Xn iidFolge von ZVen mit
Xi : (Ω, F, P) → (G, B(G))
4
und
PXi (A) = P({Xi ∈ A}) =
Dann ist
λd (A)
,
λd (G)
A ∈ B(G).
n
Sn (f ) =
1X
f (Xi )
n i=1
die klassische Monte-Carlo Methode zur Approximation von
Z
1
S(f ) =
f (x)dx.
λd (G)
G
Satz 1.5 (Starkes Gesetz der
R großen Zahlen)
Sei f : G → R mit kf k1 = G |f (x)| λddx
(G) < ∞. Dann gilt
P( lim Sn (f ) = S(f )) = 1.
n→∞
Beweis: Siehe [4, Thm. 8.3.5].
Dies zeigt, dass die klassische Monte-Carlo Methode wohldefiniert ist. Die Approximation stellt einen
„konsistenten“ Schätzer dar.
Fragen:
(i) Fehlerbegriff?
(ii) Fehlerabschätzung?
Sei An (f ) irgendeine Methode/Formel zur Approximation von S(f ), die n ZVen (nicht notwendigerweise iid) und n Funktionswerte von f benutzt. Dann betrachten wir folgenden Fehler:
Definition 1.6
Der Fehler von An (f ) für f : G → R mit kf k1 < ∞ ist durch
e(An , f ) := (E|An (f ) − S(f )|2 )1/2
gegeben (E bezüglich der Verteilung von ZV An (f )).
Bemerkung 1.7
(i) e(An , f )2 ist der mittlere quadratische Fehler (mean square error, MSE).
(ii) Es gilt
e(An , f )2 = Var(An (f )) + Bias(An (f )),
wobei
Var(An (f )) = E|An (f ) − EAn (f )|2 ,
Bias(An (f )) = |EAn (f ) − S(f )|2 .
Die klassische Monte-Carlo Methode Sn (f ) hat wegen Bias(Sn (f )) = 0 die Eigenschaft der Erwartungstreue, d.h. E(Sn (f )) = S(f ) und es gilt e(Sn , f ) = Var(Sn (f ))1/2 .
5
Satz 1.8 (Klassische Monte-Carlo Methode, Fehlerabschätzung)
R
1/2
Sei f : G → R mit kf k2 := G |f (x)|2 λddx
< ∞. Dann gilt
(G)
e(Sn , f ) = E|Sn (f ) − S(f )|2
1/2
kf k2
= (Var Sn (f ))1/2 ≤ √ .
n
Beweis (siehe [17, Abschnitt 3.2.2]):
Da X1 , . . . , Xn iid ist, ist auch f (X1 ), . . . , f (Xn ) iid. Dann gilt


n
X
1
e(Sn , f )2 = Var 
f (Xj )
n j=1
=
n
1 X
Var(f (Xi )) Formel von Bienayme, siehe [4, Thm. 8.1.2]
n2 i=1
=
n
1 X
kf − S(f )k22
n2 i=1
=
kf k22
1
kf − S(f )k22 ≤
.
n
n
Wobei die letzte Ungleichung wegen Folgendem gilt:
Z
dx
kf − S(f )k22 = |f (x) − S(f )|2
λd (G)
G
Z
=
f (x)2 − 2f (x)S(f ) + S(f )2
dx
λd (G)
G
= kf k22 − S(f )2 ≤ kf k22 .
Jetzt: Modifizierte Informationen
1) Funktionsauswertungen von f
2’) G ⊂ Rd ist durch ein Orakel, ein „Membership“-Orakel, gegeben, d.h. für x ∈ Rd können wir nach
(
1 x∈G
1G (x) =
0 x∈
/G
fragen und erhalten eine Antwort.
Regularitätsannahme an G: Sei r ≥ 1 und Bd ⊆ G ⊆ rBd , wobei rBd = {x ∈ Rd : |x| ≤ r} und
Bd = 1Bd .
Bemerkung 1.9
D.h. G ist nicht beliebig groß und nicht beliebig klein. Um die klassische Monte-Carlo Methode
anzuwenden müssen wir einen Zufallszahlengenerator konstruieren, der unter 2’) die Gleichverteilung
erzeugt.
6
Erster Ansatz: Verwerfungsmethode
i) Erzeuge z ∈ rBd gleichverteilt (siehe Bemerkung 1.4).
ii) Falls z ∈ G, dann akzeptiere z, d.h. gebe z aus, ansonsten gehe zu i).
• Was ist die Akzeptanzwahrscheinlichkeit in ii)?
• Was ist die erwartete Anzahl von Orakelaufrufen von 1G ?
Im schlechtesten Fall ist G = Bd . Dann ist die Akzeptanzwahrscheinlichkeit:
λd (Bd )
= r−d ,
λd (rBd )
da λd (rBd ) = rd λd (Bd ). Die erwartete Anzahl von Aufrufen von 1G ist rd .
Für große Dimensionen d ist r−d exponentiell klein und rd exponentiell groß. =⇒ Ansatz schlägt fehl.
Vielleicht ist Schritt i) zu schlecht?
√
Sei nun r = d und G = Bd Zweiter Ansatz: Ersetze i) durch i’):
i’) Erzeuge x ∈ [−1, 1]d gleichverteilt.
D.h. i’) ausführen und dann zu ii), falls x verworfen wurde, gehe wieder zu i’).
Skizze:
[−1, 1]d
Bd
rBd
[−1, 1]d ist die kleinste Box, die Bd enthält. Die Akzeptanzwahrscheinlichkeit ist
√ (d+1)/2
λd (Bd )
2e
eπ
, was sich asymptotisch wie
verhält.
2d
π
d+1
Dies ist exponentiell klein in d. Die erwartete Anzahl von Aufrufen von 1G ist
nentiell groß in d. =⇒ Der Ansatz schlägt fehl.
2d
λd (Bd ) .
Das ist expo-
Neue Idee: Markovketten Monte-Carlo Methoden
Mittels einer Folge von nicht notwendigerweise unabhängigen ZVen, einer Markovkette,
X1 , X2 , . . . , Xn+n0
wollen wir die Gleichverteilung in G approximieren, d.h.
P(Xn0 ∈ A) −→
n0 →∞
Dann berechne für f : G → R
λd (A)
,
λd (G)
A ∈ B(G).
(1)
n
Sn,n0 (f ) =
1X
f (Xj+n0 )
n j=1
7
als Approximation von S(f ) =
1
λd (G)
R
G
f (x)dx.
Wir stellen uns nun u.a. folgende Fragen:
1) Konvergenz in (1) quantifizierbar, d.h. Abschätzung zw. dem „Abstand“ der Verteilungen?
2) Ist dieser Algorithmus / diese Methode Sn,n0 konsistent? (Starkes Gesetz der großen Zahlen)
3) Fehlerabschätzung?
4) Information ausreichend? Bzw. welche Information wird benötigt?
8
2
Bedingter Erwartungswert und bedingte Verteilung
Dieser Abschnitt basiert im Wesentlichen auf [4, Kap. 10], siehe auch [1].
2.1
Bedingter Erwartungswert
Seien (Ω, F, P) Wraum, X : (Ω, F, P) → (R, B(R)), X ∈ L1 , d.h.
R
Ω
|X|dP < ∞.
Sei B ∈ F fest und P(B) > 0.
Zur Erinnerung: Für A ∈ F ist
P(A | B) =
P(A ∩ B)
.
P(B)
Es gilt
Z
P(A | B) =
1A∩B (ω)
Ω
d.h.
dP(ω)
,
P(B)
1B (ω)
dP(· | B)
(ω) =
.
dP
P(B)
Dann ist der elementare bedingte EW (bedingt auf B) gegeben durch
Z
Z
X(ω)
E[1B X]
E[X | B] = X(ω)dP(ω | B) =
· 1B (ω)dP(ω) =
.
P(B)
P(B)
Ω
Ω
Verallgemeinerung: Erwarungswert unter der Bedingung E, wobei E ⊆ F σ-Algebra.
Intuitiv: Der Erwartungswert von X unter E sollte die „beste Vorhersage“ von X sein, wenn uns
die Information E zur Verfügung steht (EW als „Vorhersage“, zusätzliche Information sollte diesen
präzisieren).
Definition 2.1 (Bedingter Erwartungswert unter E)
Eine ZV Z : (Ω, F, P) → (R, B(R)) heißt Version des bedingten EW von X unter E, falls
I) Z ist E-messbar (d.h. ∀A ∈ B(R) : Z −1 (A) ∈ E),
R
R
II) ∀A ∈ E : A ZdP = A XdP.
Notation: E[X | E] = Z
Satz 2.2 (Existenz und Eindeutigkeit
P-f.s., [4, Thm. 10.1.1])
R
Für beliebige X ∈ L1 (d.h. Ω |X|dP < ∞) und E ⊆ F σ-Algebra existiert eine Version des
bedingten EW unter E und alle Versionen Z1 , Z2 erfüllen Z1 = Z2 P-f.s.
Beispiel 2.3
B ∈ F, 0 < P(B) < 1, E := {∅, B, B c , Ω} die von B erzeugte σ-Algebra,
Z(ω) := 1B (ω)E[X | B] + 1B c (ω)E[X | B c ].
Wir zeigen, dass E[X | E] = Z P-f.s.
9
2.1
Bedingter Erwartungswert
I) Z ist E-messbar: Sei A ∈ B(R),

Ω



∅
Z −1 (A) =

B


 c
B
II) ∀A ∈ E :
R
A
ZdP =
R
A
E[X
E[X
E[X
E[X
| B], E[X | B c ] ∈ A
| B], E[X | B c ] ∈
/A
.
| B] ∈ A, E[X | B c ] ∈
/A
| B] ∈
/ A, E[X | B c ] ∈ A
XdP: A = ∅X, BX, B c X, ΩX.
Spezialfall: Sei Y : (Ω, F, P) → (G, G) ZV, (G, G) messbarer Raum und sei
Y −1 (G) = {Y −1 (B)|B ∈ G}.
Desweiteren sei σ(Y ) die kleinste σ-Algebra bezüglich der Y messbar ist, d.h.
σ(Y ) := σ(Y −1 (G)).
Es gilt
1) Y −1 (G) ⊆ F σ-Algebra,
2) σ(Y ) = Y −1 (G).
Verallgemeinerung: Seien Yi : (Ω, F, P) → (G, G) für i = 1, . . . , r ZVen. Es sei σ(Y1 , . . . , Yr ) die
kleinste σ-Algebra bezüglich der Y1 , . . . , Yr messbar, d.h.
σ(Y1 , . . . , Yr ) := σ({Yi−1 (A)|A ∈ G, i = 1, . . . , r}).
Definition 2.4
E[X | Y ] := E[X | σ(Y )] und E[X | Y1 , . . . , Yr ] = E[X | σ(Y1 , . . . , Yr )].
Satz 2.5 (Faktorisierungslemma, siehe [4, Thm. 4.2.8])
Sei Z : (Ω, F, P) → (R, B(R)) eine σ(Y )-messbare ZV. Dann ex. eine messbare Abbildung F :
(G, G) → (R, B(R)), s.d. Z = F (Y ).
(Ω, F, P)
Y
(G, G)
∃F : Z = F (Y )
Z
(R, B(R))
Folgerung 2.6
Es existiert eine Abbildung F : (G, G) → (R, B(R)), s.d. E[X | Y ] = F (Y ) P-f.s. und wenn F1 , F2
solche Funktionen sind, dann folgt F1 (Y ) = F2 (Y ) P-f.s..
Definition 2.7
Die Zahl F (y) heißt bedingter Erwartungswert von X unter der Bedingung Y = y, y ∈ G.
Notation: E[X | Y = y] := F (y), E[X | Y = ·] = F (·).
Vergleich der Definitionen
10
2.1
Z : (Ω, F, P) → (R, B(R))
Z = E[X | E], E ⊆ F σ-Algebra
Bedingter Erwartungswert
F : (G, G) → (R, B(R))
F (y) = E[X | Y = y]
Zusammenhang: F (Y ) = Z und E = σ(Y ).
Satz 2.8
Gegeben sei eine messbare Abbildung F : (G, G) → (R, B(R)). Dann gilt
F (y) = E[X | Y = y] PY -f.s. (y ∈ G)
genau dann, wenn
Z
∀B ∈ G :
Z
F (y)PY (dy) =
B
XdP.
Y
−1 (B)
Beweis: Zunächst sei bemerkt, dass F (Y ) : (Ω, F, P) → (R, B(R)) σ(Y )-messbar ist. Dies gilt, weil
∀A ∈ B(R) : (F (Y ))−1 (A) = Y −1 (F −1 (A)) ∈ σ(Y ).
| {z }
∈G
Weiter gilt
F (y) = E[X | Y = y] PY -f.s. ⇐⇒ F (Y ) = E[X | Y ] P-f.s.
(
I) F (Y )
σ(Y ) − messbar
R
R
⇐⇒
II) ∀A ∈ σ(Y ) A F (Y )dP = A XdP
A=Y −1 (B), B∈G
Z
Z
⇐⇒
Z
XdP =
Y −1 (B)
F (Y )dP =
Y −1 (B)
F (y)dPY (y) (Transf.-Satz)
B
Eigenschaften des bedingten EW:
E1 E = {∅, Ω} =⇒ E[X | E] = E[X]
Beweis: klar
E2 Sei X E-messbar =⇒ E[X | E] = X P-f.s.
Beweis: Z = X erfüllt I) und II).
E3 E1 ⊆ E2 ⊆ F σ-Algebren =⇒ E[X | E2 ] ∈ L1 und E[E[X | E2 ] | E1 ] = E[X | E1 ] P-f.s.
Beweis: Sei Z = E[X | E1 ]. Wir prüfen I) und II).
Zu I: Z ist E1 messbar: klar, da Z = E[X | E1 ].
Zu II: ∀A ∈ E1 gilt
Z
Z
A
II
Z
E[X | E1 ] =
ZdP =
A
XdP
A
=⇒ Z ist bedingter EW von E[X | E2 ] unter E1 .
E1 ⊆E2 , A∈E2 ,II
Z
E[X | E2 ]dP.
=
A
11
2.1
Bedingter Erwartungswert
E4 X ∈ L1 , Y E-messbar, XY ∈ L1 =⇒ E[XY | E] = Y · E[X | E] P-f.s.
Beweis (Idee):
Zu I: Klar, da Y · E[X | E] als Komposition von E-messbaren Funktionen wieder E-messbar.
Zu II: (Standardbeweistechnik „Algebraische Induktion“)
• Zeige Behauptung für Y = 1A , A ∈ F.
• Zeige Behauptung für Y einf. Funktion.
• Zeige Behauptung für Y pos. Funktion.
• Zeige Behauptung für Y allg. Funktion.
E5 E[X | E] ∈ L1 und E[E[X | E]] = E[X]
Beweis: Sei E1 = {∅, Ω}, E2 = E. Dann gilt
E1
E3
E1
E[E[X | E]] = E[E[X | E2 ] | E1 ] = E[X | E1 ] = E[X].
E6 Sei Y : (Ω, F, P) → (G, G) ZV und U : (Ω, F, P) → (E, E) ZV. Weiter
sei Y (stochastisch)
R
unabhängig von U und F : (G × E, G ⊗ E) → (R, B(R)) messbar, mit Ω |F (Y, U )|dP < ∞. Dann
gilt
Z
E[F (Y, U ) | Y ](ω) =
F (Y (ω), U (ω))dP(ω) P-f.s.
Ω
Anders ausgedrückt: E[F (Y, U ) | Y = y] = E[F (y, U )], y ∈ G.
Beweis (Idee):
R
Setze Z(ω) = Ω F (Y (ω), U (ω))dP(ω) und prüfe I) und II).
Zu I: Z ist σ(Y )-messbar, wegen F (Y (ω), U (ω)) ist σ(Y )⊗σ(U )-messbar und dem Satz von Fubini
(siehe [9, Satz 14.16]).
Zu II: Sei A ∈ σ(Y ), d.h. A = Y −1 (B), B ∈ G. Wir erhalten
Z
Z
Z
ZdP = 1A (ω)Z(ω) dP(ω) = 1B (Y (ω))Z(ω) dP(ω)
A
Ω
Ω
Z Z
1B (Y (ω)) · F (Y (ω), U (ω)) dP(ω)dP(ω)
=
Ω Ω
= E[1B (Y ) · F (Y, U )] (Unabhängigkeit von Y, U )
Z
= 1B (Y (ω))F (Y (ω), U (ω)) dP(ω)
Ω
Z
1A (ω)F (Y (ω), U (ω)) dP(ω)
=
Ω
Z
F (Y, U ) dP.
=
A
12
2.2
Bedingte Verteilung
2.2
Bedingte Verteilung
Sei (Ω, F, P) Wraum, E ⊆ F σ-Algebra. Es gilt: 1A ∈ L1 ∀A ∈ F, d.h. E[1A | E] ist definiert.
Definition 2.9
P(A | E)(ω) := E[1A | E](ω) heißt Version der bedingten Wahrscheinlichkeit von A unter der
Bedingung E.
Satz 2.10 (siehe [4, Abschnitt 10.2])
1) P(Ω | E) = 1 P-f.s.
2) ∀A ∈ F: P(A | E) ≥ 0 P-f.s.
3) ∀An ∈ F paarweise disjunkt gilt P (
S∞
n=1
An | E) =
P∞
n=1
P(An | E) P-f.s.
Bemerkung 2.11
Die Ausnahmemengen für die Gültigkeit von 2) und 3) sind abhängig von A bzw. An . Insbesondere
ist P(A | E)(ω) im Allgemeinen kein Wmaß für alle ω ∈ Ω, auch nicht P-f.s..
Definition 2.12 (reguläre bedingte Wahrscheinlichkeit)
Eine reguläre Version der bedingten Wahrscheinlichkeit von P unter E ist eine Funktion P E :
Ω × F → [0, 1], s.d.
(i) ∀A ∈ F : P E (·, A) E-messbar,
(ii) ∀ω ∈ Ω : P E (ω, ·) Wmaß auf (Ω, F),
(iii) ∀A ∈ F : P E (·, A) = P(A | E) P-f.s.
(Aus (i) und (iii) folgt, dass P E eine Version der bedingten Wkt. von A unter E ist.)
e A
e ∈ G.
Sei Y : (Ω, F, P) → (G, G) Zufallsvariable. Betrachte Ereignisse A = {Y ∈ A},
Definition 2.13 (reguläre bedingte Verteilung von Y unter E)
Eine reguläre Version der bedingten Verteilung von Y unter der Bedingung E ist eine Funktion
P Y |E : Ω × G → [0, 1], s.d.
(i) ∀A ∈ G : P Y |E (· | A) E-messbar,
(ii) ∀ω ∈ Ω : P Y |E (ω, ·) Wmaß auf (G, G),
(iii) ∀A ∈ G : P Y |E (·, A) = P({Y ∈ A} | E) P-f.s.
(Aus (i) und (iii) folgt, dass P Y |E eine Version der bedingten Wkt. von {Y ∈ A} unter E ist.)
Im Folgenden betrachten wir nur noch reguläre Versionen von bedingten Verteilungen. Es folgt: Zusammenhang zwischen regulären, bedingten Verteilungen und dem bedingten Erwartungswert.
Satz 2.14 (siehe [4, Thm 10.2.5])
Seien (Ω, F, P) Wraum, E ⊆ F σ-Algebra, Y : (Ω, F, P) → (G, G) ZV und P Y |E ist eine reguläre
Version der bedingten Verteilung von Y unter E. Dann gilt für alle messbaren Funktionen g :
13
2.2
Bedingte Verteilung
(G, G) → (R, B(R)) mit
R
G
|g(y)|PY (dy) < ∞, dass g integrierbar bezüglich P Y |E (ω, ·) ist, d.h.
Z
|g(y)|P Y |E (ω, dy) < ∞
G
für P-fast alle ω ∈ Ω, und
Z
E[g(Y ) | E](·) =
g(y)P Y |E (·, dy),
P-f.s..
G
(D.h.
R
G
g(y)P Y |E (·, dy) ist eine Version des bedingten EW von g(Y ) unter E.)
Frage: Existieren reguläre Versionen der bedingten Verteilung? Eindeutigkeit?
Satz 2.15 (siehe [4, Thm 10.2.2])
Seien (Ω, F, P) Wraum, E ⊆ F σ-Algebra, G metrischer, vollständiger, separablera Raum, B(G)
die σ-Algebra der Borelmengen, und Y : (Ω, F, P) → (G, B(G)) ZV. Dann existiert eine reguläre
Version der bedingten Verteilung von Y unter E und diese ist P-f.s. eindeutig.
a Ein topologischer Raum heißt separabel, wenn es eine höchstens abzählbare Teilmenge gibt, die in diesem Raum
dicht liegt.
Für uns ist der folgende Spezialfall am wichtigsten:
X : (Ω, F, P) → (H, H),
wobei (H, H) messbarer Raum und E = σ(X).
Definition 2.16 (reguläre bedingte Verteilung von Y unter X)
Seien X : (Ω, F, P) → (H, H), Y : (Ω, F, P) → (G, G) ZVen. Dann heißt
P Y |X : H × G → [0, 1]
eine reguläre Version der bedingten Verteilung von Y unter X, falls
(i) ∀A ∈ G : P Y |X (·, A) H-messbar,
(ii) ∀x ∈ H : P Y |X (x, ·) Wmaß auf (G, G),
(iii) ∀A ∈ G : P Y |X (·, A) = P({Y ∈ A} | X = ·) PX -f.s..
Satz 2.17 (teilweise Folgerung von Satz 2.14)
Sei P Y |X eine reguläre Version der bedingten
Verteilung von Y unter X. Dann gilt für alle messbaren
R
Funktionen g : (G, G) → (R, B(R)) mit G |g(x)|PY (dy) < ∞, dass g integrierbar bzgl. P Y |X (x, ·)
ist und
Z
E[g(Y ) | X = x] = g(y)P Y |X (x, dy)
(2)
G
für PX -fast alle x ∈ H.
Desweiteren gilt für alle messbaren Funktionen
ge : (H × G, H ⊗ G) → (R, B(R))
14
2.3
Bedingte Verteilungsdichten
mit E|g(X, Y )| < ∞, dass
(i) ge(x, ·) integrierbar bzgl. P Y |X (x, ·) für PX -fast alle x ∈ H,
R R
(ii) E[e
g (X, Y )]] = H G ge(x, y)P Y |X (x, dy)PX (dx),
R
(iii) E[e
g (X, Y ) | X = x] = G ge(x, y)P Y |X (x, dy) für PX -fast alle x ∈ H.
Beweis (Idee, bzw. Referenzen):
• (2) ist Folgerung aus Satz 2.14.
• (iii) folgt aus [8, Thm 5.4], siehe auch [1].
• (ii) folgt aus (iii) und
E[e
g (X, Y )] = E[E[e
g (X, Y ) | X]]


Z
X|Y
= E  g(X(·), y)P
(X(·), dy)
Z ZG
=
g(x, y)P X|Y (x, dy)PX (dx) (Transformationssatz)
H G
• (i) folgt aus (ii)
Bei der Anwendung von Satz 2.17 ist meist (H, H) = (G, G), d.h. X, Y : (Ω, F, P) → (G, G) ZVen.
2.3
Bedingte Verteilungsdichten
Sei (Ω, F, P) Wraum, r ∈ N und Z : (Ω, F, P) → (Rr , B(Rr )) ein zufälliger Vektor, d.h. Z =
(Z1 , . . . , Zr ) mit Zufallsvariablen Zi : (Ω, F, P) → (R, B(R)) für i = 1, . . . , r.
Definition 2.18 (Verteilungsdichte)
Eine messbare Funktion fZ : Rr → [0, ∞) heißt Verteilungsdichte (VD) von Z bzgl. (dem Lebesguemaß) λr , falls
Z
P(Z ∈ A) =
fZ (z)λr (dz),
A ∈ B(Rr ).
A
Seien s, d ∈ N und
Y : (Ω, F, P) → (Rs , B(Rs ))
X : (Ω, F, P) → (Rd , B(Rd ))
mit Y = (Y1 , . . . , Ys ), X = (X1 , . . . , Xd ) und Xi , Yj : (Ω, F, P) → (R, B(R)) ZVen für i = 1, . . . , d und
j = 1, . . . , s, d.h. X und Y zufällige Vektoren.
Satz 2.19
Sei f(Y,X) : Rs × Rd → [0, ∞) eine VD von (Y, X) bezüglich λs+d . Dann gilt
15
2.3
Bedingte Verteilungsdichten
(i)
Z
fY (y) =
f(Y,X) (y, x)λd (dx),
y ∈ Rs
f(Y,X) (y, x)λs (dy),
x ∈ Rd
Rd
ist VD von Y bzgl. λs und
Z
fX (x) =
Rs
ist VD von X bzgl. λd .
(ii) Eine reguläre Version der bedingten Verteilung P X|Y von X unter Y ist gegeben durch
Z
P X|Y (y, A) = fX|Y =y (x)λd (dx)
A
für y ∈ Rs , A ∈ B(Rd ) mit
(f
(Y,X) (y,x)
fX|Y =y (x) =
fY (y)
fX (x)
fY (y) > 0
fY (y) = 0
,
x ∈ Rd .
Definition 2.20
fX|Y =y heißt bedingte Dichte von X unter Y = y.
Beweis: Wir beweisen die Aussagen von (i) und (ii) separat:
(i) fY : Rs → [0, ∞), fX : Rd → [0, ∞) sind messbar. Dies folgt aus der Messbarkeit von f(Y,X)
und dem Satz von Fubini (siehe [9, Satz 14.16]). Für A ∈ B(Rs ) gilt
Z Z
f(Y,X) (y, x)λd (dx) λs (dy).
P(Y ∈ A) = P((Y, X) ∈ A × Rd ) =
A Rd
|
{z
=fY (y)
}
Also fY VD von Y .
Analog zeigt man mit dem Satz von Fubini für B ∈ B(Rd ), dass
Z
P(X ∈ B) = fX (x)λd (dx).
B
(ii) Ansatz:
P X|Y (y, B) :=
Z
fX|Y =y (x)λd (dx),
y ∈ Rs , B ∈ B(Rd )
B
X|Y
Zu zeigen: P
aus Def. 2.13.
ist reguläre Version der bedingten Verteilung von X unter Y . D.h. prüfe (i)-(iii)
(i) P X|Y (·, B) ist B(Rs )-messbar: Dies gilt, da fX|Y =y (x) in (y, x) eine B(Rs × Rd ) messbare
Funktion ist (als Komposition messbarer Funktionen). Nun folgt (i) mit Fubini-Argument.
(ii) P X|Y (y, ·) Wmaß:
16
2.3
Bedingte Verteilungsdichten
a) Sicher ist P X|Y (y, ·) Maß, wegen Eigenschaften des Integrals.
R
 Rd f(Y,X) (y,x)λd (dx)
fY (y) 6= 0 = 1
fY (y)
b) P X|Y (y, Rd ) = R

f (x)λd (dx)
fY (y) = 0
Rd X
(da fX VD und Def. von fY ).
(iii) Zu zeigen: ∀B ∈ B(Rd ) : P X|Y (·, B) = P({X ∈ B} | Y = ·) PY -f.s. Idee: Wende Satz 2.8
an. Sei dazu A ∈ B(Rs ). Es gilt
Z
Z Z
P X|Y (y, B)PY (dy) =
fX|Y =y (x)λd (dx)PY (dy)
A
A B
Z Z
fX|Y =y (x)λd (dx) · fY (y)λs (dy)
=
A B
Z
Z
=
f(Y,X) (y, x)
λd (dx)fY (y)λs (dy)
fY (y)
A∩{fY 6=0} B
Z
f(Y,X) (y, x)λs+d (d(y, x))
=
(Fubini)
(A∩{fY 6=0})×B
= P((Y, X) ∈ (A ∩ {fY 6= 0}) × B)
= P((Y, X) ∈ A × B)
Z
=
1X∈B dP
(da P((Y, X) ∈ {fY = 0} × B) = 0)
{Y ∈A}
=⇒ Mit Satz 2.8 folgt (iii).
17
3
Markovketten
3.1
Definitionen und erste Eigenschaften
Seien (Ω, F, P) ein Wahrscheinlichkeitsraum und (G, G) ein messbarer Raum. (G ist eigentlich immer
metrischer, separabler, vollständiger Raum und G = B(G). Zum Beispiel G ⊆ Rd und B(G) bezüglich
der euklidischen Metrik.)
Definition 3.1 (Übergangskern, Markovkern)
Die Funktion K : G × G → [0, 1] heißt Markovkern oder Übergangskern, falls
(i) ∀x ∈ G ist K(x, ·) Wmaß auf (G, G),
(ii) ∀A ∈ G ist K(·, A) eine G-messbare Abbildung.
Interpretation: Sei x1 ∈ G. Dann können wir eine Folge x1 , . . . , xn ∈ G, n ∈ N, erzeugen, indem für
i ≥ 2, i ∈ N, xi bezüglich der Verteilung K(xi−1 , ·) gewählt wird. K beschreibt in diesem Sinne den
Übergang zwischen xi−1 und xi .
Definition 3.2 (Markovkette, Startverteilung)
Eine Folge von Zufallsvariablen (Xn )n∈N mit
Xi : (Ω, F, P) → (G, G)
für i ∈ N heißt Markovkette mit Übergangskern K, falls für alle n ∈ N und A ∈ G gilt:
(i) P(Xn+1 ∈ A | X1 , . . . , Xn ) = P(Xn+1 ∈ A|Xn ) P-f.s.,
(ii) P(Xn+1 ∈ A | Xn = x) = K(x, A) PXn -f.s. in x ∈ G.
[Anders formuliert: P(Xn+1 ∈ A | Xn ) = K(Xn , A) P-f.s.]
Die Verteilung ν(A) = P(X1 ∈ A), für A ∈ G, heißt Startverteilung.
Bemerkung 3.3
a) Die Gleichheit (i) wird Markoveigenschaft (ME) genannt. Intuitiv: Die Verteilung der Zukunft (also von Xn+1 ) hängt nur von der Gegenwart (Xn ) ab und nicht von der Vergangenheit
(X1 , . . . , Xn−1 ).
b) Die Gleichheit in (ii) besagt, dass P Xn+1 |Xn (x, A) = K(x, A), n ∈ N, eine reguläre Version der
bedingten Verteilung von Xn+1 unter Xn ist, und zwar für alle n ∈ N.
c) Eine Folge von Zufallsvariablen (Xn )n∈N , die (i) und (ii) erfüllt, heißt auch homogene Markovkette (da in (ii) K unabhängig von n).
d) Die Indexmenge in (Xn )n∈N sind die natürlichen Zahlen. Solche Markovketten werden auch zeitlich diskret genannt (zeitlich stetig (Indexmenge z.B. R, [0, 1]) heißt Markovprozess).
Frage: Nun sei K ein Übergangskern und ν ein Wahrscheinlichkeitsmaß auf (G, G). Wie können wir
eine Markovkette (Xn )n∈N mit Übergangskern K und Startverteilung ν bekommen?
Definition 3.4 (Update-Funktion)
Sei (E, E, µ) Wraum. Dann heißt eine messbare Abbildung Φ : G × E → G Update-Funktion von
K, falls
P(Φ(x, U ) ∈ A) = K(x, A), ∀x ∈ G, A ∈ G,
18
3.1
Definitionen und erste Eigenschaften
wobei U : (Ω, F, P) → (E, E) eine Zufallsvariable mit Verteilung µ ist (PU = µ).
Nun sei (Un )n∈N eine iid-Folge von Zufallsvariablen mit Ui : (Ω, F, P) → (E, E) und PUi = µ, i ∈ N.
Φ sei eine Update-Funktion von K. Desweiteren sei die ZV X1 : (Ω, F, P) → (G, G) mit PX1 = ν,
unabhängig von (Un )n∈N . Dann ist (Xn )n∈N gegeben durch
n≥2
Xn = Φ(Xn−1 , Un ),
eine Markovkette mit Übergangskern K und Startverteilung ν.
Beweis (Begründung):
Die Verteilung von X1 ist ν nach Definition. Zu (ii):
P(Xn ∈ A|Xn−1 = x) = P(Φ(Xn−1 , Un ) ∈ A|Xn−1 = x)
= E(1A (Φ(Xn−1 , Un ))|Xn−1 = x)
E6
= E[1A (Φ(x, Un ))]
= K(x, A)
(Xn−1 , Un unabh., Abschn. 2.1)
(Φ Update-Funktion).
Zu (i): Intuitiv klar. Man argumentiert ähnlich, wie im Beweis von (ii).
Frage: Wann existiert eine Update-Funktion von K?
Satz 3.5
Sei G = R und G = B(R). Für jeden Übergangskern K auf (G, G) existiert eine Update-Funktion
Φ : G × [0, 1] → G
mit
P(Φ(x, U ) ∈ A) = K(x, A),
∀x ∈ G, A ∈ B(G),
wobei U : (Ω, F, P) → ([0, 1], B([0, 1])) Zufallsvariable und PU ist die Gleichverteilung in [0, 1].
Beweis: Seien x, y ∈ R. Setze F (x, y) := K(x, (−∞, y]) (mit x parametrisierte Verteilungsfunktion).
Eigenschaften von F :
1) F (x, y) = limh→0+ F (x, y + h) (rechtsseitig stetig)
2) ∀y < z gilt F (x, y) ≤ F (x, z) (Monotonie)
Nun sei Φ : G × [0, 1] → G gegeben durch
Φ(x, u) = inf{y ∈ R : F (x, y) ≥ u}
(verallg. Inverse von F (x, ·))
gegeben. Es gilt
Φ(x, u) ≤ y
⇐⇒
F (x, y) ≥ u.
(3)
„⇐“ X
„⇒“ z > y =⇒ Φ(x, u) < z =⇒ F (x, z) ≥ u und
F (x, y) = lim+ F (x, y + h) ≥ u.
{z
}
h→0 |
≥u
19
3.1
Definitionen und erste Eigenschaften
Zur Messbarkeit von Φ: Es genügt zu zeigen, dass Φ−1 ((−∞, r]) ∈ B(R) ⊗ B([0, 1]) für alle r ∈ R gilt.
Dazu sei g : R × [0, 1] → R gegeben durch g(x, u) := F (x, r) − u. Diese Abbildung ist messbar, als
Komposition messbarer Abbildungen (F (·, r) messbar). Wir erhalten
Φ−1 ((−∞, r]) = {(x, u) ∈ R × [0, 1] : Φ(x, u) ≤ r}
(3)
= {(x, u) ∈ R × [0, 1] : F (x, r) ≥ u}
= {(x, u) ∈ R × [0, 1] : g(x, u) ≥ 0}
= g −1 ([0, ∞)) ∈ B(R) ⊗ B([0, 1])
(g messbar)
Nun zeigen wir, dass P(Φ(x, U ) ≤ r) = K(x, (−∞, r]) für alle r ∈ R gilt, d.h. die Verteilungsfunktion
von Φ(x, U ) und K(x, (−∞, r]) stimmen überein. Damit folgt P(Φ(x, U ) ∈ A) = K(x, A), siehe dazu
[4, Thm. 9.1.1 und Thm. 3.26] („Korrespondenzsatz“). Wir haben
(3)
P(Φ(x, U ) ≤ r) = P(F (x, r) ≥ U ) = F (x, r) = K(x, (−∞, r]).
Bemerkung 3.6
Wir haben uns auf den Fall G = R beschränkt. Für G ⊆ Rd oder allgemeine metrische, separable, vollständige Räume gilt das Ergebnis analog. Man argumentiert dann zusätzlich mit einem „MessraumIsomorphismus“, siehe [9, S. 186, Def. 8.53 ff.].
Nun kommen wir zu Beispielen von Übergangskernen.
Beispiel 3.7
Sei stets G ⊆ Rd , A ∈ B(G), x ∈ G.
0.) Sei
(
0
K(x, A) := 1A (x) =
1
x∈
/A
,
x∈A
(beschreibt den Übergang, s.d. man immer in x verbleibt).
Update-Funktion: Φ(x, u) = x (für eine beliebige Menge E und alle u ∈ E).
1.) Sei π Wmaß auf (G, G) = (G, B(G)). Wir nehmen an, dass wir einen Zufallszahlengenerator für
die Verteilung π zur Hand haben. Nun sei K(x, A) = π(A).
(D.h. eine Folge von iid ZVen mit Verteilung π bilden eine Markovkette.)
Update-Funktion: Φ(x, u) = u, wobei u eine Realisierung einer π-verteilten ZV U : (Ω, F, P) →
(G, G) darstellt.
2.) Jede Konvexkombination zweier Übergangskerne ist wieder ein Übergangskern. D.h. für alle λ ∈
[0, 1] ist
K(x, A) = λK1 (x, A) + (1 − λ)K2 (x, A)
ein Übergangskern für Übergangskerne K1 und K2 .
Update-Funktion: Sei Φ1 eine Update-Fkt von K1 und Φ2 eine solche für K2 , wobei Φ1 , Φ2 :
G × E → G. Sei u ∈ E und r ∈ [0, 1] zufällig gewählt, d.h. u ist eine Realisierung einer ZV mit
20
3.1
Definitionen und erste Eigenschaften
Verteilung µ (siehe Def. 3.4) und unabhängig davon sei r eine Realisierung der Gleichverteilung
in [0, 1]. Wir definieren
Φ(x, (u, r)) = 1[0,λ) (r)Φ1 (x, u) + 1[λ,1] (r)Φ2 (x, u).
Algorithmus
Input:
Output:
x ∈ G, Φ1 , Φ2
y∈G
a) Wähle r ∈ [0, 1] gleichverteilt und unabhängig von u bezüglich µ.
b) Falls λ ≤ r, setze y := Φ2 (x, u), sonst y := Φ1 (x, u).
3.) Sei δ > 0 und B(x, δ) = {y ∈ Rd : |x − y| ≤ δ}. Dann ist der Übergangskern des δ ball walk
gegeben durch
λd (A ∩ B(x, δ))
λd (G ∩ B(x, δ))
Wδ (x, A) =
+ 1A (x) · 1 −
.
λd (B(0, δ))
λd (B(0, δ))
Skizze: 1. Fall: B(x1 , δ) ⊆ G, 2. Fall: B(x2 , δ) ∩ Gc 6= ∅.
G
δ
x1
δ
x2
Sei U : (Ω, F, P) → (Bd , B(Bd )) (wobei Bd = B(0, 1)) gleichverteilte ZV in Bd .
Update-Funktion: u sei Realisierung von U (Lemma 1.3 ff.)
(
x + δu x + δu ∈ G
Φ(x, u) =
.
x
x + δu ∈
/G
4.) Sei R ∈ Rd×d eine symmetrische, positiv definite Matrix und v ∈ Rd .
Dann ist die d-dimensionale Normalverteilung mit Erwartungswertvektor v und Kovarianzmatrix
R gegeben durch
Z
1
1
−1
exp − hz − v, R (z − v)i dz, A ∈ B(Rd ),
N (v, R)(A) :=
2
(2π)d/2 |R|d/2
A
wobei |R| = det R und h·, ·i das euklidische Skalarprodukt bezeichnet.
Sei nun ein sogenannter „Gaußian random walk“ Kern auf G = Rd gegeben durch
N (x, A) = N (x, s2 R)(A),
wobei s > 0 ein zusätzlicher Parameter ist.
√
e (x, A) = N ( 1 − s2 x, s2 R)(A) ist denkbar.)
(Auch etwas wie N
21
3.1
Definitionen und erste Eigenschaften
Update-Funktion: Sei U : (Ω, F, P) → (Rd , B(Rd )) standardnormalverteilte ZV, d.h. mit Verteilung N (0, Id ). (Diese kann aus d eindimensionalen standardnormalverteilten ZVen konstruiert
werden.) Sei u ∈ Rd eine Realisierung von U . Dann ist
Φ(x, u) = x + sR1/2 u
eine Update-Funktion von N (Transformationssatz der Normalverteilung).
Es gibt auch unendlich-dimensionale Hilberträume, die mit Gaußmaßen ausgestattet werden
e ) auch denkbar.
können. In diesem Fall ist ein Übergangskern der Form N (oder N
5.) „Hit-and-run“ Übergangskern. Sei G ⊆ Rd beschränkt, Sd−1 Einheitssphäre und θ ∈ Sd−1 , x ∈ G.
Wir definieren
L(x, θ) = {s ∈ R : x + sθ ∈ G} ⊆ R.
(D.h. L(x, θ) ist die Menge der Parameter s, in der die Gerade x + sθ, mit Stützvektor x und
Richtungsvektor θ, das Gebiet G schneidet.)
Dann ist der Übergangskern vom „hit-and-run“-Algorithmus gegeben durch
Z
1
Hθ (x, A)dσ(θ)
H(x, A) =
σ(Sd−1 )
Sd−1
mit
Z
Hθ (x, A) =
1A (x + sθ)
ds
.
λ1 (L(x, θ))
L(x,θ)
Skizze:
x + sθ
G
x+θ
x
x + sθ für s ∈ L(x, θ)
Update-Funktion als Algorithmus:
Input: x ∈ G
Output: y ∈ G
a) Wähle θ ∈ Sd−1 gleichverteilt (siehe Lemma 1.3).
b) Wähle s ∈ L(x, θ) gleichverteilt und setze y = x + sθ.
6.) „Random scan Gibbs sampler“
Sei G ⊆ Rd beschränkt und {e1 , . . . , ed } ⊆ Rd sei die Euklidische Einheitsbasis, d.h.
ei = (0, . . . , 0,
1
↑
i-te Stelle
, 0, . . . , 0).
Dann ist
d
R(x, A) =
22
1X
He (x, A)
d i=1 i
3.1
Definitionen und erste Eigenschaften
ein Übergangskern, wobei die in 5.) eingeführte Notation benutzt wird.
Update-Funktion als Algorithmus:
Input: x ∈ G
Output: y ∈ G
a) Wähle i ∈ {1, . . . , d} bezüglich der Gleichverteilung.
b) Wähle s ∈ L(x, ei ) gleichverteilt und setze y = x + sei .
Bemerkung (zu 5.) und 6.))
Die Gleichverteilung in L(x, θ) für θ ∈ Sd−1 zu erzeugen, ist im Allgemeinen nicht einfach. Für konvexe
Mengen G ⊆ Rd ist es unter Verwendung von Auswertungen von 1G („Membership-Orakel“) mittels
einer eindimensionalen Verwerfungsmethode „vernünftig“ implementierbar.
An den Beispielen haben wir unter anderem gesehen, wie man aus Übergangskernen neue Übergangskerne konstruiert. Beispielsweise ist die Konvexkombination von Übergangskernen wieder ein
Übergangskern. Das ist soetwas wie die „Summe“ von Übergangskernen.
Frage: Was ist mit dem Produkt?
Im Folgenden seien K, K1 , K2 Übergangskerne auf (G, G).
Definition 3.8 (Produkt von Übergangskernen)
Sei x ∈ G, B ∈ G. Dann ist
Z
K1 K2 (x, B) = K2 (y, B)K1 (x, dy)
G
das Produkt von K1 und K2 .
Eigenschaften
E1 K1 K2 ist wieder ein Übergangskern.
(Die Wmaß-Eigenschaft folgt aus K2 und Eigenschaften des Integrals. Zur Messbarkeit siehe [9,
Satz 14.20])
E2 Setze
K 0 (x, A) := 1A (x)
und
K n (x, A) :=
Z
G
K n−1 (y, A)K(x, dy) =
Z
K(y, A)K n−1 (x, dy).
G
E3 Sei (Xn )n∈N eine Markovkette mit Übergangskern K und Startverteilung ν. Dann gilt für alle
x ∈ G, A ∈ G und n, k ∈ N, dass
P(Xk+n ∈ A | Xk = x) = K n (x, A) PXk -f.s..
D.h. K n ist eine reguläre Version der bedingten Verteilung von Xk+n unter Xk .
23
3.1
Definitionen und erste Eigenschaften
Beweis: Induktion über n ∈ N. Für n = 1, siehe Definition einer Markovkette. X Sei nun n ≥ 1.
Es gilt
P(Xk+n+1 ∈ A | Xk ) = E[1{Xk+n+1 ∈A} | Xk ]
= E[E[1{Xk+n+1 ∈A} | X1 , . . . , Xk+n ] | Xk ]
(E3 in 2.1)
= E[P[Xk+n+1 ∈ A | X1 , . . . , Xk+n ] | Xk ]
= E[P[Xk+n+1 ∈ A | Xk+n ] | Xk ]
|
{z
}
(Markoveigenschaft).
K(Xk+n ,A)
Mit dem Faktorisierungslemma folgt
P(Xk+n+1 ∈ A | Xk = x) = E[K(Xk+n , A) | Xk = x].
Nun betrachten wir
E[K(Xk+n , A) | Xk = x]
Satz 2.17
Z
=
K(y, A)P Xk+n |Xk (x, dy).
G
Nach Induktionsvoraussetzung ist P Xk+n |Xk = K n und wir erhalten
Z
E[K(Xk+n , A) | Xk = x] = K(y, A)K n (x, dy) = K n+1 (x, A).
G
E4 Für alle m ∈ N gilt
Z
P(Xm+1 ∈ A) =
K m (y, A)ν(dy).
G
Beweis: Es gilt
P(Xm+1 ∈ A) = E[1{Xm+1 ∈A} ]
= E[E[1{Xm+1 ∈A} | X1 ]]
(E5 in Abschnitt 2.1)
= E[P(Xm+1 ∈ A | X1 )]
= E[K m (X1 , A)] (E3 in Abschnitt 2.1)
Z
= K m (y, A)PX1 (dy)
G
Z
=
K m (y, A)ν(dy).
G
m
E5 Für jedes bezüglich K (x, ·) integrierbare f : G → R gilt
Z
E[f (Xm+1 ) | X1 = x] = f (y)K m (x, dy).
G
24
3.1
Definitionen und erste Eigenschaften
Beweis: Wir haben
E[f (Xm+1 ) | X1 = x]
Satz 2.17
Z
=
E3
f (y)P Xm+1 |X1 (x, dy) =
G
Z
f (y)K m (x, dy).
G
E6 Für alle A, B ∈ G gilt
Z
P(X1 ∈ A, X2 ∈ B) =
K(x, B)ν(dx).
A
Beweis: Es gilt
P(X1 ∈ A, X2 ∈ B) = E[1{X1 ∈A,X2 ∈B} ]
Z Z
=
1A×B (x, y)P X2 |X1 (x, dy)PX1 (dx)
(Satz 2.17 (iii))
G G
Z
K(x, B)ν(dx)
=
(E3).
A
Bemerkung
Nn 3.9
Sei A ∈ i=1 G = G n , n ∈ N. Dann ist
P((X1 , . . . , Xn ) ∈ A)
ein Wahrscheinlichkeitsmaß auf (Gn , G n ) mit
Z
Z
P(X1 ∈ A1 , . . . , Xn ∈ An ) = . . . K(xn−1 , dxn ) . . . K(x1 , dx2 )ν(dx1 ).
A1
An
Für einen Übergangskern K und ein Wmaß ν auf (G, G) benutzen wir im Folgenden für alle m ∈ N
die Notation
Z
m
νP (A) := K m (x, A)ν(dx).
G
(Dies ist das Produkt von K
wieder ein Wmaß auf (G, G).
m
und ν, nach Definition 3.8.) Wegen E4 (oder auch Def. 3.8) ist νP m
Lemma 3.10
Sei f : G → R integrierbar bezüglich νP m . Dann gilt
Z
Z Z
f (y)νP m (dy) =
f (y)K m (x, dy)ν(dx).
G
G G
Beweis (Standardbeweistechnik der Integrationstheorie, Algebraische Induktion):
1. Schritt: f = 1A , für A ∈ G:
Es gilt
νP m (A) =
Z
G
K m (x, A)ν(dx) =
Z Z
1A (y)K(x, dy)ν(dx).X
G G
25
3.1
Definitionen und erste Eigenschaften
2. Schritt: f =
Pn
i=1 ci 1Ai ,
für n ∈ N, ci ≥ 0, Ai ∈ G (also f nichtnegative, einfache Funktion):
Die Behauptung gilt, wegen der Linearität des Integrals und dem 1. Schritt.
3. Schritt: f ≥ 0 messbar: Dann existiert Folge (fn ) von einfachen Funktionen mit fn ↑ f , d.h.
fn ≤ fn+1 punktweise und limn→∞ fn = f .
Nun gilt die Behauptung unter Anwendung des Satzen von der monotonen Konvergenz
(Satz von Beppo Levi1 , siehe [4, Thm. 4.3.2]) und dem 2. Schritt.
Es gilt
Z
f (x)νP m (dx) := lim
Z
n→∞
G
fn (x)νP m (dx)
ZG Z
= lim
n→∞
(Def. aus Integrationstheorie)
fn (y)K m (x, dy)ν(dx)
(2. Schritt)
G G
Z Z
=
f (y)K m (x, dy)ν(dx)
(2-mal Satz von Beppo Levi).
G G
4. Schritt: f integrierbar bezüglich νP m : Dann f = f + − f − , wobei f + (x) = max{f (x), 0} und
f − (x) = max{−f (x), 0}.
Die Behauptung gilt mit dem 3. Schritt und der Linearität des Integrals. D.h.
Z
Z
Z
f (x)νP m (dx) = f + (x)νP m (dx) − f − (x)νP m (dx)
G
G
G
Z Z
+
Z Z
m
f (y)K (x, dy)ν(dx) −
=
G G
f − (y)K m (x, dy)ν(dx)
G G
Z Z
m
f (y)K (x, dy)ν(dx).
=
G G
Bemerkung 3.11
Lemma 3.10 rechtfertigt die Notation
Z
K m (x, dy)ν(dx) = νP m (dy).
G
Definition 3.12 (stationäre Verteilung)
Sei π ein Wmaß auf (G, G). Dann heißt π stationäre Verteilung (von Übergangskern K), falls
Z
π(A) = K(x, A)π(dx), A ∈ G.
G
1g
n
≥ 0 monoton wachsende Folge messbarer Funktionen. Dann gilt
Z
Z
lim gn (x)ν(dx) = lim
n→∞
26
n→∞
gn (x)ν(dx).
3.1
Definitionen und erste Eigenschaften
(D.h. π(A) = πP (A) für alle A ∈ G, also π = πP .)
Interpretation: Sei (Xn )n∈N eine Markovkette mit Übergangskern K und Startverteilung π. Desweiteren sei π eine stationäre Verteilung von K. Dann gilt
P(X1 ∈ A) = π(A)
stationär
=
E4
πP (A) = P(X2 ∈ A).
D.h. falls eine Markovkette eine stationäre Verteilung als Startverteilung besitzt, bleibt die Verteilung
der Xi , für i ∈ N, gleich.
Folgerung 3.13
Sei f : G → R integrierbar bezüglich einer stationären Verteilung π von K. Dann gilt für alle m ∈ N
Z Z
Z
f (y)π(dy) =
f (y)K m (x, dy)π(dx).
G
G G
Bemerkung 3.14
Folgerung 3.13 rechtfertigt die Notation
Z
π(dy) =
K m (x, dy)π(dx)
G
für alle m ∈ N, falls π eine stationäre Verteilung von K ist.
Beispiel 3.15
Seien K1 und K2 Übergangskerne mit stationärer Verteilung π. Dann gilt
a) π ist eine stationäre Verteilung bezüglich λK1 + (1 − λ)K2 für alle λ ∈ [0, 1],
b) π ist eine stationäre Verteilung von K1 K2 und K2 K1 (man benutzt Folgerung 3.13).
Definition 3.16 (Reversibilität)
Sei π eine Verteilung auf (G, G). Dann heißt der Übergangskern K reversibel bezüglich π, falls
Z
Z
K(x, B)π(dx) = K(x, A)π(dx), ∀A, B ∈ G.
A
B
Interpretation: Sei (Xn )n∈N MK mit Übergangskern K und Startverteilung π. K sei reversibel
bezüglich π. Dann gilt
Z
Z
E6
rev.
E6
P(X1 ∈ A, X2 ∈ B) =
K(x, B)π(dx) =
K(x, A)π(dx) = P(X1 ∈ B, X2 ∈ A).
A
B
(Dies ist eine „gewisse“ Symmetrie in A und B.)
Lemma 3.17
(i) Ist K reversibel bezüglich π, dann ist π eine stationäre Verteilung von K.
(ii) Gilt für alle A, B ∈ G mit A ∩ B = ∅
Z
Z
K(x, B)π(dx) = K(x, A)π(dx),
A
B
27
3.1
Definitionen und erste Eigenschaften
folgt bereits, dass K reversibel bezüglich π ist.
(iii) Seien K reversibel bezüglich π, F : G × G → R und m ∈ N. Dann gilt
Z Z
Z Z
F (x, y)K m (x, dy)π(dx) =
F (y, x)K m (x, dy)π(dx),
G G
(4)
G G
falls eines der beiden Integrale existiert.
Beweis:
(i) Wir haben
Z
πP (A) =
rev
K(x, A)π(dx) =
G
A
(ii) Seien A, B ∈ G. Dann folgt
Z
Z
Z
K(x, B)π(dx) =
K(x, B)π(dx) +
A
Z
Z
(K(x, B \ A) + K(x, A ∩ B))π(dx)
B
A∩B
Z
Z
Z
K(x, A \ B)π(dx) +
B
K(x, A ∩ B)π(dx) +
K(x, A ∩ B)π(dx)
A∩B
B\A
Z
Z
K(x, A \ B)π(dx) +
=
=1
K(x, B)π(dx)
K(x, A \ B)π(dx) +
=
K(x, G) π(dx) = π(A).
| {z }
A∩B
A\B
=
Z
B
K(x, A ∩ B)π(dx)
B
Z
K(x, A)π(dx).
=
B
(iii) Wir beschränken uns auf den Fall m = 1. Für m ∈ N folgt die Aussage per Induktion. Wir
wenden wieder die Standardbeweistechnik der Integrationstheorie und einen „kleinen Trick“
an.
0. Schritt Seien A, B ∈ G, F (x, y) = 1A (x) · 1B (y) = 1A×B (x, y). Dann gilt (4) wegen der
Reversibilität, d.h.
Z
Z
K(x, B)π(dx) = K(x, A)π(dx).
A
B
e ∈ G ⊗ G = σ(G × G). Ziel: Zeige, dass (4) für F (x, y) = 1 (x, y) gilt. Wir
1. Schritt Sei A
e
A
benutzen folgenden Satz.
Satz 3.18 (Dynkin-System-Theorem)
Seien Ω 6= ∅ und C, D ∈ P(Ω) = 2Ω . Weiterhin gelte
a) A, B ∈ C =⇒ A ∩ B ∈ C
(D.h. C ist ein π-System),
28
3.1
Definitionen und erste Eigenschaften
b) 1) Ω ∈ D
2) A ∈ D =⇒ Ac ∈ D
3) An ∈ D paarweise disjunkt =⇒
S∞
n=1
An ∈ D
(D.h. D ist ein Dynkin-System.),
c) C ⊆ D.
Dann folgt: σ(C) ⊆ D.
Wir setzen C = G × G = {A × B : A, B ∈ G} und


Z Z
Z Z


D = A ∈ σ(C) :
1A (x, y)K(x, dy)π(dx) =
1A (y, x)K(x, dy)π(dx) .


G G
G G
Bedingung c) aus Satz 3.18 gilt wegen dem 0. Schritt. Weiterhin kann man zeigen,
dass C ein π-System und D ein Dynkin-System ist. Dann folgt σ(C) = D und (4) ist
für F (x, y) = 1Ae(x, y) erfüllt.
C ist ein π-System: A, B ∈ C =⇒ A = A1 ×A2 und B = B1 ×B2 für A1 , A2 , B1 , B2 ∈ G.
Also
A ∩ B = (A1 ∩ B1 ) × (A2 × B2 ) ∈ C.
| {z }
| {z }
∈G
∈G
D ist ein Dynkin-System:
1) G × G ∈ D, da 1 = 1.
2) Sei A ∈ D. Nutze 1Ac = 1 − 1A . Es folgt Ac ∈ D.
3) Seien An ∈ D paarweise disjunkt. Benutze 1S∞
P∞
1An (da disjunkt)
S∞
und Satz von Lebesgue zum Vertauschen von Limes und Integral. Es folgt n=1 An ∈
D.
n=1
An
=
n=1
2. Schritt Für positive, einfache Funktionen F nutze man Schritt 1.
3. Schritt Sei F ≥ 0 messbar. Es existieren Fn ↑ F , einfache Funktionen mit Fn ≥ 0. Dann gilt
mit Satz von B. Levi und 2. Schritt die Behauptung.
4. Schritt F messbar und integrierbar. Zerlege F = F + −F − und Schritt 3 liefert die gewünschte
Aussage.
Bemerkung 3.19
Lemma 3.17 (iii) rechtfertigt die Notation
K(x, dy)π(dx) = K(y, dx)π(dy),
für K reversibel bezüglich π.
Beispiel 3.20
Sei G ⊆ Rd , G = B(Rd ), A ∈ G, x ∈ G.
29
3.1
Definitionen und erste Eigenschaften
0.) Sei K(x, A) = 1A (x). Dann ist K reversibel bezüglich jedem Wmaß ν auf (G, G): Es gilt
Z
Z
K(x, B)ν(dx) = ν(A ∩ B) = K(x, A)ν(dx).
A
B
Insbesondere ist jedes ν eine stationäre Verteilung von K.
1.) Sei π ein Wmaß auf (G, G) und K(x, A) = π(A). Dann ist K reversibel bezüglich π: Es gilt
Z
Z
K(x, B)π(dx) = π(A)π(B) = K(x, A)π(dx).
A
B
2.) Seien K1 und K2 reversibel bezüglich π. Dann gilt
λK1 + (1 − λ)K2
ist reversibel bezüglich π (Beweis leicht). Im Allgemeinen sind aber K1 K2 und K2 K1 nicht wieder
reversibel bezüglich π. Reversibilität gilt jedoch, falls K1 K2 = K2 K1 (sogar gdw., dazu vielleicht
später mehr).
Beweis: Für A, B ∈ G erhalten wir
Z
Z Z
K1 K2 (x, B)π(dx) =
K2 (y, B)K1 (x, dy)π(dx)
A
A G
Z Z
1A (x)K2 (y, B)K1 (x, dy)π(dx)
=
G G
Z Z
1A (y)K2 (x, B)K1 (x, dy)π(dx) (Lemma 3.17 (iii))
=
G G
Z
K2 (x, B)K1 (x, A)π(dx)
=
G
..
.
zurück argumentieren
Z
= K2 K1 (x, A)π(dx).
B
3.) „δ-ball-walk“. Sei δ > 0 und G ⊆ Rd beschränkt. Zur Erinnerung:
λd (G ∩ B(x, δ))
λd (A ∩ B(x, δ))
+ 1A (x) 1 −
Wδ (x, A) =
λd (B(0, δ))
λd (B(0, δ))
x3
x2
x1
30
G
3.2
Grundlegende Algorithmen
Ansatz: Raten (bzw. Intuition) liefert einen Kandidaten für π: Die stationäre Verteilung π ist
gegeben durch die Gleichverteilung
Beweis: Um dies zu beweisen, prüfen wir Reversibilität bezüglich dieser. Seien dazu C, D ∈ B(G)
und C ∩ D = ∅. Dann gilt
Z
Z
λd (D ∩ B(x, δ))
Wδ (x, D)dx =
dx
λd (B(0, δ))
C
C
Z Z
1
=
1B(x,δ) (y)dydx
λd (B(0, δ))
C D
Z Z
1
1B(y,δ) (x)dydx (da 1B(x,δ) (y) = 1B(y,δ) (x))
=
λd (B(0, δ))
C D
Z
= Wδ (x, C)dx (Fubini).
D
Mit Lemma 3.17 (ii) erhalten wir, dass der δ-ball-walk reversibel bezüglich der Gleichverteilung
ist.
3.2
Grundlegende Algorithmen
Sei stets G ⊆ Rd , G = B(G) und ρ : G → (0, ∞) eine nichtnormierte Dichtefunktion, d.h.
Z
0 < ρ(x)dx < ∞.
G
(Es gilt nicht notwendigerweise
R
G
ρ(x)dx = 1.)
Mittels der nichtnormierten Dichtefunktion kann man ein Wmaß πρ mittels
R
ρ(x)dx
πρ (A) = RA
, A∈G
ρ(y)dy
G
definieren.
Ziel: Konstruktion von Übergangskernen, die reversibel bezüglich πρ sind, oder zumindest πρ als
stationäre Verteilung besitzen.
Fragen:
1) Warum soll πρ eine stationäre Verteilung sein?
Intuitiv: Falls πρ eine stationäre Verteilung von Übergangskern K ist, dann gilt πρ = πρ P . Das
heißt πρ ist eine Art „Fixpunkt“ der linearen Abbildung ν 7→ νP (ν Wmaß). Die Idee ist dann, ein
iteratives Verfahren zu starten, um den Fixpunkt zu approximieren. Falls P gewisse Eigenschaften
erfüllt (Kontraktionseigenschaft, wie beim Banachschen Fixpunktsatz gefordert) sollte νP n , für
n groß genug, „nahe“ an πρ liegen.
2) Warum gerade diese Form πρ ?
31
3.2
Grundlegende Algorithmen
a) Sei G beschränkt und ρ(x) = 1G (x). Dann ist πρ die Gleichverteilung in G. Dies umfasst das
Szenario aus Abschnitt 1.
b) statistische Physik: Boltzmannverteilung,
µβ (A) =
1
Zβ
Z
e−βH(x) dx,
A
wobei β die inverse Temperatur, H die Hamiltonfunktion und Zβ die Normierungskonstante
ist.
c) Bayesscher Ansatz in der Statistik:
Gegeben sind Beobachtungen / Daten / Messwerte b ∈ Rs , s ∈ N. Diese werden modelliert, als
Realisierungen von einem zufälligen Vektor
B : (Ω, F, P) → (Rs , B(Rs )),
B = (B1 , . . . , Bs ),
Bi zufällige Größen.
In der klassischen Statistik ist B ein zufälliger Vektor mit Verteilungsdichte
f : Rs × Rd → (0, ∞),
(b, x) 7→ f (b, x)
wobei x ∈ Rd als fester Parameter angenommen wird. Die VD heißt auch Likelihoodfunktion.
Zum Beispiel: Wir machen den Ansatz Bi iid N (m, v 2 ) (1-dimensional) normalverteilt, d.h.
s
Y
(bi −m)2
1
e− 2v2 .
f (b1 , . . . , bs , m, v 2 ) =
2
s/2
| {z }
(2πv )
i=1
=x
Nun bestimmt die Maximum-Likelihood-Methode der klassischen Statistik den Parametervektor
x, sodass f (b, x) maximal wird. Intuitiv hat die Realisierung b die höchste Wahrscheinlichkeit.
Bei dem Bayesschen Ansatz werden B und der Parameter X, hier gegeben durch
X : (Ω, F, P) → (Rd , B(Rd )),
als zufällige Vektoren mit gemeinsamer Verteilungsdichte
p(x) · f (b, x) x ∈ Rd , b ∈ Rs
angesehen. Dabei ist p : Rd → [0, ∞) die W-dichte der sogenannten Priorverteilung von X und
f (·, x) die bedingte Dichte von B unter X = x.
Nun ist die Posteriorverteilung von X unter B = b durch die bedingte Dichte (siehe Satz 2.16)
πX|B=b (x) = R
p(x) · f (b, x)
p(y)f (b, y)dy
Rd
gegeben und von Interesse.
In Kurzschreibweise (wobei ∝ für „proportional“ steht):
posterior
32
∝
prior × likelihood
3.2
Grundlegende Algorithmen
Frage: Warum „Bayes“?
In der Situation
fX|B=b
fB|X=x
fX
fB
bedingte Dichte von X unter der Bedingung b,
bedingte Dichte von B unter der Bedingung x,
Dichte von X,
Dichte von B,
besagt der Satz von Bayes:
fX|B=b (x)
(Elementare Form: P(A | B) =
3.2.1
=
Satz 2.16
fX (x) · fB|X=x (b)
fB (b)
P(A)P(B|A)
.)
P(B)
Metropolis-Hastings Algorithmus
Der Algorithmus wurde von Metropolis et al. in [14] vorgestellt und später von Hastings in [7]
verallgemeinert.
Spezialfall: „Unabhängiger Metropolis sampler“
Annahmen:
• G beschränkt
• Q sei die Gleichverteilung in G, d.h.
Q(A) =
λd (A)
,
λd (G)
A ∈ B(G).
Wir erläutern einen Übergangsmechanismus von xi zu xi+1 (x1 , x2 , . . . ∈ G).
Algorithmus:
Input:
Output:
xi (aktueller Zustand) und Q (Vorschlagskern)
xi+1 (nächster Zustand)
1) Wähle voneinander unabhängig y bezüglich Q und u gleichverteilt in [0, 1].
2) Setze Akzeptanzwahrscheinlichkeit
ρ(y)
α(xi , y) := min 1,
.
ρ(xi )
Falls u ≥ α(xi , y) setze xi+1 := xi , sonst setze xi+1 := y.
3) Ausgabe xi+1 .
Skizze:
33
3.2
Grundlegende Algorithmen
α = 1
ρ
α = 1
α = 3/4
G = [a, b],
y x5 x3
x4
x2
b
=
a
x1
Folge der Realisierung
Zustand in G
x1
y1
x2
y2
x3
y3
x4
y3
x5
.
y4
Warnung: Wir weisen auf den Unterschied zur Verwerfungsmethode zur Konstruktion von Zufallszahlengeneratoren hin (siehe Kapitel 1, Motivation). Bei dieser Verwerfungsmethode werden so lange
Zustände erzeugt, bis einer akzeptiert wird, und letztlich werden nur akzeptierte Zustände ausgegeben. Beim Metropolis-Hastings Algorithmus wird der akzeptierte oder der aktuelle Zustand
zurückgegeben.
Zum allgemeinen Verfahren:
Sei q : G × G → [0, ∞) und für alle x ∈ G sei q(x, ·) integrierbar (bzgl. dem Lebesguemaß) mit
Z
q(x, y)dy ≤ 1.
G
Dann ist

Z

Z
q(x, y)dy + 1A (x) 1 −
Q(x, A) =
A
q(x, y)dy 
G
ein Übergangskern in (G, B(G)).
(In der Literatur wird q auch als Übergangsdichte bezeichnet.)
Wir nennen Q Vorschlagskern.
Frage: Wie können/müssen wir Q modifizieren, sodass ein Übergangskern heraus kommt, der reversibel bezüglich πρ ist? Wir wollen also eine messbare Funktion kρ : G × G → [0, ∞) mit
Z
kρ (x, y)dy ≤ 1
G
angeben, sodass

Z
kρ (x, y)dy + 1A (x) 1 −
Kρ (x, A) =
A
ein bezüglich πρ reversibler Übergangskern ist.
34

Z
kρ (x, y)dy 
G
3.2
Grundlegende Algorithmen
Lemma 3.21
Ein Übergangskern der Form Kρ ist genau dann reversibel bezüglich πρ , wenn
kρ (x, y)ρ(x) = kρ (y, x)ρ(y) f.ü. bzgl. Lebesguemaß.
Beweis: Seien A, B ∈ B(G) beliebig und c =
Z
R
G
ρ(x)dx. Nun gilt:
Z
Kρ (x, B)πρ (dx) =
A
Kρ (x, A)πρ (dx)
B

Z Z
⇐⇒
kρ (x, y)

Z
ρ(x)dy
dx +
c
Z
1 −
A B
A∩B
Z Z
Z
kρ (x, y)dy 
ρ(x)dx
c
G


Z
ρ(x)dy
1 − kρ (x, y)dy  ρ(x)dx
dx +
c
c
B A
A∩B
G
Z Z
Z Z
⇐⇒
kρ (x, y)ρ(x)dydx =
kρ (x, y)ρ(x)dydx
=
kρ (x, y)
A B
B A
Z Z
⇐⇒
Z Z
kρ (x, y)ρ(x)dydx =
A B
kρ (y, x)ρ(y)dydx
(Fubini)
A B
Z Z
⇐⇒
(kρ (x, y)ρ(x) − kρ (y, x)ρ(y)) dydx = 0
A B
„⇐“ Trivial.
„⇒“ Sei f (x, y) = kρ (x, y)ρ(x) − kρ (y, x)ρ(y). Dann gilt für alle A, B ∈ B(G)
Z Z
f (x, y)d(x, y) = 0.
A B
e ∈ B(G × G), dass
Mit dem Satz 3.18 folgt dann für A
Z
f (x, y)d(x, y) = 0.
e
A
e = {(x, y) ∈ G × G : f (x, y) ≥ 0} = {f ≥ 0} ∈ B(G × G). Es gilt
Sei nun A
Z
Z
0=
f dλ =
f + dλ.
f ≥0
G×G
Und da f + ≥ 0 folgt f + = 0 λ-f.ü. Analog zeigt man f − = 0 λ-f.ü.
Dies impliziert f = f + − f − = 0 λ-f.ü. Also kρ (x, y)ρ(x) = kρ (y, x)ρ(y) λ-f.ü. und die Behauptung ist
bewiesen.
35
3.2
Grundlegende Algorithmen
Das heißt, falls für alle x, y ∈ G gilt, dass q(x, y)ρ(x) = q(y, x)ρ(y), wären wir bereits fertig. Dann
würden wir nämlich kρ = q setzen.
Wie sollen wir kρ wählen, s.d. Kρ reversibel bezüglich πρ ist? Wir können o.B.d.A. annehmen, dass
x, y ∈ G exisitieren, mit
q(x, y)ρ(x) > q(y, x)ρ(y),
(5)
denn sonst läge bereits Reversibilität vor.
(Intuitiv: Übergang von x zu y tritt zu häufig auf und Übergang von y zu x zu selten.)
Wir modifizieren die linke Seite von (5): Wir multiplizieren diese mit einer Akzeptanzwahrscheinlichkeit 0 ≤ α(x, y) < 1, s.d. wir Gleichheit erhalten, d.h.
α(x, y)q(x, y)ρ(x) = q(y, x)ρ(y) = α(y, x)q(y, x)ρ(y),
(6)
mit α(y, x) = 1. Wir setzen nun
kρ (x, y) := α(x, y)q(x, y)
und erhalten Reversibilität, da
kρ (x, y)ρ(x) = α(x, y)q(x, y)ρ(x) = q(y, x)ρ(y) = α(y, x)q(y, x)ρ(y) = kρ (y, x)ρ(y).
Wie sieht jetzt α(x, y) aus? Aus der ersten Gleichheit in (6) folgt
α(x, y) =
q(y, x)ρ(y)
<1
q(x, y)ρ(x)
und α(y, x) = 1.
Zusammen ergibt das für alle x, y ∈ G
n
o
(
q(y,x)ρ(y)
min 1, q(x,y)ρ(x)
q(x, y)ρ(x) 6= 0
α(x, y) :=
.
1
sonst
Der Übergangskern ist dann


Z
Z
α(x, y)q(x, y)dy + 1A (x) 1 −
Kρ (x, A) =
A
α(x, y)q(x, y)dy 
G

Z

Z
α(x, y)Q(x, dy) + 1A (x) 1 −
=
A
denn Q(x, dy) = q(x, y)dy + 1dy (x)(1 −
α(x, y)Q(x, dy) ,
G
R
q(x, z)dz) und Einsetzen in rechte Seite liefert


Z
Z
Z
α(x, y)q(x, y)dy + 1A (y) 1 − q(x, z)dz  α(x, y) 1dy (x)
| {z }
=1
A
G
G




Z
Z
Z
+ 1A (x) 1 − α(x, y)q(x, y)dy − 1 − q(x, z)dz  α(x, y) 1dy (x)
| {z }
=1
G
G
G


Z
Z
= α(x, y)q(x, y)dy + 1A (x) 
−q(x, z)dz 
1
A
G




Z
Z
+ 1A (x) 1 − α(x, y)q(x, y)dy  − 1A (x) 
−q(x, z)dz .
1
G
36
G
G
3.2
Grundlegende Algorithmen
Im Prinzip folgt aus der vorherigen Betrachtung (der Wahl von α(·, ·)) folgendes Lemma.
Lemma 3.22
Der Übergangskern Kρ , mit kρ = α(x, y)q(x, y), ist reversibel bezüglich πρ .
Beweis: Da kρ (x, y)ρ(x) = kρ (y, x)ρ(y) und unter Verwendung von Lemma 3.21.
Metropolis-Hastings-Algorithmus: Übergangskern von xi zu xi+1 :
Eingabe:
Ausgabe:
xi (aktueller Zustand)
Q (Vorschlagskern)
xi+1 (nächster Zustand)
1) Wähle unabhängig voneinander y mit Verteilung Q(xi , ·) und u ∈ [0, 1] gleichverteilt.
n
o
q(y,xi )ρ(y)
2) Berechne α(x, y) = min 1, q(x
.
,y)ρ(x
)
i
i
Falls u ≥ α(xi , y), setze xi+1 = xi , sonst setze xi+1 = y.
3) Ausgabe xi+1 .
Bemerkung 3.23
a) Falls q(x, y) = q(y, x), für alle x, y ∈ G, so heißt Kρ Metropolis-Algorithmus oder MetropolisKern.
b) Falls q(x, y) = η(y), für alle x, y ∈ G und eine Funktion η : G → (0, ∞), so heißt Kρ unabhängiger
Metropolis-Algorithmus (siehe Spezialfall am Anfang dieses Abschnittes).
c) Für den Metropolis-Hastings-Algorithmus benötigt man (bezüglich ρ) eigentlich nur Quotienten ρ(x)
ρ(y) . Insbesondere spielt die Normierungskonstante keine Rolle. Funktionsauswertungen von ρ
(bis auf einen konstanten Vorfaktor) reichen vollkommen aus.
Bemerkung 3.24
In der Herleitung und Darstellung von Kρ taucht α : G × G → [0, 1] auf. Um das Ziel „Kρ reversibel
bezüglich πρ “ zu erreichen, kann man α auch durch andere Funktionen, bzw. Akzeptanzwahrscheinlichkeiten α
e : G × G → [0, 1] ersetzen. Zum Beispiel falls q(x, y) = q(y, x), für alle x, y ∈ G, führt die
Wahl
ρ(y)
α
e(x, y) =
ρ(x) + ρ(y)
auch zum Ziel. Der daraus resultierende Übergangsschritt, bzw. Übergangskern, wird Barkers Algorithmus genannt, siehe Barker [3].
Es gilt allerdings α
e(x, y) ≤ α(x, y), für alle x, y ∈ G.
3.2.2
Hit-and-run-Algorithmus
Seien G ⊆ Rd , G = B(G), ρ gegeben wie zu Beginn von Abschnitt 3 und
R
ρ(x)dx
πρ (A) = RA
.
ρ(x)dx
G
(Für ρ(x) = 1G (x) siehe fünftes Beispiel in Bsp 3.7)
37
3.2
Grundlegende Algorithmen
Idee: Übergang von x zu y durch hit-and-run-Schritt.
ρ|{x+sθ∈G:s∈R}
x + sθ
G
x+θ
x
Übergangskern ist
Hρ (x, A) =
Z
1
σ(Sd−1 )
Hθ,ρ (x, A)dσ(θ),
Sd−1
wobei
R∞
Hθ,ρ (x, A) =
−∞
1A (x + sθ)ρ(x + sθ)ds
`ρ (x, θ)
,
mit
Z∞
`ρ (x, θ) =
1G (x + sθ)ρ(x + sθ)ds.
−∞
Wir schreiben auch kurz Hθ für Hθ,ρ und ` für `ρ .
Bemerkung 3.25
Hθ (x, ·) ist die Verteilung, die durch ρ eingeschränkt auf {x + sθ ∈ G : s ∈ R} gegeben ist.
Hit-and-run-Algorithmus: Übergang von xi zu xi+1 durch hit-and-run-Algorithmus:
Input:
Output:
xi (aktueller Zustand)
xi+1 (nächster Zustand)
1) Wähle θ ∈ Sd−1 gleichverteilt.
2) Wähle xi+1 bezüglich Hθ (xi , ·).
Bemerkung 3.26
Schritt 2) ist unter Umständen nicht einfach zu implementieren. Idee für G beschränkt und konvex:
• Bestimme Schnittpunkte von Gerade mit G (falls G konvex erhalte ich 2 Punkte), z.B. mit
Bisektionsverfahren.
• Bestimme globales Maximum ρmax von ρ, eingeschränkt auf dem Segment {x + sθ ∈ G : s ∈ R}.
• Verwerfungsmethode auf 1-dimensionales Problem anwenden. Die auf dem Rechteck gleichverteilten Punkte werden unterhalb des Graphen von ρ akzeptiert und darüber verworfen. Wir
haben die Hoffnung, dass hier die Akzeptanzwahrscheinlichkeit nicht zu klein ist.
Skizze:
38
3.2
Grundlegende Algorithmen
ρmax
ρ|{x+sθ∈G:s∈R}
Es ist nicht mehr so klar, welche Informationen wir von ρ für den Algorithmus benötigen. Reichen
Funktionswerte von ρ aus? Unter weiteren Voraussetzungen kann man dies bejahen, z.B. G konvex,
beschränkt und ρ log-konkav (d.h. log ρ konkav auf G).
Frage: Ist Hρ reversibel bezüglich πρ ?
Lemma 3.27
(i) ∀θ ∈ Sd−1 ist Hθ reversibel bezüglich πρ .
(ii) Hρ ist reversibel bezüglich πρ .
Beweis:
(i) Sei c :=
R
G
ρ(x)dx, A, B ∈ B(G). Dann gilt
Z Z∞
Z
1A (x)1B (x + sθ)ρ(x + sθ)
ρ(x)dsdx
c · `(x, θ)
Hθ (x, B)πρ (dx) =
Rd −∞
Z∞ Z
A
1A (y − sθ)1B (y)ρ(y)ρ(y − sθ)
dyds
c · `(y − sθ, θ)
=
−∞ Rd
Z Z∞
1A (y − sθ)ρ(y − sθ)
dsπρ (dy),
`(y, θ)
=
B −∞
denn `(y − sθ, θ) = `(y, θ), wegen
Z∞
`(y, θ) =
Z∞
1G (y + (t − s)θ)ρ(y + (t − s)θ)dt = `(y − sθ, θ).
1G (y + t̃θ)ρ(y + t̃θ)dt̃ =
−∞
−∞
Mit −s = t folgt
Z Z∞
Z
Hθ (x, B)πρ (dx) =
1A (y + tθ)ρ(y + tθ)
dtπρ (dy)
`(y, θ)
B −∞
A
Z
Hθ (y, A)ρ(dy).
=
B
39
3.2
Grundlegende Algorithmen
(ii) Seien A, B ∈ B(G) beliebig. Dann gilt
Z
Z
Hρ (x, B)πρ (dx) =
A
Z
σ(dθ)
πρ (dx)
σ(Sd−1 )
Hθ (x, B)
A Sd−1
Z
1
σ(Sd−1 )
=
Z
Hθ (x, B)πρ (dx) σ(dθ)
Sd−1 A
|R
B
{z
}
Hθ (x,A)πρ (dx)
..
.
zurückargumentieren
Z
= Hρ (x, A)πρ (dx).
B
3.2.3
Gibbs sampler
Wir unterscheiden 2 verschiedene Algorithmen, die unter dem Begriff „Gibbs sampler“ zusammengefasst werden.
a) Random scan Gibbs sampler
Dieser ist konzeptionell sehr ähnlich zum hit-and-run Algorithmus.
Skizze: d = 2, G ⊆ R2 beschränkt
e2
0
x x + e1
y
Übergangskern
d
Rρ (x, A) =
1X
He (x, A),
d i=1 i
wobei {e1 , . . . , ed } Euklidische Einheitsbasis im Rd ist, d.h.
ei = (0, . . . , 0,
und Hei (x, A) wie beim hit-and-run.
40
1
↑
i-te Stelle
, 0, . . . , 0),
e1
3.2
Grundlegende Algorithmen
Lemma 3.28
Rρ ist reversibel bezüglich πρ .
Beweis: Analog zu Lemma 3.27 (ii).
b) Deterministic scan Gibbs sampler
Dρ (x, A) = He1 He2 . . . Hed (x, A),
das Produkt von Übergangskernen. Skizze: d = 2, Übergang von x zu y:
y
x x + e1
Lemma 3.29 (a) Dρ ist im Allgemeinen nicht reversibel bezüglich πρ .
(b) πρ ist eine stationäre Verteilung von Dρ .
Beweis:
(a) Konstruktion von nicht reversiblem Beispiel. Sei G = [0, 2] × [0, 2] \ [1, 2] × [1, 2] und ρ = 1G ,
d.h. d = 2 und πρ die Gleichverteilung auf G. Setze A1 = [0, 1] × [1, 2] und A2 = [1, 2] × [0, 1].
2
A1
1
A2
0
Wir haben
Z
1
dx
He1 He2 (x, A2 )
6=
λ2 (G)
A1
G
2
Z
He1 He2 (x, A1 )
dx
,
λ2 (G)
A2
denn für alle x ∈ A1 gilt He1 He2 (x, A2 ) = 0 und für alle x ∈ A2 gilt He1 He2 (x, A1 ) =
=⇒ Im Allgemeinen ist Dρ nicht reversibel bezüglich πρ .
1
4.
(b) Aus Lemma 3.27 (i) wissen wir, dass Hej für j = 1, . . . , d reversibel bezüglich πρ ist. Insbesondere ist πρ eine stationäre Verteilung von Hej (siehe Lemma 3.17 (i)). Deshalb gilt
πρ D(A) = πρ He1 He2 . . . Hed (A) = πρ He2 . . . Hed (A) = . . . = πρ (A).
| {z }
πρ
41
3.2
Grundlegende Algorithmen
Bemerkung 3.30
Die „Kosten“ von einem Schritt bezüglich Dρ sind so hoch wie von d Schritten mit Rρ . Ein Vorteil
des „deterministic scan Gibbs sampler“ ist, dass er Produktstrukturen ausnutzen kann. Zum Beispiel:
G = [0, 1]d und ρ = 1G . Bei diesem Beispiel erzeugt ein Schritt von Dρ genau die richtige Verteilung.
3.2.4
Slice sampling
In diesem Abschnitt stellen wir den „simple slice sampler“ und eine Verallgemeinerung, den „hybriden
Slice sampler“ vor.
R
Zur Erinnerung: G ⊆ Rd , ρ : G → (0, ∞) mit 0 < G ρ(x)dx < ∞. Und zusätzlich sei kρk∞ =
supx∈G |ρ(x)| < ∞. Dies ist nicht unbedingt nötig, schließt aber im Wesentlichen nur uninteressante
Beispiele aus.
Für t ∈ [0, kρk∞ ] sei
G(t) = {x ∈ G : ρ(x) ≥ t}
die Levelmenge von ρ bezüglich t. Es sei bemerkt, dass aus
Z∞
Z∞ Z
λd (G(t))dt =
0
Z∞ Z
1G(t) (x)dxdt =
0 G
Z
ρ(x) < ∞
1[0,ρ(x)] (t)dxdt =
0 G
G
folgt, dass λd (G(t)) < ∞ für fast alle t ∈ [0, kρk∞ ] bezüglich dem Lebesguemaß. Für t = 0 gilt dies
nicht unbedingt, aber {0} ist Nullmenge.
a) Simple slice sampler
Idee/Skizze: Übergang von x zu y mit einem Schritt vom simple slice sampler: d = 1
ρ(x)
t
ρ
x
Skizze für d = 2:
42
y
G(t)
3.2
Grundlegende Algorithmen
ρ(x)
t
y
x
G(t)
Übergangskern: Dazu sei
Ut (A) =
λd (A ∩ G(t))
,
λ(G(t))
A ∈ B(G),
die Gleichverteilung auf G(t). Das heißt (Ut )t>0 eine Folge von Übergangskernen/Verteilungen.
Dann ist
1
U (x, A) =
ρ(x)
ρ(x)
Z
Ut (A)dt
0
der Übergangskern vom simple slice sampler.
Algorithmus:
Eingabe:
Ausgabe:
xi (aktueller Zustand)
xi+1 (nächster Zustand)
1) Wähle t ∈ [0, ρ(xi )] bezüglich Gleichverteilung.
2) Wähle xi+1 gleichverteilt in G(t).
Lemma 3.31
U ist reversibel bezüglich πρ .
43
3.2
Grundlegende Algorithmen
Beweis: Seien A, B ∈ B(G) beliebig und c =
Z
R
ρ(x)dx
=
U (x, A)
c
G
ρ(x)dx. Dann gilt
Z ρ(x)
Z
dx
Ut (A)dt
c
B 0
Z∞ Z
B
=
0 B
dxdt
1[0,ρ(x)] (t) Ut (A)
c
| {z }
1G(t) (x)
Z∞
Ut (A)Ut (B)
=
λd (G(t))dt
c
0
..
.
zurück argumentieren
Z
= U (x, B)πρ (dx).
A
Bemerkung 3.32
Schritt 2) des Algorithmus ist problematisch, da das Erzeugen/ die Simulation der Gleichverteilung
in einer d-dimensionalen Menge (zum Beispiel der Levelmenge G(t)) im Allgemeinen nicht in
vernünftiger Zeit realisierbar ist.
Diese Bemerkung stellt die Hauptmotivation für die folgende Verallgemeinerung, bzw. Modifikation, dar.
b) Hybride slice sampler
Sei (Qt )t>0 eine Folge von Übergangskernen, wobei Qt ein Übergangskern auf G(t) ⊆ Rd definiert.
Wir setzen Qt auf (G, B(G)) wie folgt für x ∈ G, A ∈ B(G) fort:
(
0
x∈
/ G(t)
Qt (x, A) :=
.
Qt (x, A ∩ G(t)) x ∈ G(t)
Allerdings schreiben wir für Qt einfach Qt .
Die Idee ist 2) im Algorithmus vom simple slice sampler durch einen Schritt mittels Qt zu ersetzen.
Als Übergangskern ergibt sich für den hybriden slice sampler (basierend auf Qt )
1
Q(x, A) =
ρ(x)
ρ(x)
Z
Qt (x, A)dt,
x ∈ G, A ∈ B(G).
0
Frage: Unter welchen Voraussetzungen an Qt ist der Übergangskern Q reversibel bezüglich πρ ?
Zur Vereinfachung der Notation definieren wir die Dichtefunktion
λd (G(s))
w(s) = R ∞
,
λd (G(t))dt
0
auf ([0, ∞), B([0, ∞))).
44
s ∈ [0, ∞),
3.3
Markovoperator
Lemma 3.33 (siehe auch [10])
Q ist genau dann reversibel bezüglich πρ , wenn
Z∞ Z
Z∞ Z
Qt (x, A)Ut (dx)w(t)dt =
0 B
Qt (x, B)Ut (dx)w(t)dt
0 A
für alle A, B ∈ B(G).
Insbesondere: Falls Qt reversibel bezüglich Ut ist (für fast alle t bezüglich w), dann ist Q reversibel bezüglich πρ .
Wir bemerken: Die Gleichung aus Lemma 3.33, welche die Reversibilität bezüglich πρ charakterisiert, kann als Reversibilitätsbedingung von Qt bezüglich Ut in Erwartung interpretiert werden,
da
Ew (Q(x, dy)U (dx)) = Ew (Q(y, dx)U (dy)),
∀x, y ∈ G.
Beweis (des Lemmas):
R
R∞
Wir wissen bereits, dass c := G ρ(x)dx = 0 λd (G(t))dt (siehe Anfang dieses Abschnittes). Seien
A, B ∈ B(G) beliebig. Wir erhalten
Z
Z ρ(x)
Z
dx
Q(x, A)πρ (dx) =
Qt (x, A)dt
c
B
B
0
Z Z∞
=
dx
1G(t) (x) Qt (x, A)dt
c
| {z }
B 0 =1[0,ρ(x)] (t)
Z∞ Z
Qt (x, A)
=
0 B
dx
λd (G(t))
dt
λd (G(t)) | {z
c }
=w(t)
Z∞ Z
Qt (x, A)Ut (dx)w(t)dt.
=
0 B
Analog gilt
Z∞ Z
Z
Q(x, B)πρ (dx) =
A
Qt (x, B)Ut (dx)w(t)dt.
0 A
Damit folgt die Behauptung.
3.3
Markovoperator
Sei (G, G) ein messbarer Raum und π Wmaß auf (G, G). Für p ∈ [1, ∞) sei


Z
1/p


Lp = f : G → R messbar : kf kp :=
|f (x)|p π(dx)
<∞


G
45
3.3
Markovoperator
und für p = ∞ sei
L∞ = {f : G → R messbar : kf k∞ := ess supx∈G |f (x)| < ∞}.2
Die Funktionenräume Lp sind gegeben durch den Raum der Äquivalenzklassen [f ] =: f der π-f.ü.
gleichen Funktionen mit Repräsentanten der Norm kf kp < ∞. Lp ist ein Banachraum für p ∈ [1, ∞],
d.h. ein normierter, vollständiger Raum.
Wichtige Ungleichungen:
• Jensensche Ungleichung
Sei f ∈ L1 und ϕ : [0, ∞) → R konvex, d.h.
∀λ ∈ [0, 1] ∀x, y ∈ [0, ∞) : ϕ(λx + (1 − λ)y) ≤ λϕ(x) + (1 − λ)ϕ(y).
Dann gilt


Z
Z
ϕ  |f |dπ  ≤ ϕ(|f |)dπ.
G
G
• Höldersche Ungleichung
Seien p, q ∈ [1, ∞] mit
1
p
+
1
q
= 1 und f ∈ Lp , g ∈ Lq . Dann gilt f · g ∈ L1 und
kf · gk1 ≤ kf kp kgkq .
Für den Beweis der Jensenschen Ungleichung siehe [9, Satz 7.9] und für den Beweis der Hölderschen
Ungleichung siehe [9, Satz 7.16].
Folgerung 3.34
Sei p ≤ q und p, q ∈ [1, ∞]. Dann gilt kf kp ≤ kf kq und damit Lq ⊆ Lp .
Beweis: Wir haben
1/q
1/p·q/q 

Z
Z
≤  |f |p·q/p dπ 
= kf kq .
kf kp =  |f |p dπ 
G
G
Nun sei K Übergangskern auf (G, G) und π eine stationäre Verteilung von K.
Definition 3.35 (Markovoperator)
Sei f ∈ Lp . Dann ist der Markovoperator P gegeben durch
Z
P f (x) = f (y)K(x, dy).
G
2 Zur
Erinnerung:
ess supx∈G |f (x)| =
46
inf
sup |f (x)|.
N ∈G
G\N
π(N )=0
3.3
Markovoperator
Interpretation: Sei (Xn )n∈N Markovkette mit Übergangskern K. Dann gilt
E[f (X2 ) | X1 = x]
E5 aus 3.1
Z
f (y)K(x, dy).
=
G
Eigenschaften:
E0 P ist ein linearer Operator.
E1 Ist f (x) = 1 für alle x ∈ G, folgt P f (x) = f (x) = 1.
(Gilt auch allgemeiner für konstante Funktionen.)
R
E2 f ∈ L1 =⇒ P(|f |)(x) = G |f (y)|K(x, dy) ist π-f.s. überall endlich.
Beweis: Nach Folgerung 3.13 ist
Z Z
Z
|f (y)|K(x, dy)π(dx) =
G G
|f (y)|π(dy) < ∞.
G
Wegen der Positivität des Integranden folgt die Behauptung.
E3 P : Lp → Lp und k P kLp →Lp = 1 für p ∈ [1, ∞].
Beweis: Sei f ∈ Lp . Es gilt
k P f kpp =
Z
| P f (x)|p π(dx)
G
=
p
Z Z
f (y)K(x, dy) π(dx)
G
G
Z Z
≤
=
G G
kf kpp
|f (y)|p K(x, dy)π(dx)
(E2 und Jensen, da p ≥ 1)
(Folg. 3.13)
Also haben wir k P f kp ≤ kf kp und k P 1kp = 1. k P kLp →Lp ist gegeben durch
k P kLp →Lp := sup k P f kp .
kf kp ≤1
Also
1 ≤ k P kLp →Lp ≤ 1
und somit k P kLp →Lp = 1.
R
E4 f ∈ Lp =⇒ Pm f (x) = G f (y)K m (x, dy), für m ∈ N.
47
3.3
Markovoperator
Beweis: Induktion über m ∈ N. m = 1 X nach Definition. Weiterhin gilt
Pm+1 f (x) = P(Pm f )(x)


Z
= P  f (z)K m (·, dz) (x)
G
Z Z
=
f (z)K m (y, dz)K(x, dy)
G G
Z
=
f (z)KK m (x, dz)
(Lemma 3.10 mit ν = K)
G
Z
=
f (z)K m+1 (x, dz).
G
E5 K reversibel bezüglich π ⇐⇒ P selbstadjungiert auf L2 .
Bemerkung 3.36
L2 ist ein Hilbertraum mit Skalarprodukt
Z
hf, gi := f · g dπ,
f, g ∈ L2
G
und P heißt selbstadjungiert, falls hP f, gi = hf, P gi.
Beweis (von E5):
„⇐“ klar, da für A, B ∈ G gilt
Z
Z
K(x, B)π(dx) = h1A , P 1B i = hP 1A , 1B i = K(x, A)π(dx).
A
B
„⇒“ Es gilt für beliebige f, g ∈ L2
Z Z
hP f, gi =
f (y)g(x)K(x, dy)π(dx)
G G
Z Z
f (x)g(y)K(x, dy)π(dx)
=
(Lemma 3.17 (iii))
G G
= hf, P gi.
Zur Erinnerung: Sei µ Wmaß auf (G, G). Dann ist
Z
m
µ P (A) := K m (x, A)µ(dx),
G
auch ein Wmaß auf (G, G) (siehe E4 in Abschnitt 3.1).
48
A∈G
3.4
Klassische Konvergenzeigenschaften von Markovketten
Bemerkung 3.37
Diese „Art von Operation“ lässt sich auch auf signierte Maße µ anwenden, d.h. auf σ-additive Funktionen µ : G → R mit µ(∅) = 0. Man kann dann auch zeigen, dass µ 7→ µ P einen linearen, beschränkten
dµ
dµ
∈ Lp induziert ( dπ
verallgemeinerte Dichte von µ und π).
Operator auf den signierten Maßen mit dπ
Lemma 3.38
dµ
∈ L2 (d.h. µ besitzt eine Dichte bezüglich π und ∀A ∈ G gilt
Sei µ ein Wmaß (G, G) mit dπ
R dµ
dµ
µ(A) = A dπ dπ, sowie dπ ≥ 0).
Sei K reversibel bezüglich π. Dann gilt
d(µ Pm )
= Pm
dπ
dµ
dπ
π-f.s.
Beweis: Wir zeigen
m
Z
µ P (A) =
m
P
dµ
dπ
dπ,
A ∈ G.
A
Dann gilt mit dem Satz von Radon-Nycodim, siehe [9, Korollar 7.34], dass
d Pm
dπ
= Pm
dµ
dπ
π-f.s.
Wir haben
µ Pm (A) =
Z
K m (x, A)µ(dx)
G
Z Z
=
1A (y)K m (x, dy)
dµ
(x)π(dx)
dπ
G G
Z Z
dµ
(y)K m (x, dy)π(dx)
dπ
G G
Z
dµ
m
(x)π(dx).
= P
dπ
1A (x)
=
(Lemma 3.17 (iii))
A
3.4
Klassische Konvergenzeigenschaften von Markovketten
Sei (Xn )n∈N eine Markovkette (auf messbaren Raum (G, G)) mit Übergangskern K und Startverteilung ν. Desweiteren sei K reversibel bezüglich π.
n→∞
Frage: Wann gilt P(Xn ∈ A) −→ π(A), für beliebiges A ∈ G?
| {z }
=ν Pn−1 (A)
Definition 3.39
Der totale Variationsabstand von 2 Wmaßen µ, ν auf (G, G) ist gegeben durch
kµ − νktv = sup |ν(A) − µ(A)|.
A∈G
(m(µ, ν) = kµ − νktv ist in der Tat eine Metrik auf der Menge der Wmaße.)
49
3.4
Klassische Konvergenzeigenschaften von Markovketten
Lemma 3.40
dµ
Seien µ, ν Wmaße auf (G, G). Wir nehmen an, dass die Dichtefunktionen dπ
(Dichte von µ bezüglich
dν
π) und dπ (Dichte von ν bezüglich π) existieren. Dann gilt
Z
1 dν
1
dµ .
kµ − νktv =
sup f (x)(µ(dx) − ν(dx)) = −
2 kf k∞ ≤1
2 dπ dπ 1
G
Beweis: Wir zeigen zunächst kµ − νktv =
1
2
R
supkf k∞ ≤1 G f (x)(µ(dx) − ν(dx)). Es gilt
kµ − νktv = sup |µ(B) − ν(B)|
B∈G
Z dµ
dν
= sup (x) −
(x) π(dx)
dπ
dπ
B∈G B
Z
dµ
dν
= sup 1B (x)
(x) −
(x) π(dx) .
dπ
dπ
B∈G G
|
{z
}
(∗)
dµ
dν
Für welches B wird (∗) am größten? Antwort: Der Integrand dπ
− dπ
sollte immer positiv oder negativ
sein. D.h. für
dµ
dν
A := x ∈ G :
(x) ≥
(x)
dπ
dπ
oder Ac wird obiges Integral betragsmäßig am größten und
kµ − νktv = max{µ(A) − ν(A), ν(Ac ) − µ(Ac )} = µ(A) − ν(A).
Andererseits ist aber
Z
Z
1
dν
1
dµ
sup f (x)(ν(dx) − µ(dx) =
sup f (x)
(x) −
(x) π(dx) .
2 kf k∞ ≤1 dπ
dπ
2 kf k∞ ≤1 G
G
|
{z
}
(?)
Für welche Funktion f wird (?) am größten? Antwort: Mit obiger Menge A für
(
−1 x ∈ A
g(x) :=
.
1
x ∈ Ac
Wir erhalten
Z
Z
1
dµ
1
dν
sup f (dν − dµ) = g(x)
(x) −
(x) π(dx)
2 kf k∞ ≤1 dπ
dπ
2
G
G
1
= |ν(A) − µ(A) − ν(Ac ) + µ(Ac )|
2
= |ν(A) − µ(A)| = µ(A) − ν(A).
50
3.4
Klassische Konvergenzeigenschaften von Markovketten
Damit haben wir die erste Gleichheit und eigentlich auch die zweite aufgrund der Wahl von A. Wir
bieten aber auch eine alternative Argumentation an:
R
dν
dµ Wir zeigen nun supkf k∞ ≤1 | G f (x)(µ(dx) − ν(dx))| = dπ
− dπ
. Es gilt
1
Z dν
dν
dµ
dµ sup f
−
dπ ≤ sup kf k∞ −
dπ dπ .
dπ dπ
kf k∞ ≤1 kf k∞ ≤1
1
G
Weiterhin gilt
Z dν
dν
dµ dµ dπ − dπ = dπ (x) − dπ (x) π(dx)
1
G
Z
dν
dµ
= g(x)
(x) −
(x) π(dx)
dπ
dπ
G
Z dν
dµ
≤ sup f
−
dπ .
dπ dπ
kf k∞ ≤1 G
n
Nun fragen wir nach einer Abschätzung von kν P −πktv (K reversibel bezüglich π, ν Startverteilung).
Es gilt
d(ν Pn )
n
2 kν P −πktv = − 1
(Lemma 3.40)
dπ
1
n dν
(Lemma 3.38)
=
− 1
P
dπ
1
n
dν
=
(P −S) dπ − 1 ,
1
R
dν
wobei S, wie üblich, der Solution Operator ist (S(f ) = G f dπ, S( dπ
− 1) = 0, P 1 = 1.) Damit
erhalten wir folgenden Satz.
Satz 3.41
Sei ν ein Wmaß und K reversibel bezüglich π. Weiter sei
Dann gilt
dν
(i) kν Pn −πktv ≤ 12 kPn −SkL1 →L1 dπ
− 11 ,
| {z }
dν
dπ
∈ L1 und für f ∈ L1 sei S(f ) =
R
G
f dπ.
≤2
n
(ii) kν P −πktv ≤
1
2
n
dν
kP −SkL2 →L2 dπ
− 12 , falls
dν
dπ
∈ L2 .
Beweis: Siehe Vorbetrachtung für Formel (i). Um (ii) zu zeigen wende Folgerung 3.34 an.
Bemerkung 3.42
Falls ν = π, dann gilt kν Pn −πktv = 0, da π Pn = π.
dν
dν
In den Abschätzungen (i) und (ii) kommt dπ
− 11 bzw. dπ
− 12 vor. Diese Größen sind für ν = π
ebenfalls Null. D.h. falls die Markovkette bezüglich einer stationären Verteilung startet, dann sind die
Abschätzungen scharf, also erhalten wir Gleichheit.
51
3.4
Klassische Konvergenzeigenschaften von Markovketten
Sei nun
L02 = {f ∈ L2 : S(f ) = 0}.
Dies ist ein abgeschlossener Unterraum von L2 und somit wieder ein Hilbertraum (gleiches Skalarprodukt und Norm in L2 .) Desweiteren erhalten wir
L2 = L02 ⊕ (L02 )⊥ ,
mit
(L02 )⊥ = {f ∈ L2 : f ≡ c, c ∈ R},
das orthogonale Komplement von L02 , da für f ∈ L2 gilt
f = f − S(f ) + S(f ) .
| {z } | {z }
∈L02
∈(L02 )⊥
Lemma 3.43
Sei π eine stationäre Verteilung von K, n ∈ N. Dann gilt P : L02 → L02 und
kPn kL0 →L0 = kPn −SkL2 →L2 .
2
2
Beweis: Wir zeigen P : L02 → L02 . Sei g ∈ L02 . Zu zeigen: S(P g) = 0. Dies gilt wegen
Z Z
Z
Folg. 3.13
g(y)K(x, dy)π(dx)
=
gdπ = S(g) = 0.
G G
G
Nun zu den Operatornormen: Zunächst bemerken wir
Z
Pn (S(f )) =
S(f )K n (x, dy) = S(f ),
(7)
G
kf −
2
S(f )k2
=
2
kf k2
2
− (S(f )) ≤ kf k22
=⇒
kf − S(f )k2 ≤ 1, falls kf k2 ≤ 1.
(8)
Damit erhalten wir
k Pn −SkL2 →L2 = sup k(Pn −S)f k2
kf k2 ≤1
(7)
= sup k Pn (f − S(f ))k2
kf k2 ≤1
(8)
≤ sup k Pn gk2
kgk2 ≤1
g∈L02
= k Pn kL02 →L02 .
Andererseits
k Pn kL02 →L02 = sup k Pn g − S(g) k ≤ sup k Pn f − S(f )k2 = k Pn −SkL2 →L2 .
|{z}
kf k2 ≤1
kgk2 ≤1
g∈L02
=0
Damit folgt die Behauptung.
Diese Vorbetrachtungen motivieren folgende Konvergenzbegriffe:
52
3.4
Klassische Konvergenzeigenschaften von Markovketten
Definition 3.44 (L1 -exponentielle Konvergenz)
Sei α ∈ [0, 1), M ∈ (0, ∞). Dann heißt der Übergangskern K L1 -exponentiell konvergent mit
(α, M ), falls
k Pn −SkL1 →L1 ≤ αn M, ∀n ∈ N.
Eine Markovkette mit Übergangskern K heißt L1 -exponentiell konvergent, falls α und M existieren,
s.d. K L1 -exponentiell konvergent mit (α, M ) ist.
Definition 3.45 (absolute Spektrallücke)
Ein Übergangskern K und der dazugehörige Markovoperator P haben eine (absolute) Spektrallücke
(1 − β), falls
β = k P kL02 →L02 < 1.
Bemerkung 3.46
a) Falls P eine (absolute) Spektrallücke (1 − β) besitzt, folgt mittels der zwei vorherigen Resultate
(Lemma 3.43 und Satz 3.41)
n
n dν
− 1
kν P −πktv ≤ β .
dπ
2
b) Der Name Spektrallücke kommt von folgender Interpretation: Mit σ(P | L2 ) bezeichnen wir das
Spektrum von P : L2 → L2 . Wir wissen, dass P 1 = 1, d.h. 1 ∈ σ(P | L2 ) ist ein Eigenwert und
(L02 )⊥ der dazugehörige Eigenraum. Desweiteren gilt
k P kL02 →L02 = sup{|λ| : λ ∈ σ(P | L02 )},
siehe dazu [21, Korollar VII.1.2].
Da σ(P | L2 ) = σ(P | L02 ) ∪ {1} wissen wir, dass 1 ein isolierter Punkt im Spektrum von P ist, falls
β = k P kL02 →L02 < 1. D.h. es existiert eine Lücke zur 1, die mit Spektrallücke bezeichnet wird.
Satz 3.47
Seien α ∈ [0, 1) und M ∈ (0, ∞). Sei Übergangskern K reversibel bezüglich π. Dann gilt
k Pn −SkL1 →L1 ≤ αn M =⇒ k P kL02 →L02 ≤ α
Beweis (Idee. Details in [18, Prop. 3.24]):
Es gilt wegen der Dualität von L1 und L∞ (wegen Reversibilität ist der duale Operator von P wieder
P) die Äquivalenz
k Pn −SkL1 →L1 ≤ αn M
⇐⇒
k Pn −SkL∞ →L∞ ≤ αn M.
Der Satz von Riesz-Thorin liefert
k Pn −SkL2 →L2 ≤ αn M.
(9)
(siehe Bennet, Sharpley „Interpolation of operators“ Kap. 4, Corollary 1.8)
Weiter gilt Pn −S = (P −S)n , n ∈ N. Dies zeigen wir per Induktion. n = 1 X.
Für n ≥ 1 gilt dann
(P −S)n+1 = (P −S)(P −S)n = (P −S)(Pn −S) = Pn+1 − P S − S Pn +S 2 .
53
3.4
Klassische Konvergenzeigenschaften von Markovketten
Wegen
P(S(f )) = S(f ),
S(S(f )) = S(f ) und
Z
Z Z
n
n
S(P (f )) = P f (x)π(dx) =
f (y)K n (x, dy)π(dx) = S(f )
G
G G
gilt Pn −S = (P −S)n für alle n ∈ N. Nun folgt die Behauptung mit der Spektralradiusformel
k P kL02 →L02 = k P −SkL2 →L2
=
(Lemma 3.43)
1/n
lim k(P −S)n kL2 →L2
n→∞
(P − S selbstadj., siehe [21, VI.1.7])
(9)
≤ lim (αn M )1/n = α.
n→∞
Wir führen einen weiteren Konvergenzbegriff ein.
Definition 3.48 (gleichmäßig ergodisch)
Seien α ∈ [0, 1) und M ∈ (0, ∞). Dann heißt der Übergangskern gleichmäßig ergodisch mit (α, M ),
falls π-f.ü.
kK n (x, ·) − πktv ≤ αn M, n ∈ N
gilt.
Eine Markovkette heißt gleichmäßig ergodisch, falls Konstanten α und M existieren, sodass der dazugehörige Übergangskern K gleichmäßig ergodisch mit (α, M ) ist.
Unter der Annahme der Reversibilität bezüglich π haben wir folgende Beziehungen zwischen den
verschiedenen Konvergenzbegriffen:
kK n (x, ·) − πktv ≤ αn M
(glm. erg. mit (α, M ))
β = k P kL02 →L02 ≤ α
Satz 3.47
k Pn −SkL∞ →L∞ ≤ 2M αn
(L2 - Spektrallücke ≥ 1 − α)
(
n
k P −SkLp →Lp ≤
54
22/p−1 α2n(p−1)/p
2(p−2)/p α2n/p
p ∈ (1, 2)
p ∈ [2, ∞)
k Pn −SkL1 →L1 ≤ 2M αn
(L1 -exp. konv. mit (α, 2M ))
4
Markovketten Monte-Carlo Methoden
Sei (Xn )n∈N eine Markovkette mit Übergangskern K und Startverteilung ν. Sei K reversibel bezüglich
π.
Ziel: Approximation von
Z
S(f ) =
f (y)π(dy),
G
für integrierbare Funktionen f .
Approximationsschema: Markovketten Monte-Carlo Methode (MCMC)
n
Sn,n0 (f ) =
1X
f (Xj+n0 ),
n j=1
mit n, n0 ∈ N (n0 wird burn-in genannt).
Zunächst wollen wir das „starke Gesetz der großen Zahlen“ für Sn,n0 formulieren und damit zeigen,
dass Sn,n0 wohldefiniert ist.
Satz 4.1 (Starkes Gesetz der großen Zahlen für MCMC)
Sei K π-irreduzibel, d.h.
∀A ∈ G ∀x ∈ G ∃n ∈ N : π(A) > 0 =⇒ K n (x, A) > 0.
Dann gilt für alle f ∈ L1
Px
lim Sn,n0 (f ) = S(f ) = 1,
n→∞
für π-fast alle x ∈ G. Hier zeigt x in Px an, dass die Startverteilung ν = δx (deterministischer Start
in x ∈ G).
Bemerkung 4.2
(i) Einen Beweis dieses Satzes findet man in [15, Thm. 17.17]. Einen einfacheren Zugang bietet die
Arbeit [2] von Asmussen u. Glynn „Harris-Recurrence and MCMC: A simplified approach“.
(ii) Die Bedingung der π-Irreduzibilität besagt, dass man jede Menge A ∈ G von jedem x ∈ G in
endlich vielen Schritten mit der Markovkette erreichen kann, falls π(A) > 0 ist. Im Sinne von
lokalen „random walks“ (z.B. ball walk) ist dies eine Zusammenhangsbedingung an G.
Skizze: δ-ball walk
δ - ball walk nicht
π-irreduzibel
Abstand > δ
55
4.1
Fehlerabschätzungen
wahrscheinlich
„schlechte
Konvergenz“, ist aber
π-irreduzibel
Wir sind allerdings an expliziten Abschätzungen des Fehlers von Sn,n0 interessiert.
4.1
Fehlerabschätzungen
Wir betrachten den Fehlerbegriff aus Definition 1.6 für Sn,n0 . Zur Erinnerung:
eν (Sn,n0 , f ) = Eν,K |Sn,n0 (f ) − S(f )|2
1/2
.
Hier zeigt ν, K in Eν,K an, dass in Sn,n0 eine Markovkette mit Übergangskern K und Startverteilung
ν vorliegt.
Folgendes Lemma wird bei der Fehlerananlyse hilfreich sein.
Lemma 4.3
Sei (Xn )n∈N eine Markovkette mit Übergangskern K und Startverteilung ν, sowie i, j ∈ N mit
j ≥ i.
(a) Dann gilt
Z
Pj−1 (f · Pi−j f )dν(x).
Eν,K [f (Xi )f (Xj )] =
G
(b) Insbesondere gilt
Eπ,K [f (Xi )f (Xj )] = hf, Pi−j f i.
Beweis: (b) folgt aus (a), der Stationarität (siehe Folgerung 3.13) und
Z
Z
Lemma 3.10
Pj−1 (f Pi−j f )dν
=
f (x) Pi−j f (x)ν Pj−1 (dx).
G
G
Es verbleibt also (a) zu zeigen. Es gilt
Z Z
Eν,K [f (Xi )f (Xj )] =
f (x)f (y) PXi |Xj (x, dy)PXj (dx)
(Satz 2.16 (ii))
G G
Z Z
=
f (x)f (y)K i−j (x, dy) ν Pj−1 (dx)
G G
|
Z
=
G
56
{z
f (x) Pi−j f (x)
Pj−1 (f Pi−j f )dν
}
(Lemma 3.10).
(E3, E4 aus Abschnitt 3.1)
4.1
Fehlerabschätzungen
Frage: Wie verhält sich der Fehler, falls die Startverteilung eine stationäre Verteilung ist?
Satz 4.4
Sei (Xn )n∈N eine Markovkette mit Übergangskern K und Startverteilung π. Desweiteren sei π eine
stationäre Verteilung von K und β = k P kL02 →L02 < 1. Dann gilt für f ∈ L2
√
1+β
eπ (Sn,n0 , f ) ≤ √ √
kf k2 .
n 1−β
Lemma 4.5
Sei α ∈ [0, 1). Dann gilt
1+2
n
X
αj ≤
j=1
2
1+α
≤
.
1−α
1−α
Beweis: Mittels der geometrischen Reihe erhalten wir
n
X
1+α
1
1+2
−1 =
.
αj ≤ 1 + 2
1
−
α
1−α
j=1
Beweis (von Satz 4.4):
Sei g(x) = f (x) − S(f ). Dann gilt
eπ (Sn,n0 , f )2 = Eπ,K |Sn,n0 (f ) − S(f )|2
2
X
1 n
= Eπ,K g(Xj+n0 )
n j=1
=
n
n−1
n
1 X
2 X X
2
E
[g(X
E[g(Xj+n0 )g(Xk+n0 )]
)
]
+
π,K
j+n
0
n2 j=1
n2 j=1
k=j+1
n
n−1
n
2 X X
1 X
kgk22 + 2
hg, Pk−j gi (Lemma 4.3 (b))
= 2
n j=1
n j=1
k=j+1
n−1
X n−j
X
1
2
kgk22 + 2
kgk2 k Pl gk2
n
n j=1
l=1


n−1
n
XX
1
2
≤ kgk22 1 +
βl
n
n j=1
l=1
!
n
X
1
2
l
β
≤ kgk2 1 + 2
n
≤
(Hölderungl. und Indexverschiebung)
l=1
1+β 1
≤
kgk22
1−β n
(Lemma 4.5).
Da kgk22 = kf k22 − (S(f ))2 ≤ kf k22 , gilt die Behauptung.
57
4.1
Fehlerabschätzungen
Bemerkung 4.6
0.) n0 können wir Null setzen.
1.) Sei K(x, A) = π(A) für alle x ∈ G und für alle A ∈ G (siehe Beispiel 3.7.1)). Dann betrachten wir
die klassische Monte-Carlo Methode (Erzeuge Xi iid bezüglich π). Es gilt
Z
P f (x) =
f (y)π(dy) = S(f )
G
und daraus folgt
β := k P kL02 →L02
Lemma 3.43
=
k P −SkL2 →L2 = 0.
Wir erhalten
eπ (Sn,n0 , f )2 ≤
1
kf k22
n
wie in Satz 1.8 (Fehlerabschätzung der klassischen Monte-Carlo Methode).
Frage: Wie verhält sich der Fehler eν (Sn,n0 , f ) in Abhängigkeit von der Startverteilung?
Idee: Zerlege den Fehler eν (Sn,n0 , f ) in geeigneter Form. Beispielsweise: Zerlegung in Bias
|E[Sn,n0 (f )] − S(f )|
und Varianz
E|Sn,n0 (f ) − E[Sn,n0 (f )]|2 .
Den Bias kann man gut abschätzen, die Varianz macht aber Probleme.
Wir wollen Satz 4.4 benutzen und schlagen folgende Zerlegung vor:
eν (Sn,n0 , f )2 = eπ (Sn,n0 , f )2 + REST.
Das heißt
REST = eν (Sn,n0 , f )2 − eπ (Sn,n0 , f )2 .
Diesen REST wollen wir nun abschätzen.
Lemma 4.7
Sei f ∈ L2 . Sei (Xn )n∈N eine Markovkette mit Übergangskern K und Startverteilung ν mit
L∞ . Dann gilt
eν (Sn,n0 , f )2 − eπ (Sn,n0 , f )2 =
n
n−1
n
1 X
2 X X
2
L
(g
)
+
Lj+n0 −1 (g Pk−j g),
j+n0 −1
n2 j=1
n2 j=1
k=j+1
dν
wobei g = f − S(f ) und Li (h) = h(Pi −S)h, dπ
− 1i für h ∈ L1 .
58
dν
dπ
∈
4.1
Fehlerabschätzungen
Beweis: Wir haben
2
n
X
1
2
eν (Sn,n0 , f ) = Eν,K g(Xj+n0 )
n j=1
=
n
n−1
n
1 X
2 X X
2
E
[g(X
)
]
+
Eν,K [g(Xj+n0 )g(Xk+n0 )]
ν,K
j+n
0
n2 j=1
n2 j=1
k=j+1
=
1
n2
n Z
X
Pj+n0 −1 (g 2 )dν +
j=1 G
2
n2
n−1
X
Z
n
X
Pj+n0 −1 (gP k−j g)dν,
j=1 k=j+1 G
wobei die letzte Gleichheit aus Lemma 4.3 (a) folgt. Für h ∈ L1 und i ∈ N betrachten wir
Z
Z
Pi h(x)ν(dx) − Pi h(x)π(dx)
G
G
Z
dν
i
= P h(x) (x)π(dx) − Pi h(x)π(dx)
dπ
G
G
Z
dν
i
(x) − 1 π(dx)
= P h(x)
dπ
Z
G
dν
dν
=h(Pi −S)h, dπ
− 1i (da hS(h), dπ
− 1i = 0)
Nun setzen wir einmal h = g 2 ∈ L1 (wegen g ∈ L2 ) und einmal h = g Pk−j g.
Lemma 4.8
dν
Sei h ∈ L1 , dπ
∈ L∞ und K sei L1 -exponentiell konvergent bezüglich (α, M ), für ein α ∈ [0, 1) und
M ∈ (0, ∞). Dann gilt mit i ∈ N, dass
Li (h) ≤ |h(Pi −S)h,
dν
dν
− 1i| ≤ αn M khk1 k dπ
− 1k∞ .
dπ
Beweis: Die rechte Ungleichung folgt mit der Hölder-Ungleichung.
Wir haben jetzt alles beisammen um Fehlerabschätzungen für allgemeine Startverteilungen zu formulieren und zu beweisen.
Satz 4.9 (Fehlerabschätzung)
Sei f ∈ L2 . Sei (Xn )n∈N eine Markovkette mit Übergangskern K und Startverteilung ν. Sei K
reversibel bezüglich π und L1 -exponentiell konvergent bezüglich (α, M ) mit einem α ∈ [0, 1) und
M ∈ (0, ∞).
Weiterhin sei
dν
dπ
∈ L∞ . Dann gilt
|eν (Sn,n0 , f )2 − eπ (Sn,n0 , f )2 | ≤
1+α
dν
− 1k∞ αn0 kf k22 .
M k dπ
− α)2
n2 (1
59
4.1
Fehlerabschätzungen
Insbesondere gilt β ≤ α (siehe Satz 3.47), d.h. es existiert eine Spektrallücke und somit
Satz 4.4
eν (Sn,n0 , f )
≤
(1 + β)1/2
(1 + α)1/2
√
kf
k
+
2
n(1 − α)
n(1 − β)1/2
q
dν
M k dπ
− 1k∞ αn0 kf k2 .
Beweis: Sei g(x) = f (x) − S(f ). Dann gilt
|eν (Sn,n0 , f )2 − eπ (Sn,n0 , f )2 |
≤
n
n
n−1
1 X
2 X X
2
L
(g
)
+
Lj+n0 −1 (g Pk−j g)
j+n0 −1
n2 j=1
n2 j=1
(Lemma 4.7)
k=j+1
≤
n
n−1
n
X X
k−j dν
1 X j+n0 −1 j+n0 −1
dν − 1 kgk22 + 2
g P
g ,
α
M
−
1
α
M
2
2
n j=1
dπ
n
dπ
1
∞
j=1
k=j+1
wobei die letzte Ungleichung mit Lemma 4.8 folgt.
Es gilt
kg Pk−j gk1 ≤ kgk2 k Pk−j gk2 ≤ kgk22 k Pk−j kL02 →L02
Satz 3.47
≤
kgk22 αk−j .
dν
Wir erhalten mit ε0 = αn0 M dπ
− 1∞ somit
|eν (Sn,n0 , f )2 − eπ (Sn,n0 , f )2 |
n
n
n−1
2 X X
1 X
j−1
2
ε0 kgk22 αk−1
ε
α
kgk
+
0
2
n2 j=1
n2 j=1
k=j+1


n
n
n−1
2
X
X
X
ε0 kgk2 
=
αk−1 
αj−1 + 2
n2
j=1
j=1 k=j+1


n
n−1
n
2
X
X
X
kgk2 ε0 
≤
αj−1 + 2
αj−1+l 
n2
j=1
j=1 l=1
!
n
X
kgk22 ε0 1
1+2
αl
≤
n2 1 − α
≤
l=1
kgk22 ε0 1 + α
≤
n2 (1 − α)2
(Lemma 4.5.)
Und wegen kgk2 ≤ kf k2 gilt die Behauptung.
Bemerkung 4.10
1) Der burn-in von n0 Schritten ist sinnvoll um den Einfluss der Startverteilung zu verringern (oder
gar zu eleminieren). Für n hinreichend groß verhält sich eν (Sn,n0 , f ) wie eπ (Sn,n0 , f ).
2) Falls ν = π erhalten wir die Abschätzung aus Satz 4.4. Intuitiv: Je besser die Startverteilung („ν
nah an π“), umso näher ist eν (Sn,n0 , f ) an eπ (Sn,n0 , f ).
3) Falls K(x, A) = π(A), dann gilt mit α = 0 und M = 1 die L1 -exponentielle Konvergenz. Wir
erhalten
0≤β
60
Satz 3.47
≤
α ≤ 0.
4.1
Fehlerabschätzungen
D.h. β = 0 und
kf k2
eν (Sn,n0 , f ) ≤ √
n
für n0 ≥ 1 (die rechte Seite der Ungleichung ist unabhängig von ν). Damit ist Satz 4.9 eine
Verallgemeinerung der Fehlerabschätzung der klassischen Monte-Carlo Methode.
4) „Kritik“: Die Voraussetzung, dass K L1 -exponentiell konvergent ist, ist restriktiv.
Frage: Gibt es eine ähnliche Abschätzung für Markovketten deren Übergangskern/Markovoperator
„nur“ eine Spektrallücke besitzt?
Wir geben folgendes Ergebnis ohne Beweis an.
Satz 4.11 (siehe [18, Thm 3.41])
Sei (Xn )n∈N eine Markovkette mit Übergangskern K und Startverteilung ν. Sei π eine stationäre
dν
Verteilung von K und β = k P kL02 →L02 < 1. Dann gilt für f ∈ L4 und dπ
∈ L2
q
√
√
dν
− 1k2
128
β n0 k dπ
1+β
eν (Sn,n0 , f ) ≤ √ √
kf k4 +
kf k4 .
n(1 − β)
n 1−β
Bemerkung 4.12
a) Bemerkung 4.10 lässt sich sinngemäß auch für Satz 4.11 formulieren.
dν < ∞ folgt eine Fehlerabschätzung von
b) Für f ∈ Lp mit p > 2 gilt: Falls β < 1 und dπ
p/(p−2)
Sn,n0 , die für p → 2 gegen unendlich strebt.
Frage: Wie sollte man n0 wählen?
Die Fehlerabschätzungen aus Satz 4.9 und Satz 4.11 von eν (Sn,n0 , f )2 haben im Wesentlichen für
kf kp ≤ 1 (wobei p ∈ {2, 4}) mit α = β die Form
2
2Cν β n0
+ 2
,
n(1 − β) n (1 − β)2
dν
mit einer Konstante Cν , die von M und/oder von dπ
− 1p abhängt. Desweiteren sei N = n + n0
fest vorgegeben (N entspricht den verfügbaren Gesamtressourcen).
e(n, n0 ) :=
Ziel: Minimiere e(n, n0 ) unter der Nebenbedingung N = n + n0 .
Das ist äquivalent zum Bestimmen einer Minimalstelle von e(N − n0 , n0 ). Es gilt folgendes Ergebnis
(siehe [18, Lemma 2.26]):
Für alle η > 0 und N , sowie Cν groß genug, haben wir eine Minimalstelle nopt von e(N − n0 , n0 ) mit
nopt
Da log β −1 = (1 − β) +
P∞
j=2
(1−β)j
,
j
log Cν
log Cν
∈
, (1 + η)
.
log β −1
log β −1
gilt für β ∈ (0, 1]
log β −1 ≥ 1 − β
61
4.1
Fehlerabschätzungen
und für β nahe bei 1 ist
n0 =
log Cν
1−β
eine gute Wahl für n0 .
Man kann sich auch andere Fehlerbegriffe anschauen. Insbesondere sind Konfidenzabschätzungen von
Sn,n0 von Interesse. Ziel ist, für gegebenen Fehlerparameter ε ∈ (0, 1) und Konfidenzparameter δ ∈
(0, 1) eine Abschätzung der Form
Pν,K (|Sn,n0 (f ) − S(f )| ≥ ε) ≤ δ.
(K)
Man erhält leicht folgenden Zusammenhang zum mittleren quadratischen Fehler.
Satz 4.13
Sei (Xn )n∈N eine Markovkette mit Übergangskern K und Startverteilung ν. Dann gilt für ε ∈ (0, 1)
Pν,K (|Sn,n0 (f ) − S(f )| ≥ ε) ≤
eν (Sn,n0 , f )2
.
ε2
Falls K reversibel bezüglich π und L1 -exponentiell konvergent mit (α, M ), sowie f ∈ L2 und
dν
dπ ∈ L∞ erhalten wir für
dν
log(M k dπ
− 1k∞ )
n0 ≥
log α−1
und
n≥
4kf k22 −2 −1
ε δ ,
1−α
dass (K) erfüllt ist.
Folgendes Lemma ist hilfreich.
Lemma 4.14 (Markov Ungleichung)
Sei h : R → [0, ∞) monoton wachsend und Z Zufallsvariable mit Z : (Ω, F, P) → (R, B(R)). Dann
gilt
h(c)P(Z ≥ c) ≤ E[h(Z)],
für c ∈ R.
Beweis: Es gilt
Z
h(c)P(Z ≥ c) =
h(c)1{Z≥c} (ω)dP(ω)
Ω
Z
≤
h(Z)1{Z≥c} (ω)dP(ω)
(Monotonie von h)
Ω
≤ E[h(Z)].
Beweis (von Satz 4.13):
62
4.1
Fehlerabschätzungen
Wende Lemma 4.14 mit Z = |Sn,n0 (f ) − S(f )| und
(
z2
h(z) =
0
z≥0
sonst
an. Die Abschäztungen von n und n0 für L1 -exponentiell konvergente Markovketten folgen aus Satz
4.9.
Für f ∈ L∞ und die klassische Monte-Carlo Methode erhält man unter Verwendung der Hoeffding-Ungleichung (siehe [17, Abschnitt 3.4] oder Abschnitt 8.2 in „A probabilistic theory of pattern
recognition“ Devoyre et al.) ein „besseres“ Resultat.
Satz 4.15
Sei f ∈ L∞ . Sei (Xn )n∈N eine iid-Folge von Zufallsvariablen mit PXn = π, für n ∈ N. Dann gilt
P(|Sn,0 (f ) − S(f )| ≥ ε) ≤ 2e
−
nε2
2kf k2
∞
.
D.h. für n ≥ 2kf k2∞ log(2δ −1 )ε−2 erhält man (K).
Bemerkung 4.16
Satz 4.15 impliziert ein „besseres“ Resultat als Satz 4.13, da in der unteren Schranke an n, um (K)
zu garantieren, ein log δ −1 statt δ −1 auftaucht.
Frage: Gilt eine ähnliche Abschätzung für allgemeine Sn,n0 , falls f ∈ L∞ .
Satz 4.17 (Glynn und Ormoneit, siehe [6])
Sei (Xn )n∈N eine Markovkette mit Übergangskern K und Startverteilung ν. Für ein α ∈ [0, 1) und
M ∈ (0, ∞). Sei
k Pn −SkL∞ →L∞ ≤ αn M, n ∈ N ∪ {0}
erfüllt. Für f ∈ L∞ gilt dann
Pν,K (|Sn,0 (f ) − S(f )| ≥ ε) ≤ 2 exp −
2M
2
1−α kf k∞ )
2nM 2 kf k2∞
(1 − α)2 (nε −
!
(1 − α)ε nε(1 − α)
≤ 2 exp −
−2
.
M kf k∞ 2M kf k∞
D.h. für
n≥2
M2
4kf k∞ M −1
kf k2∞ ε−2 log(2δ −1 ) +
ε
(1 − α)2
1−α
ist (K) erfüllt.
Folgendes Lemma ist hilfreich, um Satz 4.17 zu beweisen.
Lemma 4.18 (Lemma 8.1 Devoyre et al.)
Sei Z Zufallsvariable, d.h. Z : (Ω, F, P) → (R, B(R)) und EZ = 0, sowie a ≤ Z ≤ b, für a, b ∈ R.
Dann gilt für c > 0
E{ecZ } ≤ e
c2 (b−a)2
8
.
63
4.1
Fehlerabschätzungen
Beweis: Die Funktion ψ(z) = ecz ist konvex, d.h. insbesondere ∀λ ∈ [0, 1] gilt
e(λa+(1−λ)b)c ≤ λeac + (1 − λ)ebc .
Die Wahl von λ =
b−z
b−a
für z ∈ [a, b] liefert
ecz ≤
b − z ac z − a bc
e +
e .
b−a
b−a
Wegen EZ = 0 gilt
E{ecZ } ≤
a bc
b ac
e −
e .
b−a
b−a
a
. Dann erhalten wir
Sei p = − b−a
E{ecZ } ≤ (1 − p)eac + pebc
= (1 − p + pe(b−a)c )eac
= (1 − p + pe(b−a)c )e−pc(b−a)
= eΦ(u) ,
mit u = c(b − a) und
Φ(u) = −pu + log(1 − p + peu ).
Wir erhalten Φ(0) = Φ0 (0) = 0 und
Φ00 (u) =
da
rs
(r+s)2
1
4
≤
e−u (1 − p)p
1
≤ ,
(p + (1 − p)e−u )2
4
für alle r, s ≥ 0 gilt.
Mit dem Satz von Taylor gilt: ∃ξ ∈ [0, u], s.d.
u
Φ(u) = Φ(0) + Φ0 (0)u + Φ00 (ξ) .
2
Daraus folgt
Φ(u) ≤
u2
c2 (b − a)2
=
8
8
und damit ist die Behauptung bewiesen.
Beweis (von Satz 4.17):
Sei g(x) = f (x) − S(x) und h(x) =
P∞
k=0
Pk g(x).
1. Schritt: (Eigenschaften von h)
Lemma 4.19
Es gilt khk∞ ≤
M kf k∞
1−α
und für M ≥ 2 auch
khk∞ ≤ 2kf k∞
64
log M/2
1
+1+
log α−1
1−α
.
4.1
Fehlerabschätzungen
Beweis (von Lemma 4.19):
Wir erhalten wegen k Pk −SkL∞ →L∞ ≤ M αk folgendes:
khk∞ ≤
∞
X
k Pk gk∞ =
k=0
∞
X
k(Pk −S)f k∞ ≤
k=0
∞
X
kf k∞ k Pk −SkL∞ →L∞ ≤
k=0
kf k∞ M
.
1−α
Für M ≥ 2 ist es besser kP k − SkL∞ →L∞ ≤ min{2, αk M } zu benutzen.
Für b
k=
l
log M/2
log α−1
m
gilt αbk M ≤ 2. Damit erhalten wir
khk∞ ≤
∞
X
k Pk −SkL∞ →L∞ kf k∞
k=0


bkX
∞
−1
X
= kf k∞ 
k Pk −SkL∞ →L∞ +
k Pk −SkL∞ →L∞ 
k=0
k=b
k


∞
X
= kf k∞ 2b
k+
αk M 
k=b
k
≤ kf k∞ 2
∞
log M/2
bk M X αk
+
2
+
α
| {z }
log α−1
k=0
≤2
| {z }
≤(1−α)−1
Nun ist nach dem Satz von Lebesgue („Integral und Limes vertauschbar“)
P h(x) =
∞
X
Pk+1 g(x) =
k=0
= g(x) − g(x) +
∞
X
Pk g(x)
k=1
∞
X
Pk g(x) = h(x) − g(x).
k=1
Also g(x) = h(x) − P h(x).
2. Schritt: (Eigenschaften von n · Sn,0 (f ) − n · S(f ))
65
4.1
Fehlerabschätzungen
Es gilt
n · Sn,0 (f ) − n · S(f )
n
X
=
g(Xj )
j=1
=
n
X
h(Xj ) − P h(Xj )
(nach Schritt 1)
j=1
= h(X1 ) +
n
X
(h(Xj ) − P h(Xj−1 )) − P h(Xn ) + h(Xn+1 ) − h(Xn+1 )
j=2
= h(X1 ) − h(Xn+1 ) +
n+1
X
h(Xj ) − P h(Xj−1 ).
j=2
Sei D(Xi−1 , Xi ) = h(Xi ) − P h(Xi−1 ). Nun gilt
n · Sn,0 (f ) − n · S(f ) ≤ 2khk∞ +
n+1
X
D(Xj−1 , Xj ).
(10)
j=2
Und für i = 2, . . . , n + 1 haben wir
E[D(Xi−1 , Xi ) | X1 , . . . , Xi−1 ] = E[D(Xi−1 , Xi ) | Xi−1 ]
= E[h(Xi ) | Xi−1 ] − P h(Xi−1 ) = 0,
sowie
|D(Xi−1 , Xi )| ≤ 2khk∞ .
Somit sind die Voraussetzungen von Lemma 4.18 für Z = D(Xi−1 , Xi ) erfüllt und wir
erhalten
2
2 khk∞
8
E[exp(cD(Xi−1 , Xi )) | X1 , . . . , Xi−1 ] ≤ ec
,
für c > 0.
3. Schritt: (Abschätzung von P(Sn,0 (f ) − S(f ) ≥ ε))
Es gilt
P(Sn,0 (f ) − S(f ) ≥ ε) = P(nc(Sn,0 (f ) − S(f )) ≥ ncε)
≤
66
E[ε(nc[Sn,0 (f ) − S(f )])]
encε
(Lemma 4.14 für h(z) = ez ),
(11)
4.1
Fehlerabschätzungen
für c > 0. Wir betrachten
E[ε(nc[Sn,0 (f ) − S(f )])]
"
n+1
#
X
≤ E exp c
D(Xi−1 , Xi )
(mit (10))
i=2
" "
##
n+1
X
= exp(2ckhk∞ )E E exp c
D(Xi−1 , Xi ) | X1 , . . . , Xn
(E5, Abs. 2.1)
i=2
"
#
X
n
= exp(2ckhk∞ )E exp c
D(Xi−1 , Xi ) · E[exp(cD(Xn , Xn+1 )) | X1 , . . . , Xn ]
{z
}
|
i=2
≤exp(c2 khk2∞ /2), wg. (11)
"
X
#
n
2
2
≤ exp(2ckhk∞ + c khk∞ /2) · E exp c
D(Xi−1 , Xi ) .
i=2
Den letzten Abschätzungsschritt können wir iterativ durchführen und bekommen
E[exp(nc[Sn,0 (f ) − S(f )])] ≤ exp(2ckhk∞ + nc2 khk2∞ /2).
Damit erhalten wir für beliebiges c > 0
P(Sn,0 (f ) − S(f ) ≥ ε) ≤ exp(−ncε + 2ckhk∞ + nc2 khk2∞ /2).
Wir minimieren den Exponenten
Φ(c) = −ncε + 2ckhk∞ + nc2 khk2∞ /2
und erhalten die Minimalstelle
cmin =
nε − 2khk∞
.
nkhk2∞
Einsetzen und Umformen liefert
(nε − 2khk∞ )2
P(Sn,0 (f ) − S(f ) ≥ ε) ≤ exp −
.
2nkhk2∞
(12)
4. Schritt: (Abschätzung von P(|Sn,0 (f ) − S(f )| ≥ ε))
Wir verwenden die Ungleichung (12) für f und −f . Es gilt
P(|Sn,0 (f ) − S(f )| ≥ ε) = P(Sn,0 (f ) − S(f ) ≥ ε oder − Sn,0 (f ) + S(f ) ≥ ε)
≤ P(Sn,0 (f ) − S(f ) ≥ ε) + P(Sn,0 (−f ) − S(−f ) ≥ ε)
(nε − 2khk∞ )2
≤ 2 exp −
.
2nkhk2∞
2
M
Da Ψ(v) = − (nε−2v)
monoton wachsend ist, folgt mit khk∞ ≤ 1−α
kf k∞ (Lemma 4.19) die
2nv 2
erste Ungleichung der Behauptung. Die zweite Ungleichung folgt aus Ψ(v) ≤ − vε ( nε
2v − 2).
Bemerkung 4.20
67
4.2
Hinreichende Bedingungen für explizite Konvergenzabschätzungen
a) Die Abschätzung von Satz 4.17 gilt für beliebige Startverteilungen, aber sie erlaubt keine Feineinstellung durch den burn-in.
b) Die Konvergenzbedingung
k Pn −SkL∞ →L∞ ≤ αn M
ist für reversible Übergangskerne zur L1 -exponentiellen Konvergenz äquivalent (Übersicht Ende
Abschnitt 3.4).
c) Die untere Schranke für n (aus Satz 4.17) um (K) zu garantieren hängt von log δ −1 anstatt von
δ −1 ab. Die Abhängigkeit von M ist im Vergleich zu Satz 4.13 „schlechter“. Diese kann verbessert
werden, indem man die Ungleichung
log M/2
1
+
khk∞ ≤ 2kf k∞ 1 +
1−α
log α−1
M
aus Lemma 4.19, statt khk∞ ≤ kf k∞ 1−α
benutzt.
Als nächstes formulieren wir ein Ergebnis (ohne Beweis) für Markovketten mit L2 -Spektrallücke.
Satz 4.21 (Miasojedow, siehe [16])
Sei (Xn )n∈N eine Markovkette mit Übergangskern K und Startverteilung ν. Sei π eine stationäre
dν
∈ L2 und β := k P kL02 →L02 < 1. Dann gilt
Verteilung von K, f ∈ L∞ , dπ
dν 1 − β ε2
·n .
Pν,K (|Sn,0 (f ) − S(f )| ≥ ε) ≤ 2 exp −
dπ
1 + β 2kf k2
2
D.h. für
n≥
∞
dν −1
1+β
δ
4kf k2∞ ε−2 log 2 dπ 1−β
2
ist (K) erfüllt.
Bemerkung 4.22
dν
In Satz 4.21 ist wieder eine Integrabilitätsbedingung der Dichtefunktion dπ
von nöten. Allerdings kann
hier auch ein „burn-in“ zur Feineinstellung benutzt werden. Man betrachte dazu ν Pn0 , die Verteilung
von Xn0 +1 , als Startverteilung. Falls der Übergangskern reversibel bezüglich π ist, gilt dann
n0
d(ν Pn0 ) ≤ 1 + d(ν P ) − 1
dπ dπ
2
2
n dν
0
=1+
(Lemma 3.38)
−1 P
dπ
2
dν
dν
dν
≤ 1 + β n0 (S( dπ
− 1) = 0, d.h. dπ
− 1 ∈ L02 ).
dπ − 1
2
Es ist schwierig explizite Konvergenzabschätzungen von Markovketten zu bekommen. Im folgenden
lernen wir 2 mögliche Techniken dafür kennen.
4.2
4.2.1
Hinreichende Bedingungen für explizite Konvergenzabschätzungen
Doeblin-Bedingung
Sei K ein Übergangskern und π eine stationäre Verteilung auf (G, G).
68
4.2
Hinreichende Bedingungen für explizite Konvergenzabschätzungen
Satz 4.23
Wir nehmen an, dass ein c > 0 existiert, s.d.
∀x ∈ G ∀A ∈ G :
K(x, A) ≥ c · π(A).
(D)
Dann gilt
k Pn −SkL1 →L1 ≤ 2(1 − c)n .
Bemerkung 4.24
Die Bedingung, die in (D) formuliert ist, wird Doeblin-Bedingung genannt. Offensichtlich ist auch
c ≤ 1.
Wir benötigen folgendes Lemma.
Lemma 4.25
Seien µ, ν Wahrscheinlichkeitsmaße auf (G, G) und für c > 0 gelte
µ(A) ≥ c · ν(A),
Dann exisitert die Dichtefunktion
dν
dµ
∀A ∈ G.
und es gilt
0≤
dν
1
≤
dµ
c
µ − f.s.
Beweis (des Lemmas):
Aus dem Satz von Radon-Nikodym (siehe [9, Korollar 7.34]) folgt, dass eine Dichtefunktion
dν
dν
existiert. Wir zeigen dµ
≤ 1c µ-f.s. Sei dazu A = { dµ
> 1c }. Angenommen µ(A) > 0, dann gilt
Z
ν(A) =
dν
dµ
≥0
1
dν
dµ > µ(A) ≥ ν(A).
dµ
c
A
Also µ(A) = 0 und es folgt die Behauptung.
Beweis (des Satzes 4.23):
Aus Lemma 4.25 folgt, dass ∀x ∈ G eine Dichtefunktion
dπ
1
≤
dK(x, ·)
c
dπ
dK(x,·)
mit
K(x, ·) − f.s.
existiert.
69
4.2
Hinreichende Bedingungen für explizite Konvergenzabschätzungen
Sei g ∈ L1 und S(g) = 0. Dann gilt
Z
| P g(x)| = g(y)K(x, dy)
G
Z
Z
dπ
= g(y)K(x, dy) − c g(y)
(y)K(x, dy) dK(x, ·)
G
G
|
{z
}
=0 da S(g)=0
Z
dπ
= g(y) 1 − c
(y) K(x, dy)
dK(x, ·)
{z
}
|
G
≥0 nach Vor.
dπ
(y) K(x, dy)
dK(x, ·)
G
Z
Z
= |g(y)|K(x, dy) − c |g(y)|π(dy).
Z
≤
|g(y)| 1 − c
G
G
Also folgt
Z
k P gk =
| P g(x)|π(dx)
π stat.
≤
kgk1 − ckgk1 = (1 − c)kgk1 .
G
n
Da S(P g) = 0, falls S(g) = 0 (weil π stationäre Verteilung ist), können wir die letzte Ungleichung
iterativ anwenden. D.h.
k Pn gk1 ≤ (1 − c)k Pn−1 gk1 ≤ . . . ≤ (1 − c)n kgk1 .
Schließlich folgt
k Pn −SkL1 →L1 = sup k(Pn −S)f k1
kf k1 ≤1
= sup k Pn (f − S(f ))k1
kf k1 ≤1
≤ sup (1 − c)n
kf k1 ≤1
kf − S(f )k1
{z
}
|
≤kf k1 +|S(f )|≤2kf k1
n
≤ 2(1 − c) .
4.2.2
Leitfähigkeit und Cheeger-Ungleichung
Sei K Übergangskern der reversibel bezüglich π ist. Desweiteren nehmen wir an, dass der Markovoperator P (von K) positiv semidefinit ist, d.h.
∀f ∈ L2 : hP f, f i ≥ 0.
Wir nutzen folgende Darstellung von k P kL02 →L02 :
70
4.2
Hinreichende Bedingungen für explizite Konvergenzabschätzungen
Lemma 4.26
Sei K reversibel bezüglich π und P positiv semidefinit. Dann gilt
k P kL02 →L02 = sup
g∈L2
g6=0
hP g, gi
.
hg, gi
Beweis: Diese Aussage folgt, da für selbstadjungierte, lineare Operatoren T : H → H (H Hilbertraum)
kT kH→H = sup
f ∈H
f 6=0
hT f, f i
hf, f i
(13)
gilt. Dies ist ein Standardresultat aus der Funktionalanalysis, siehe beispielsweise [21, Satz V.5.7]. Wir
setzen H = L02 und somit folgt die Behauptung wegen der Reversibilität, (13) und hP g, gi ≥ 0.
Nun kommen wir zu dem Begriff der Leitfähigkeit / conductance.
Definition 4.27 (Leitfähigkeit)
Sei K reversibel bezüglich π. Dann heißt
R
ϕ=
A
inf
1
π(A)∈(0, 2 ]
K(x, Ac )π(dx)
h1A , P 1Ac i
=
inf
1
π(A)
π(A)
π(A)∈(0, ]
2
Leitfähigkeit von K. Wir schreiben auch manchmal ϕ(K), statt ϕ.
In der Literatur findet der englische Begriff „conductance“ (= Leitfähigkeit) Anwendung.
Bemerkung (Interpretation)
Sei (Xn )n∈N eine Markovkette mit Übergangskern K und Startverteilung π. Dann ist für π(A) > 0:
R
A
K(x, Ac )π(dx)
= P(X2 ∈ Ac | X1 ∈ A)
π(A)
(da PX1 = π).
Es gibt in der Literatur verschiedene Versionen der Leitfähigkeit. Manche sind äquivalent zu obiger
Definition. Insbesondere unterschiedet sich der Nenner im Quotienten in ϕ und die Bedingung im
Infimum.
Wir benötigen folgendes Lemma:
Lemma 4.28
Sei (G, G, π) ein Wahrscheinlichkeitsraum und h : G → [0, ∞) eine Mengenfunktion mit h(A) =
h(Ac ), für alle A ∈ G. Dann gilt
inf
π(A)∈(0,1/2]
h(A) =
inf
h(A).
π(A)∈(0,1)
71
4.2
Hinreichende Bedingungen für explizite Konvergenzabschätzungen
Beweis: Es gilt
h(A) = inf{h(A) : π(A) ∈ (0, 1/2] ∪ [1/2, 1)}
inf
π(A)∈(0,1)
= min {inf{h(A) : π(A) ∈ (0, 1/2]}, inf{h(A) : π(A) ∈ [1/2, 1)}}
= min {inf{h(A) : π(A) ∈ (0, 1/2]}, inf{h(A) : π(Ac ) ∈ (0, 1/2]}}
= min {inf{h(A) : π(A) ∈ (0, 1/2]}, inf{h(Ac ) : π(A) ∈ (0, 1/2]}}
=
h(A)
inf
(da h(A) = h(Ac )).
π(A)∈(0,1/2]
Folgerung 4.29
Es gilt
ϕ=
inf
π(A)∈(0,1)
h1A , P 1Ac i
.
min{π(A), π(Ac )}
Beweis: Wir haben
ϕ=
inf
π(A)∈(0, 12 ]
h1A , P 1Ac i
π(A)
π(A)=min{π(A),π(Ac )}
=
inf
π(A)∈(0, 12 ]
h1A , P 1Ac i
.
min{π(A), π(Ac )}
Da K reversibel bezüglich π ist, folgt die Selbstadjungiertheit von P (siehe E5 aus Abschnitt 3.3). Für
h(A) =
h1A , P 1Ac i
min{π(A), π(Ac )}
erhalten wir damit h(A) = h(Ac ).
Somit gilt die Behauptung mit Lemma 4.29.
Nun kommen wir zum Hauptergebnis dieses Abschnittes. Es folgt die Beziehung zwischen Leitfähigkeit
und der Spektrallücke 1 − β = 1 − k P kL02 →L02 .
Satz 4.30 (Cheeger-Ungleichung)
Sei K reversibel bezüglich π und der dazugehörige Markovoperator P sei positiv semidefinit. Dann
gilt
ϕ2
≤ 1 − k P kL02 →L02 ≤ 2ϕ.
2
Bemerkung 4.31
Die Ungleichung kann man im Allgemeinen nicht (wesentlich) verbessern. Siehe dazu [11].
Lemma 4.32
Sei g ∈ L2 und π({g 2 > 0}) ≤ 12 . Dann gilt
Z Z
(g(x) − g(y))2 K(x, dy)π(dx) ≥ ϕ2 kgk22 .
G G
Bemerkung 4.33
Sei (Xn )n∈N Markovkette mit Übergangskern K und Startverteilung π (K reversibel bezüglich π).
72
4.2
Hinreichende Bedingungen für explizite Konvergenzabschätzungen
Dann gilt
Eπ,K |g(X1 ) − g(X2 )|2 =
Z Z
(g(x) − g(y))2 K(x, dy)π(dx).
G G
Beweis (von Lemma 4.32):
O.B.d.A. sei g 6= 0, d.h. kgk2 6= 0. Es gilt nach 3. binomischer Formel und Hölder-Ungleichung
(p = q = 2)
Z Z
2
2
2
|g(x) − g(y) |K(x, dy)π(dx)
G G
Z Z
≤
(g(x) − g(y))2 K(x, dy)π(dx) ·
G G
Z Z
G G
(g(x) + g(y))2 K(x, dy)π(dx).
{z
}
|
≤2g(x)2 +2g(y)2
Nebenbetrachtung: Es gilt
Z Z
2
(g(x) + g(y))2 K(x, dy)π(dx) = 2kgk22 + 2hP g 2 , 1i = 4kgk22 .
G G
Damit folgt
Z Z
(g(x) − g(y))2 K(x, dy)π(dx) ≥
G G
1
kgk22
2
|g(x)2 − g(y)2 |K(x, dy)π(dx) .
Z Z
(W)
G G
Nun sei A(t) = {x ∈ G : g(x)2 > t} die Levelmenge von g 2 (für t1 ≤ t2 gilt A(t2 ) ⊆ A(t1 )). Damit
erhalten wir für x, y ∈ G
2
g(x)
Z
2
(g(x) − g(y) )1A(g(y)2 ) (x) =
2
2
g(y)
Z
dt 1A(g(y)2 ) (x)
dt −
0
Z∞
0
(1[0,g(x)2 ) (t) − 1[0,g(y)2 ] (t))dt 1A(g(y)2 ) (x)
=
0
Z∞
=
1[g(y)2 ,g(x)2 ) (t)dt 1A(g(y)2 ) (x)
0
Z∞
1A(t) (x)1A(t)c (y)dt,
=
0
denn es gilt
(
1[g(y)2 ,g(x)2 ) (t) =
1
0
g(y)2 ≤ t, g(x)2 > t
= 1A(t) (x) · 1A(t)c (y).
sonst
Wir betrachten nun den Faktor
Z Z
|g(x)2 − g(y)2 |K(x, dy)π(dx)
G G
73
4.2
Hinreichende Bedingungen für explizite Konvergenzabschätzungen
aus (W). Es gilt
Z Z
|g(x)2 − g(y)2 |K(x, dy)π(dx)
G G
Z Z
=2
(g(x)2 − g(y)2 )1A(g(y)2 ) (x)K(x, dy)π(dx)
(Betrag symm. + Reversibilität)
G G
Z Z Z∞
1A(t) (x)1A(t)c (y)dtK(x, dy)π(dx)
=2
G G 0
Z∞ Z Z
1A(t) (x)1A(t)c (y)K(x, dy)π(dx)dt
=2
(Fubini)
0 G G
Z∞
= 2 h1A(t) , P 1A(t)c idt
0
Zt
=2
∗
h1A(t) , P 1A(t)c i
π(A(t))dt,
π(A(t))
0
wobei t∗ := sup{t : π(A(t)) > 0}. Wegen π(A(t)) ≤ π(A(0)) = π({g 2 > 0}) ≤
Abschätzung:
1
2
gilt folgende
∗
Z Z
|g(x)2 − g(y)2 |K(x, dy)π(dx) ≥ 2ϕ
Zt
π(A(t))dt
0
G G
Z∞ Z
1A(t) (x)π(dx)dt
= 2ϕ
0 G
Z
= 2ϕ
g(x)2 π(dx)
G
= 2ϕkgk22 .
Dies setzen wir in (W) ein und erhalten
Z Z
(g(x) − g(y))2 K(x, dy)π(dx) ≥
G G
2ϕ2 kgk42
= ϕ2 kgk22 .
4kgk22
Damit ist das Lemma bewiesen.
Beweis (von Satz 4.30):
1. Teil: Wir beweisen zunächst
ϕ2
≤ 1 − k P kL02 →L02 .
2
74
4.2
Hinreichende Bedingungen für explizite Konvergenzabschätzungen
Es gilt
ϕ2
≤ 1 − k P kL02 →L02
2
hP f, f i
ϕ2
⇐⇒ sup
≤1−
2
f 6=0 hf, f i
f ∈L02
⇐⇒ ∀f ∈
L02
: hP f, f i ≤
ϕ2
1−
2
hf, f i
⇐⇒ ∀f ∈ L02 : 2h(I − P)f, f i ≥ kf k22 ϕ2 .
Wobei I der Identitätsoperator (mit If = f ) ist.
Sei f ∈ L02 und r (ein) Median von f bezüglich π, d.h.
π({f > r}) ≤
1
2
und π({f < r}) ≤
1
.
2
Wir setzen für alle x ∈ G
f+ (x) := max{f (x) − r, 0},
f− (x) := min{f (x) − r, 0},
es gilt f+ ≥ 0 und f− (x) ≤ 0.
Wir erhalten
2
π({f+
> 0}) = π({f > r}) ≤
1
,
2
2
π({f−
> 0}) = π({f < r}) ≤
1
,
2
da r ein Median ist.
Somit gilt wegen Lemma 4.32
Z Z
(f+ (x) − f+ (y))2 K(x, dy)π(dy) ≥ ϕ2 kf+ k22
(14)
(f− (x) − f− (y))2 K(x, dy)π(dy) ≥ ϕ2 kf− k22 .
(15)
G G
und
Z Z
G G
75
4.2
Hinreichende Bedingungen für explizite Konvergenzabschätzungen
Nun betrachten wir
Z Z
2h(I − P)f, f i =
(f (x) − f (y))2 K(x, dy)π(dx)
(binom. Formel)
G G
Z Z
[ f (x) − r
| {z }
=
G G =f+ (x)+f− (x)
Z Z
=
−( f (y) − r )]2 K(x, dy)π(dx)
| {z }
=f+ (y)+f− (y)
[(f+ (x) − f+ (y)) + (f− (x) − f− (y))]2 K(x, dy)π(dx)
G G
Z Z
≥
ZG ZG
+
(f+ (x) − f+ (y))2 K(x, dy)π(dx)
(f− (x) − f− (y))2 K(x, dy)π(dx)
G G
(wegen 2(f+ (x) − f+ (y))(f− (x) − f− (y)) ≥ 0)
2
≥ ϕ [kf+ k22 + kf− k22 ]
(mit (14) und (15)).
Also
2h(I − P)f, f i ≥ ϕ2 [kf+ k22 + kf− k22 ]

Z
Z

= ϕ
(f (x) − r)2 π(dx) +
{f >r}
= ϕ2
Z


(f (x) − r)π(dx)
{f <r}
(f (x) − r)2 π(dx)
G
= ϕ2 (kf k22 − 2 S(f ) r + r2 )
| {z }
=0
=
≥
2
ϕ (kf k22
ϕ2 kf k22 .
2
+r )
2. Teil: Nun beweisen wir
1 − k P kL02 →L02 ≤ 2ϕ.
Es reicht zu zeigen (Lemma 4.26)
1 − 2ϕ ≤ sup
f ∈L02
f 6=0
hP f, f i
.
hf, f i
f,f i
Idee: Wähle eine einfache Funktion f ∈ L02 , f 6= 0 und berechne hP
hf,f i . Damit erhalten wir
eine untere Abschätzung von k P kL02 →L02 . Es gilt mit f = 1A − π(A) und π(A) ∈ (0, 12 ]
sup
f ∈L02
f 6=0
76
hP(1A − π(A)), 1A − π(A)i
hP f, f i
≥
.
hf, f i
h1A − π(A), 1(A) − π(A)i
4.2
Hinreichende Bedingungen für explizite Konvergenzabschätzungen
Dabei ist
h1A − π(A), 1(A) − π(A)i = π(A) − 2π(A)2 + π(A) = π(A)π(Ac ).
Das heißt
hP f, f i
hP 1A , 1A i − 2hP 1A , π(A)i + π(A)2
≥
hf, f i
π(A)π(Ac )
sup
f ∈L02
f 6=0
hP 1A , 1A i − π(A)2
π(A)π(Ac )
π(A) − hP 1Ac , 1A i − π(A)2
=
(1A = 1 − 1Ac )
π(A)π(Ac )
hP 1Ac , 1A i
=1−
π(A)π(Ac )
hP 1Ac , 1A i
≥1−2
(π(Ac ) ≥ 1/2).
π(A)
=
Wir nehmen das Supremum über alle A mit π(A) ∈ (0, 12 ] über die rechte Seite und erhalten
sup
f ∈L02
f 6=0
hP f, f i
hP 1Ac , 1A i
≥1−2
inf
= 1 − 2ϕ.
1
hf, f i
π(A)
π(A)∈(0, /2]
Somit ist der Satz bewiesen.
Bemerkung 4.34
Ist die Annahme, dass P positiv semidefinit ist, sehr restriktiv?
Falls K reversibel bezüglich π und der dazugehörige Markovoperator P nicht positiv semidefinit ist
(oder man kann die positive Semidefinitheit für P nicht zeigen), dann kann man sich den Übergangskern
1
e
K(x,
A) = (1A (x) + K(x, A)), x ∈ G, A ∈ G
2
e wird „lazy version“ von K genannt. Der Übergangskern K
e ist wieder reversibel bezüglich
anschauen. K
e Letzteres gilt wegen
π und induziert zusätzlich einen positiv semidefiniten Operator P.
Z Z
1
e fi =
hPf,
f (x) + f (y)K(x, dy) f (x)π(dx)
2
G
G
1
hP f, f i
= kf k22 +
2
2
≥ 0 (wg. |hP f, f i| ≤ k P f k2 kf k2 ≤ kf k22 ).
Folgendes Lemma ist später hilfreich.
Lemma 4.35
e
Sei K reversibel bezüglich π und K(x,
A) = 12 (1A (x) + K(x, A)) ist die „lazy version“ von K. Dann
77
4.2
Hinreichende Bedingungen für explizite Konvergenzabschätzungen
gilt
e =
ϕ(K)
1
ϕ(K).
2
e wieder reversibel bezüglich π ist. Es gilt
Beweis: Wir wissen bereits aus Bem. 4.31, dass K
e =
ϕ(K)
e Ac i
h1A , P1
π(A)
π(A)∈(0,1/2]
inf
=0
z }| {
1
1 h1Ac , 1A i +hP 1Ac , 1A i
= ϕ(K)
=
inf
1
π(A)
2
π(A)∈(0, /2] 2
Folgendes Resultat gibt eine hinreichende Bedingung für eine untere Schranke von ϕ.
Satz 4.36
Sei K reversibel bezüglich π und κ > 0. Für beliebiges A ∈ G seien
T1 (A, κ) := {x ∈ A : K(x, Ac ) < κ},
T2 (A, κ) := {y ∈ Ac : K(y, A) < κ},
T3 (A, κ) := G \ (T1 (A, κ) ∪ T2 (A, κ)).
Für ein γ > 0 gelte für alle A ∈ G
π(T3 (A, κ)) ≥ 2γ min{π(T1 (A, κ), π(T2 (A), κ))}.
Dann folgt
ϕ≥
κ
min{1, γ}.
2
Skizze:
G
Ac
A
T1 (A, κ)
T1 (A, κ)
T2 (A, κ)
T3 (A, κ)
T3 (A, κ)
T2 (A, κ)
„Übergang von T1 (A, κ) zu Ac nicht so häufig.“
„Übergang von T2 (A, κ) zu A nicht so häufig.“
„Übergang von A nach Ac und zurück häufig.“
Intuitiv: Die Leitfähigkeit ist „groß“, falls Übergang von A zu Ac (bezüglich K) sehr häufig, für alle
A ∈ G. Die Menge T3 (A, κ) sollte im Vergleich zu T1 (A, κ) und T2 (A, κ) ein großes Maß haben.
78
4.2
Hinreichende Bedingungen für explizite Konvergenzabschätzungen
Beweis: Zur Vereinfachung der Notationen sei Ti := Ti (A, κ), i = 1, 2, 3. Wir machen eine Fallunterscheidung.
1. Fall: π(T1 ) < 21 π(A)
Wir erhalten
Z
Z
K(x, Ac )π(dx) ≥
A
K(x, Ac )π(dx)
A\T1
≥ κ · π(A \ T1 ) = κ(π(A) − π(T1 ))
κ
≥ π(A) (nach Fallunterscheidung)
2
κ
≥ min{π(A), π(Ac )}.
2
2. Fall: π(T2 ) < 12 π(Ac )
Wir erhalten
Z
K(x, Ac )π(dx) =
Z
K(x, A)π(dx)
(Reversibilität)
Ac
A
Z
≥
K(x, A)π(dx)
Ac \T
2
≥ κπ(Ac \ T2 ) = κ(π(Ac ) − π(T2 ))
κ
≥ π(Ac ) (nach Fallunterscheidung)
2
κ
≥ min{π(A), π(Ac )}.
2
3. Fall: π(T1 ) ≥
π(A)
2
Wir erhalten
Z
und π(T2 ) ≥
π(Ac )
2
K(x, Ac )π(dx) =
1
2
A
Z
K(x, Ac )π(dx) +
1
≥
2
K(x, A)π(dx)
(Rev.)
Ac
A
Z
Z
c
Z
K(x, A )π(dx) +
A∩T3
K(x, A)π(dx)
Ac ∩T3
1
κ
(π(A ∩ T3 ) + π(Ac ∩ T3 )) = π(T3 )
2
2
≥ κγ min{π(T1 ), π(T2 )}
κγ
≥
min{π(A), π(Ac )}.
2
≥
Durch Zusammenfügen der Fälle und unter Verwendung von Folgerung 4.30 erhalten wir die
Behauptung.
79
4.3
Anwendungen und Beispiele
4.3
4.3.1
Anwendungen und Beispiele
Integration über konvexe Mengen
Wir betrachten die Fragestellung aus Abschnitt 1.
Sei r ≥ 1, d ∈ N und
K := {G ⊆ Rd : G konvex, Bd ⊆ G ⊆ rBd }.
Sei πG die Gleichverteilung in G ∈ K und Lp = Lp (πG ), für p ∈ [1, ∞]. Wir definieren
FK,p := {(f, G) : G ∈ K, f ∈ Lp , kf kp ≤ 1}.
Ziel: Berechne
R
S(f, G) =
f (x)dx
λd (G)
G
für beliebige (f, G) ∈ FK,p .
Der hit-and-run Algorithmus auf G scheint zur Approximation von πG geeignet zu sein (siehe Beispiel
3.7.5 und oder Abschnitt 3.2.2).
Der Übergangskern ist gegeben durch
H(x, A) =
1
σ(Sd−1 )
Z
∫L(x,θ) 1A (x + sθ)ds
σ(dθ),
λ1 (L(x, θ))
Sd−1
mit
L(x, θ) = {s ∈ R : x + sθ ∈ G}.
Gegebene Informationen:
1) Für beliebiges x erhalten wir von einem Orakel den Funktionswert f (x) auf Anfrage.
e θ) = {x + sθ : s ∈ R}, mit x ∈ rBd und θ ∈ Sd−1 ,
2) Ein weiteres Orakel kann für jede Gerade L(x,
e
die Gleichverteilung in L(x, θ) ∩ G simulieren.
(Diese Annahme passt zur Diskussion in Bemerkung 3.26.)
MCMC Methode: Sei G ∈ K, (Xn )n∈N Markovkette mit Übergangskern H (auf G) und Startverteilung πBd (Gleichverteilung in Bd ). Berechne
n
Sn,n0 (f, G) =
1X
f (Xj+n0 ).
n j=1
Für die Anwendung der Fehlerabschätzungen, stellen sich für beliebiges G ∈ K folgende Fragen:
• Ist H L1 -exponentiell konvergent mit (α, M )? Gibt es explizite Konstanten α und M ?
• Explizite Abschätzungen der Spektrallücke? Oder explizite Abschätzung der Leitfähigkeit?
• Ist H (als Markovoperator) positiv semidefinit in L2 ?
Im Folgenden bezeichnen wir mit H auch den Markovoperator vom hit-and-run Algorithmus. Aus
dem Kontext ist immer klar, ob wir Übergangskern oder Operator meinen.
80
4.3
Anwendungen und Beispiele
Lemma 4.37
Für beliebige G ⊆ Rd gilt
H(x, A) ≥
λd (G)
2
πG (A),
σ(Sd−1 ) diam(G)d
(16)
wobei x ∈ G, A ∈ B(G) und diam(G) = supx,y∈G |x − y|.
Außerdem gilt
λd (Bd )
λd (G)
(2r)λd−1 (Bd−1 )
≤ inf
≤
.
G∈K diam(G)d
(2r)d
(2r)d
(17)
D.h. für alle G ∈ K gilt
≥
H(x, A)
σ(Sd−1 )=dλd (Bd )
und
n
2
πG (A)
d(2r)d
kH − SkL1 →L1
2
≤
2 1−
d(2r)d
Satz 4.23
n
.
Beweis: Die erste Ungleichung in (17) ist sofort klar, da man Zähler und Nenner getrennt abschätzen
kann. Für die zweite Ungleichung konstruieren wir eine Menge G ∈ K, die die Ungleichung erfüllt.
Dazu betrachten wir die folgende Skizze in der Dimension d = 3:
rBd
Bd
r
G
Bd−1
Nun gilt offenbar λd (G) ≤ (2r)λd−1 (Bd−1 ) und diam(G) = 2r. Dieses Argument funktioniert analog
für beliebiges d, womit sich die zweite Ungleichung in (17) ergibt.
Es genügt nun (16) zu zeigen, denn der Rest folgt daraus.
81
4.3
Anwendungen und Beispiele
Für x ∈ G und A ∈ B(G) haben wir
H(x, A) =
1
σ(Sd−1 )
Z∞
Z
Sd−1 −∞
1
=
σ(Sd−1 )
Z
1
=
σ(Sd−1 )
Z
1A (x + sθ)
dsdσ(θ)
λ1 (L(x, θ))
Z0
Sd−1 −∞
Z∞
Z
1A (x + sθ)
dsdσ(θ) +
λ1 (L(x, θ))
Z∞
1A (x + sθ)
dsdσ(θ)
λ1 (L(x, θ))
Sd−1 0
1A (x − sθ)
dsdσ(θ) +
λ1 (L(x, θ))
Sd−1 0
Z∞
Z
1A (x + sθ)
dsdσ(θ) ,
λ1 (L(x, θ))
Sd−1 0
wobei wir die letzte Umformung durch Substitution t = −s im ersten Summanden erhalten. Mit
Polarkoordinatentransformation
Z
Z
f (y)dy
Rd
Z∞
=
y=rθ
y
y=|y| |y|
Sd−1 0
f (rθ)rd−1 drσ(dθ)
erhalten wir weiterhin
1
H(x, A) =
σ(Sd−1 )
Z
1A (x − y)
dy +
λ1 (L(x, y/|y|))|y|d−1
Rd
Z
1A (x + y)
dy .
λ1 (L(x, y/|y|))|y|d−1
Rd
Wir substituieren im ersten Integral x − y = z, d.h. y = x − z, und im zweiten Integral x + y = z, d.h.
y = z − x. Also
Z
Z
1
1A (z)
1A (z)
H(x, A) =
dz
+
dz
.
x−z
z−x
σ(Sd−1 )
λ1 (L(x, |x−z|
λ1 (L(x, |z−x|
))|x − z|d−1
))|z − x|d−1
Rd
Rd
Wegen λ1 (L(x, θ)) = λ1 (L(x, −θ)) folgt schließlich
Z
1A (z)
2
dz
H(x, A) =
x−z
σ(Sd−1 )
λ1 (L(x, |x−z| ))|x − z|d−1
Rd
≥
2
λd (G)
πG (A)
σ(Sd−1 ) diam(G)d
(wegen λ1 (L(x, θ)) ≤ diam(G)).
Aus dem Lemma folgt, dass wir eigentlich jede Fehlerabschätzung für die MCMC Methode anwenden
2
können, da H L1 -exponentiell konvergent bezüglich (α, M ) mit α = 1 − d(2r)
d und M = 2. Allerdings
2
taucht letztlich dann in jeder Abschätzung die Größe 1 − α = d(2r)d im Nenner auf. D.h. Lemma 3.27
impliziert Fehlerabschätzungen, die in der Dimension d exponentiell „schlecht“ sind. Im Wesentlichen
wie bei der Konstruktion von Zufallszahlengeneratoren für πG , siehe Abschnitt 1.
Folgendes Ergebnis ist hilfreich um Satz 4.36 anwenden zu können.
82
4.3
Anwendungen und Beispiele
Lemma 4.38 (siehe [19])
(i) Der Markovoperator H ist positiv semidefinit auf L2 .
(ii) Für alle G ∈ K und A ∈ B(G) beliebig gilt
λd (T3 (A, 10−3 )) ≥
1
min{λd (T1 (A, 10−3 )), λd (T2 (A, 10−3 ))}.
d · 2r · 4000
Beweis: Wir beweisen nur Aussage (i) aus Lemma 4.38. Für die Aussage (ii) verweisen wir auf [12,
Theorem 4.2].
e = Prθ (G) × R, mit
Sei also g : G → R, G ⊆ Rd , θ ∈ Sd−1 , G
Prθ (G) = {x ∈ Rd : x⊥θ, L(x, θ) ∩ G 6= ∅}.
Skizze:
x + sθ
G
x⊥θ
θ
Prθ (G)
e
G
Z
Z
g(y)dy =
G
1G (y)g(y)dy
e
G
Z
Z
1G (x + sθ)g(x + sθ)dsdx
=
Prθ (G) R
nach Substitution y = x + sθ und Satz von Fubini. Damit folgt die Gleichung
Z
Z
Z
g(y)dy =
g(x + sθ)dsdx.
G
(18)
Prθ (G) L(x,θ)
83
4.3
Anwendungen und Beispiele
Für f ∈ L2 bekommt man nun
Z Z
dy
hf, Hf i =
f (y)f (z)H(y, dz)
λd (G)
G G
Z Z
Z
1
f (y + tθ)
=
dr · f (y) dσ(θ)dy
λd (G)σ(Sd−1 )
λ1 (L(y, θ))
G Sd−1 L(y,θ)
{z
|
=
1
λd (G)σ(Sd−1 )
Z
}
=:gθ (y)
Z
gθ (y)dydσ(θ)
Sd−1 G
1
=
λd (G)σ(Sd−1 )
Z
Z
gθ (x + sθ)ds dxdσ(θ).
Prθ (G) L(x,θ)
|
{z
Ziel: ≥0
}
Wir betrachten
Z
Z
f (x + (r + s)θ)
drf (x + sθ)ds
λ1 (L(x + sθ, θ))
L(x,θ) L(x+sθ,θ)
1
=
λ1 (L(x, θ))
?
Z
Z
f (x + (r + s)θ)drf (x + sθ)ds
L(x,θ) L(x+sθ,θ)
1
=
λ1 (L(x, θ))
Z
Z
f (x + sθ)
L(x,θ)
L(x+e
r θ)
2

≥
1


λ1 (L(x, θ))
de
rds (mit r + s = re)
Z

f (x + sθ)ds
L(x,θ)
≥ 0.
Wir begründen ?: Wir wollen λ(L(x + sθ, θ)) = λ1 (L(x, θ)) zeigen. Dies gilt wegen
L(x + sθ, θ) = {r ∈ R : x + (r + s)θ ∈ G} + s
= {r + s ∈ R : x + (r + s)θ ∈ G}
= {t ∈ R : x + tθ ∈ G}
= L(x, θ).
Damit folgt die Behauptung hHf, f i ≥ 0.
Nun sind wir in der Lage die Fehlerabschätzung vom mittleren quadratischen Fehler aus Satz 4.11
anzuwenden.
Satz 4.39
Sei G ∈ K beliebig und (Xn )n∈N eine Markovkette mit Übergangskern H (auf (G, B(G))) und
Startverteilung πBd (Gleichverteilung in Bd ). Mit
n0 = d3 r2 log r · 2, 56 · 1014
84
4.3
erhalten wir
sup
(f,G)∈FK,4
eπBd (Sn,n0 , (f, G)) ≤
Anwendungen und Beispiele
3, 2 · 107 dr 5, 8 · 1015 d2 r2
√
+
.
n
n
Beweis:
1. Schritt: (Untere Schranke für Spektrallücke)
Wegen Lemma 3.39 (i) wissen wir, dass H positiv semidefinit ist. D.h. wir können die
Cheeger-Ungleichung (Satz 4.31) anwenden und erhalten mit β = kHkL02 →L02 , dass
1−β ≥
ϕ2
2
gilt.
Aus Satz 4.36 und Lemma 3.38 (ii) folgt
ϕ
≥
κ=10−3
1
1
min{1, (dr · 16000)−1 } =
2000
dr · 32 · 106
Somit gilt
1−β ≥
d−2 r−2
.
5, 12 · 10−14
2. Schritt: (Startverteilung)
Es gilt
Z
πBd (A) =
1Bd (x)
λd (G)πG (dx)
λd (Bd )
A
und wir haben
dπBd
1Bd (x)λd (G)
.
(x) =
dπG
λd (Bd )
Daraus folgt
2 Z
dπBd
=
−
1
dπG
2
2
1Bd (x)
dx
λ
(G)
−
1
λd (Bd ) d
λd (G)
G
=
λd (G)
−1
λd (Bd )
≤ rd − 1
≤r
(binom. Formel)
(G ⊆ rBd )
d
3. Schritt: (Fehlerabschätzung)
Bestimme n0 :
dπ
log dπBGd − 1
2
1−β
d
log rd2 r2 · 5, 12 · 1014
2
= 2, 56 · 1014 d3 r2 log r
≤
=: n0 .
85
4.3
Anwendungen und Beispiele
Dann gilt unter Verwendung vom 1. Schritt und Satz 4.11
√
√
2
128
sup eπBd (Sn,n0 (G, f )) ≤ √
+
.
1/2
n(1 − β)
n(1 − β)
(f,G)∈FK,4
Damit folgt die Behauptung.
Bemerkung 4.40
(a) Bemerkenswert ist, dass die Fehlerabschätzung polynomiell von d abhängt.
(b) Die vorkommenden Konstanten sind sehr groß. D.h. aus praktischer Sicht ist die Fehlerabschätzung wahrscheinlich nicht brauchbar. Für kleine d und r ist es unter Umständen besser die aus
Lemma 4.37 resultierende Fehlerabschätzung zu benutzen.
(c) Man könnte auch Satz 4.21 anwenden und damit für f ∈ L∞ Aussage über Konfidenzabschätzungen treffen.
4.3.2
Integration bezüglich nichtnormierter Dichtefunktion
Sei G ⊂ Rd beschränkt und fest. Sei ρ : G → (0, ∞) integrierbar bezüglich dem Lebesguemaß und
R
ρ(x)dx
, A ∈ B(G).
πρ (A) := RA
ρ(x)dx
G
Für c ≥ 1 definieren wir die Menge
Rc = {ρ : G → (0, ∞) | ∀x ∈ G : 1 ≤ ρ(x) ≤ c}.
Lemma 4.41
(a) Für alle ρ ∈ Rc gilt
supx∈G ρ(x)
sup ρ
=
≤ c.
inf ρ
inf x∈G ρ(x)
(b) Sei ρ : G → (0, ∞) und
sup ρ
inf ρ
≤ c. Dann folgt
c · ρ(x)
∈ Rc .
kρk∞
Beweis: Teil (a) ist klar und Teil (b) folgt aus
sup ρ
≤ ρ(x) ≤ sup ρ = kρk∞ .
c
Wir betrachten
FRc ,2 = {(f, ρ) : ρ ∈ Rc , f ∈ L2 (πρ ), kf k2 ≤ 1}.
Ziel: Berechne
R
S(f, ρ) =
86
f (x)ρ(x)dx
=
ρ(x)dx
G
GR
Z
f (x)πρ (dx)
G
4.3
Anwendungen und Beispiele
für ein beliebiges Paar (f, ρ) ∈ FRc ,2 .
Wir betrachten hier den/einen „unabhängigen Metropolis sampler“ (siehe Spezialfall in Abschnitt
3.2.1). Dabei bezeichnen wir mit µG die Gleichverteilung in G. Der Vorschlagskern ist gegeben durch
Q(x, A) = µG (A) :=
λd (A)
,
λd (G)
A ∈ B(G)
und der Metropoliskern mit Vorschlag µG und ρ ∈ Rc hat die Form


Z
Z
ρ(y)
dy
ρ(y)
dy
.
Mρ (x, A) = min 1,
+ 1A (x) 1 − min 1,
ρ(x) λd (G)
ρ(x) λd (G)
A
G
Mit Mρ bezeichnen wir den Übergangskern und den zugehörigen Markovoperator.
Gegebene / benötigte Informationen:
1) Wir haben ein Orakel welches uns Funktionswerte von f liefert und ein Orakel welches Funktionswerte von ρ liefert.
2) Wir haben einen Zufallszahlengenerator für die Gleichverteilung in G („G hat eine recht einfache
Struktur.“).
MCMC Methode
Sei (f, ρ) ∈ FRc ,2 und (Xn )n∈N eine Markovkette mit Übergangskern Mρ und Startverteilung µG .
Berechne
n
1X
f (Xj+n0 )
Sn,n0 (f, ρ) =
n j=1
als Approximation von S(f, ρ).
Lemma 4.42
Für alle ρ ∈ Rc , x ∈ G und A ∈ B(G) gilt
Mρ (x, A) ≥ c−1 πρ (A).
Das heißt unter Verwendung von Satz 4.23 folgt
kMρn − SkL1 →L1 ≤ 2(1 − c−1 )n .
Beweis: Wir haben
ρ(y)
dy
min 1,
ρ(x) λd (G)
A
R
Z
ρ(z)dz
−1
−1
= min{ρ(y) , ρ(x) } G
πρ (dy)
λd (G)
A


Z
Z
Z ρ(z)
ρ(z)  πρ (dy)
= min
dz,
dz
 ρ(y)
ρ(x)  λd (G)
Z
Mρ (x, A) ≥
A
−1
≥c
G
πρ (A)
(wegen
G
ρ(z)
ρ(y)
≥
inf ρ
sup ρ
≥ c−1 ).
87
4.3
Anwendungen und Beispiele
Wir können nun die Fehlerabschätzung aus Satz 4.9 für L1 -exponentiell konvergente Übergangskerne
anwenden.
Satz 4.43
Für ρ ∈ Rc sei (Xn )n∈N eine Markovkette mit Übergangskern Mρ und Startverteilung µG . Mit
n0 = c log(2c)
erhalten wir
√
2c 2c
e(Sn,n0 , (f, ρ)) ≤ √ + .
n
n
(f,ρ)∈FRc ,2
sup
Beweis:
1. Schritt: (Konvergenzeigenschaft)
Aus Lemma 4.42 und wegen der Reversibilität folgt, dass Mρ L1 -exponentiell konvergent
mit (1 − c−1 , 2) und
β := kMρ kL02 →L02 ≤ 1 − c−1 .
2. Schritt: (Startverteilung)
Es gilt
Z
µG (A) =
ρ(x)
ρ(y)dy
G
R
A
und somit
Wir erhalten
R
G
ρ(y)dy dx
=
ρ(x) λd (G)
Z R
ρ(y)dy
G
πρ (dx)
ρ(x)λd (G)
A
R
ρ(y)dy
dµG
.
(x) = G
dπρ
ρ(x)λd (G)
dµG
dπρ − 1
∞
≤
a,b≥0
|a−b|≤max{a,b}
dµG max 1, ≤ c.
dπρ ∞
3. Schritt: (Fehlerabschätzung bzw. Wahl von n0 )
G
log M dµ
−
1
dπρ
∞
1−α
≤
c · log(2c) =: n0
M =2
1−α=c−1
und wir erhalten mit Satz 4.9 die Behauptung.
Bemerkung 4.44
Prinzipiell ist für ρ ∈ Rc die Konstruktion von Zufallszahlengeneratoren möglich. Die Akzeptanzwahrscheinlichkeit einer einfachen Verwerfungsmethode zur Simulation von πρ liegt mindestens bei
c−1 .
88
4.3
Anwendungen und Beispiele
Die obige Fehlerabschätzung ist unabhängig von d, solange c unabhängig von der Dimension ist. Dies
ist unter Umständen nicht der Fall. Zum Beispiel mit α ≥ 0 und
ρ(x) = exp(−α|x|2 ).
√
Für G = Rd wäre dies die nichtnormierte Dichte von N (0, 2α−1 ).
Falls nun G = Bd , dann ist exp(α) · ρ ∈ Rexp(α) , da
supx∈Bd ρ(x)
1
≤ −α = eα
inf x∈Bd ρ(x)
e
und wegen Lemma 4.41.
Und falls G = [−1, 1]d , so ist exp(αd) · ρ ∈ Rexp(αd) , da
supx∈[−1,1]d ρ(x)
1
≤ −αd = eαd .
inf x∈[−1,1]d ρ(x)
e
Das heißt, abhängig von G, kann c exponentiell von d abhängen. Die Abhängigkeit von α ist in beiden
Fällen exponentiell.
Wir hätten gerne eine Fehlerabschätzung, die nur von log c abhängt. Für die Klasse FRc ,2 ist dies
leider nicht möglich.
Wir formulieren folgendes Ergebnis ohne Beweis.
Satz 4.45 (folgt aus [13, Theorem 1])
Für jeden randomisierten Algorithmus Sn (Zufallsvariable Sn ), der n Funktionswerte von f und ρ
benutzen darf, gilt
sup
e(Sn , (f, ρ)) =
(f,ρ)∈FRc ,2
E|Sn (f, ρ) − S(f, ρ)|2
sup
2
(f,ρ)∈FRc ,2
√
2
≥
·
6
(p
c
2n
3c
c+2n−1
2n ≥ c − 1
.
2n < c − 1
Die Klasse FRc ,2 ist einfach zu groß! Wir betrachten daher eine kleinere Klasse von Dichtefunktionen.
Wir beginnen dafür mit der Definition und Eigenschaften von log-konkaven Funktionen auf Rd .
Definition 4.46
Eine Funktion f : Rd → [0, ∞) heißt log-konkav, falls
f (x) = exp(−ϕ(x)),
x ∈ Rd
für eine konvexe Funktion ϕ : Rd → (−∞, ∞] gilt (e−∞ = 0).
Weiterhin sei für f : Rd → [0, ∞)
supp f = {x ∈ Rd : f (x) > 0}
der Support von f , auch Träger genannt, und für t ≥ 0 sei
Lt (f ) = {x ∈ Rd : f (x) ≥ t}
die Levelmenge von f zu Level t.
89
4.3
Anwendungen und Beispiele
Beispiel 4.47
1) Sei G ⊆ Rd konvex und λd (G) ∈ (0, ∞), dann ist die Dichte der Gleichverteilung in G
ρe(x) =
1G (x)
λd (G)
log-konkav.
Dies gilt, da
ρe(x) = exp(−ϕ(x) − λd (G))
(
0
für ϕ(x) =
∞
x∈G
x∈
/G
gilt.
2
2) Die Funktion ρe(x) = exp(−α|x|
), mit α > 0 ist log-konkav. Diese Funktion ist die nichtnormierte
√
−1
Dichtefunktion von N (0, 2α Id ). Allgemein gilt, dass die Dichte von N (µ, R) mit µ ∈ Rd und
Kovarianzmatrix R (symmetrisch, positiv definit) log-konkav ist.
Eigenschaften log-konkaver Funktionen
Sei ρe : Rd → [0, ∞).
E1: Folgende Aussagen sind äquivalent
(a) ρe ist log-konkav,
(b) ∀x, y ∈ Rd ∀λ ∈ [0, 1] : ρe(λx + (1 − λ)y) ≥ ρe(x)λ ρe(y)1−λ ,
(c) supp ρe ist konvex und log ρe ist konkav auf supp ρe.
Beweis:
(a) =⇒ (b) ρe log-konkav =⇒ ∃ϕ : Rd → (−∞, ∞] konvex mit ρe(x) = exp(−ϕ(x)). Wegen ϕ
konvex, folgt
∀x, y ∈ Rd ∀λ ∈ [0, 1] : −ϕ(λx + (1 − λ)y) ≥ −λϕ(x) − (1 − λ)ϕ(y).
Damit folgt
∀x, y ∈ Rd ∀λ ∈ [0, 1] : e−ϕ(λx+(1−λ)y) ≥ e−λϕ(x) e−(1−λ)ϕ(y) .
Die Definition von ϕ liefert schließlich die Behauptung
∀x, y ∈ Rd ∀λ ∈ [0, 1] : ρe(λx + (1 − λ)y) ≥ ρe(x)λ ρe(y)(1−λ) .
(b) =⇒ (c) supp ρe ist konvex, da ∀x, y ∈ supp ρe ∀λ ∈ [0, 1] gilt
(b)
ρe(λx + (1 − λ)y) ≥ ρe(x)λ ρe(y)1−λ
>
0.
x,y∈supp e
ρ
D.h. reicht es genügt nun zu zeigen, dass log ρe konkav auf supp ρe. Seien x, y ∈ supp ρe und
λ ∈ [0, 1] beliebig. Dann haben wir
(b)
ρe(λx + (1 − λ)y ) ≥ ρe(x)λ ρe(y)1−λ
|
{z
}
∈supp e
ρ
⇐⇒ log ρe(λx + (1 − λ)y) ≥ λ log ρe(x) + (1 − λ) log ρe(y).
Also ist log ρe konkav auf supp ρe.
90
4.3
Anwendungen und Beispiele
(c) =⇒ (a) Da log ρe konkav auf supp ρe, ist − log ρe konvex auf supp ρe und wir haben
ρe(x) = exp(−ϕ(x))
mit
(
− log ρe(x) x ∈ supp ρe
.
ϕ(x) =
∞
sonst
E2: Falls ρe auf supp ρe konkav und supp ρe konvex, dann ist ρe log-konkav.
Beweis: Wegen E1 g.z.z., dass log ρe konkav auf supp ρe. Wir haben für x, y ∈ supp ρe, λ ∈ [0, 1]
log{e
ρ(λx + (1 − λ)y)} ≥ log{λe
ρ(x) + (1 − λ)e
ρ(y)} ≥ λ log ρe(x) + (1 − λ) log ρe(y).
E3: Sind ρe1 , ρe2 log-konkav, dann auch ρe1 ρe2 und min{e
ρ1 , ρe2 } log-konkav.
E4: ρe log-konkav =⇒ ∀t ∈ (0, ∞) ist Lt (e
ρ) konvex.
Beweis: Seien x, y ∈ Lt (e
ρ), λ ∈ [0, 1]. Dann ist
E1 (b)
ρe(λx + (1 − λ)y)
≥
ρe(x)λ ρe(y)1−λ ≥ tλ t1−λ = t.
Bemerkung 4.48
Funktionen deren Levelmengen konvex sind, werden auch quasi-konvex gennant.
quasi-konvex
nicht quasi-konvex
In der Literatur wird manchmal der Begriff „unimodal“ für quasi-konvex verwendet. „Unimodal“ kann
aber auch heißen, dass eine Funktion nur einen „Mode = Modus“ hat, d.h. die Maximalstellen bilden
eine konvexe Menge.
Folgendes Resultat über log-konkave Funktionen ist sehr wichtig.
91
4.3
Anwendungen und Beispiele
Lemma 4.49 (Isoperimetrische Ungleichung)
Sei ρe log-konkav mit beschränktem „support“, d.h.
D := diam(supp ρe) =
sup
|x − y| < ∞.
x,y∈supp e
ρ
Es sei T1 , T2 , T3 eine paarweise disjunkte Zerlegung von supp ρe. Dann gilt
2
πe
(T ) ≥
inf |x − y| · min{πe
(T ), πe
(T )}.
ρ 3
ρ 1
ρ 2
D x∈T1
y∈T2
Beweis: Einen Beweis dieser Behauptung und weitere Ungleichungen dieser Art finden sich in [20]. Wir wollen nun eine Klasse von Dichtefunktionen mit zusätzlicher Struktur einführen.
Für α ≥ 0 sei nun ρ ∈ Lα , falls ρ : G → (0, ∞), G ⊆ Rd , G konvex und
a) Es existiert eine log-konkave Funktion ρe mit supp ρe = G und ρ = ρe auf G.
b) Für alle x, y ∈ G gilt | log ρ(x) − log ρ(y)| ≤ α|x − y|.
Wir können sagen, dass Lα log-konkave, log-Lipschitz-stetige, nichtnormierte Dichtefunktionen enthält.
Nun betrachten wir
FLα ,4 = {(f, ρ) : ρ ∈ Lα , f ∈ L4 (πρ ), kf k4 ≤ 1}.
Ziel: Für geeignete Mengen G berechne
R
Z
f (x)ρ(x)dx
S(f, ρ) = GR
= f (x)πρ (dx),
ρ(x)dx
G
(f, ρ) ∈ FLα ,4 .
G
Gegebene Informationen:
1) Funktionsauswertungen von f
2) Funktionsauswertungen von ρ
3) Zufallszahlengenerator zur Simulation der Gleichverteilung in G.
Als Markovkette schauen wir uns einen Metropolis-Algorithmus mit „ball-walk“ als Vorschlagskern
an. Für δ > 0 ist der Übergangskern des „δ-ball-walk“ gegeben durch
λd (G ∩ B(x, δ))
λd (A ∩ B(x, δ))
+ 1A (x) 1 −
Wδ (x, A) =
,
λd (B(0, δ))
λd (B(0, δ))
für A ∈ B(G) und x ∈ G. Und der Metropoliskern mit Vorschlag Wδ sowie ρ ∈ Lα hat die Form


Z
Z
ρ(y)
ρ(y)
Mρ,δ (x, A) = min{1, ρ(x) }Wδ (x, dy) + 1A (x) 1 − min{1, ρ(x) }Wδ (x, dy) .
A
G
fρ,δ (x, A) = 1 (Mρ,δ (x, A) + 1A (x)) die „lazy version“ von Mρ,δ . Wir betrachten M
fρ,δ um sicherSei M
2
zustellen, dass der Markovoperator positiv semidefinit ist, siehe Abschnitt 4.2.2.
92
4.3
Anwendungen und Beispiele
fρ,δ , Markovoperator und Markovkern.
Notation: Wir bezeichnen mit Mρ,δ , bzw. M
fρ,δ und
MCMC Methode: Sei (f, ρ) ∈ FLα ,4 und (Xn )n∈N eine Markovkette mit Übergangskern M
Startverteilung µG (Gleichverteilung in G). Berechne
n
Sn,n0 (f, ρ) =
1X
f (Xj+n0 )
n j=1
als Approximation von S(f, ρ).
Wir benötigen einige Eigenschaften des „δ-ball-walk“. Zunächst definieren wir für einen allgemeinen
Übergangskern K die local conductance
`K (x) = K(x, G \ {x}).
Interpretation: Sei r > 0 und Ar := {x ∈ G : `K (x) < r} und π(Ar ) ∈ (0, 1/2], sowie K reversibel
bezüglich π. Dann gilt
R
K(x, Ac )π(dx)
A
ϕ=
inf
π(A)
π(A)∈(0,1/2)
Z
π(dx)
≤ K(x, Acr )
π(Ar )
Ar
Z
π(dx)
≤ K(x, G \ {x})
|
{z
} π(Ar )
Ar
<r
< r.
D.h., intuitiv, falls r sehr klein und Ar genügend groß (π(Ar ) ∈ (0, 1/2]) haben wir kleines ϕ. Für Wδ
und Mρ,δ zeigen wir die Umkehrung.
Ziel: Falls `Wδ (x) ≥ ` > 0, bzw. `Mρ,δ (x) ≥ ` > 0 für alle x ∈ G und ` „groß genug“, dann ist auch
ϕ(Wδ ) bzw. ϕ(Mρ,δ ) „groß“.
Für Wδ erhalten wir
`Wδ (x) =
λd (G ∩ B(x, δ))
λd (B(0, δ))
Skizze:
x1
G
x2
`Wδ (x1 ) = 1/4
`Wδ (x2 ) = 1/2
`Wδ (x3 ) = 1
x3
93
4.3
Anwendungen und Beispiele
Im Folgenden benutzen wir öfter die Ungleichung
λd−1 (Bd−1 )
≤
λd (Bd )
r
d+1
.
2π
(BR)
Für einen Beweis siehe [13]. (BR) folgt im Wesentlichen aus
Γ(z + 1/2) √
≤ z,
Γ(z)
mit der Gammafunktion Γ(z) =
R∞
0
tz−1 e−t dt.
Lemma 4.50
1
, G = Bd und beliebiges x ∈ Bd gilt `Wδ (x) ≥ 0, 3.
Für δ ≤ √d+1
Beweis: Es ist klar, dass `Wδ (x) für x ∈ ∂Bd = Sd−1 am kleinsten ist. Es sei H die Tangentialhyperebene am Punkt x bezüglich Bd . Dann halbiert diese B(x, δ).
H
Bd
h
x
B(x, δ)
Grundfläche von Z(h)
Z(h)
Desweiteren sei Z(h) der d-dimensionale Zylinder, der die (d − 1)-dimensionale Kugel um x mit
Radius δ als Grundfläche besitzt. Die Höhe h des Zylinders ergibt sich aus dem Abstand von H und
der Hyperebene, die durch Bd ∩ ∂B(x, δ) bestimmt wird. Diese Höhe h kann man mit Hilfe des Satzes
von Thales und einer Verhältnisungleichung genau berechnen.
94
4.3
Anwendungen und Beispiele
C
δ
h
A
h
B
δ
2
∆ABC liegt auf Halbkreis über AB, also liegt nach Satz des Thales bei C ein rechter Winkel. Damit
2
folgt mit Strahlensatz hδ = 2δ , also h = δ2 .
Nach obiger Konstruktion können wir `Wδ (x) wie folgt abschätzen
`Wδ (x) =
Da λd (Z(h)) = h · λd−1 (δBd−1 ) =
λd (B(x, δ) ∩ Bd )
1
λd (Z(h))
≥ −
.
λd (B(0, δ))
2 λd (B(0, δ))
δ 2 d−1
λd−1 (Bd−1 )
2 δ
erhalten wir
1 hδ d−1 λd−1 (Bd−1 )
−
2
δ d λd (B2 )
1
λd−1 (Bd−1 )
=
1−δ
2
λd (Bd )
√
1
d−1
≥
(mit (BR)).
1−δ √
2
2π
`Wδ (x) ≥
√
Für δ ≤ 1/
d+1
bekommen wir
`Wδ (x) ≥
1
√
(1 − 1/ 2π) ≥ 0, 3
2
und damit ist der Beweis abgeschlossen.
Wir benötigen 2 weitere technische Lemmas über Wδ .
Lemma 4.51
Sei δ > 0 und t ∈ (0, 1). Falls x, y ∈ Rd mit
√
2π
|x − y| ≤ tδ √
,
d+1
dann ist
λd (B(x, δ) ∩ B(y, δ)) ≥ (1 − t)λd (B(0, δ)).
95
4.3
Anwendungen und Beispiele
Beweis: Sei u := |x − y|. Falls u < δ, dann ist das Lebesguemaß von B(x, δ) ∩ B(y, δ) gleich dem
Lebesguemaß von B(0, δ) ohne Mittelstück mit Breite u.
Skizze:
Mittelstück der Breite u
x
u
u
y
Grundfläche δBd−1 von Zylinder
Das Lebseguemaß dieses Mittelstückes kann nach oben durch das Lebesguemaß des Zylinders mit
Grundfläche δBd−1 und Höhe u abgeschätzt werden.
Das führt zu
λd (B(x, δ) ∩ B(y, δ) ≥ λd (δBd ) − uλd−1 (Bd−1 )
λd−1 (δBd−1 )
= λd (δBd ) 1 − u
λd (δBd )
u λd−1 (Bd−1 )
= λd (Bd ) 1 −
δ λd (Bd )
!
r
u d+1
(mit (BR))
≥ λd (δBd ) 1 −
δ
2π
p
≥ (1 − t)λd (δBd ) (wegen u ≤ tδ 2π/d+1).
Lemma 4.52
Für Wδ auf G sei `Wδ (x) ≥ ` für alle x ∈ G, d.h. ` ist eine untere Schranke für die „local conductance“. Dann gilt für alle t ∈ (0, `) und A ∈ B(G) mit
A1 := {x ∈ A : Wδ (x, Ac ) < 21 (` − t)},
A2 := {y ∈ Ac : Wδ (y, A) < 21 (` − t)},
dass
r
inf |x − y| > tδ
x∈A1
y∈A2
2π
.
d+1
Beweis: Der Beweis erfolgt mit Widerspruch. Angenommen
r
2π
.
inf |x − y| ≤ tδ
x∈A1
d+1
y∈A2
96
4.3
Anwendungen und Beispiele
Dann folgt
λd (B(x, δ) ∩ B(y, δ) ∩ G)
= λd ((B(x, δ) ∩ G) \ (B(x, δ) \ B(y, δ)))
(wegen A ∩ B ∩ C = (A ∩ B) \ (A \ C))
≥ λd (B(x, δ) ∩ G) − λd (B(x, δ) \ B(y, δ))
= λd (B(x, δ) ∩ G) − λd (B(x, δ) \ [B(x, δ) ∩ B(y, δ)])
= λd (B(x, δ) ∩ G) − λd (δBd ) − λd (B(x, δ) ∩ B(y, δ)).
Wegen obiger Annahme können wir Lemma 4.51 anwenden. Damit und mit
`Wδ =
λd (Bx,δ ∩ G)
≥`
λd (B(0, δ))
folgt
λd (B(x, δ) ∩ B(y, δ) ∩ G) ≥ `λd (δBd ) − λd (δBd ) + (1 − t)λd (δBd ).
(19)
Nun gilt
B(x, δ) ∩ B(y, δ) ∩ G ⊆ (B(x, δ) ∩ G ∩ Ac ) ∪ (B(y, δ) ∩ G ∩ A).
Wegen der Monotonie des Lebesguemaßes folgt
λd (B(x, δ) ∩ B(y, δ) ∩ G) ≤ λd ((B(x, δ) ∩ G ∩ Ac ) ∪ (B(y, δ) ∩ G ∩ A))
≤ λd (B(x, δ) ∩ G ∩ Ac ) + λd (B(y, δ) ∩ G ∩ A).
Aus (19) und vorheriger Ungleichung folgt
λd (B(x, δ) ∩ G ∩ Ac ) + λd (B(y, δ) ∩ G ∩ A) ≥ (` − t)λd (B(0, δ)).
Mit Wδ ausgedrückt ergibt dies
Wδ (x, Ac ) + Wδ (y, A) ≥ ` − t.
Dies ist ein Widerspruch zu x ∈ A1 und y ∈ A2 .
Folgerung 4.53
Für Wδ auf G ⊆ Rd sei `Wδ (x) ≥ ` für alle x ∈ G. Desweiteren sei G konvex und beschränkt mit
D := diam(G). Dann gilt für d ≥ 6 und δ < D
r
2π
`2 δ
.
ϕ(Wδ ) ≥
8 D d+1
Insbesondere ist für G = Bd
√
9 π
1
.
ϕ(Wd ) ≥ √
2 · 800 d + 1
Beweis: Mit Lemma 4.52 erhalten wir für t ∈ (0, `) und A ∈ B(G) mit
A1 := {x ∈ A : Wδ (x, Ac ) < 21 (` − t)},
A2 := {y ∈ Ac : Wδ (y, A) < 21 (` − t)},
dass
r
inf |x − y| > tδ
x∈A1
y∈A2
2π
.
d+1
97
4.3
Anwendungen und Beispiele
Da G konvex ist, folgt
1G (x)
λd (G)
ist log-konkav. Aus Lemma 4.49 folgt mit A3 = G \ (A1 ∪ A2 )
2
µG (A3 ) ≥ tδ
D
r
2π
min{µG (A1 ), µG (A2 )}.
d+1
Da Wδ reversibel bezüglich µG erhalten wir mit Satz 4.36
(
)
r
tδ
2π
`−t
min 1,
ϕ(Wδ ) ≥
2
D d+1
r
` − t tδ
2π
≥
(wegen d ≥ 6, δ ≤ D, t ∈ (0, `) und ` ≤ 1)
2 D d+1
r
`2 δ
2π
=
(wegen t(` − t) maximal für tmax = `/2).
8 D d+1
Der letzte Teil folgt mit Lemma 4.35.
Wir können dieses Ergebnis verallgemeinern.
Satz 4.54
Sei ρ ∈ Lα , `Wδ ≥ ` für alle x ∈ G und D := diam(G) < ∞. Dann gilt
r
π
`δ
`e−αδ
√
min 1,
.
ϕ(Mρ,δ ) ≥
8
2 D d+1
Beweis:
1. Schritt: (Mρ,δ (x, A) ≥ e−αδ Wδ (x, A))
Für ρ ∈ Lα gilt mit y ∈ B(x, δ), x ∈ G und der log-Lipschitz-Bedingung, dass
n
o
ρ(y)
min 1,
= min 1, elog ρ(y)−log ρ(x) ≥ e−α|x−y| ≥ e−αδ .
ρ(x)
Damit erhalten wir
Z
Mρ,δ (x, A) ≥
ρ(y)
min 1,
ρ(x)
Wδ (x, dy)
A
ρ(y)
Wδ (x, dy)
= 1B(x,δ) (y) min 1,
ρ(x)
A
Z
≥ e−αδ 1B(x,δ) (y)Wδ (x, dy)
Z
A
= e−αδ Wδ (x, A).
2. Schritt: (Verallgemeinerung von Lemma 4.52)
Für A ∈ B(G) sei
T1 :=
98
c
−αδ
x ∈ A : Mρ,δ (x, A ) < e
`
4
,
4.3
T2 :=
y ∈ Ac : Mρ,δ (y, A) < e−αδ
Anwendungen und Beispiele
`
4
.
Wir zeigen T1 ⊆ A1 und T2 ⊆ A2 mit A1 , A2 aus Lemma 4.52 und t = `/2.
Es gilt für x ∈ T1 , dass
`
4
eαδ Mρ,δ (x, Ac ) <
≤
Wδ (x, A)
x∈T1
1. Schritt
und somit x ∈ A1 . Analog zeigt man T2 ⊆ A2 . Damit erhalten wir
r
Lemma 4.52 `
2π
inf |x − y| ≥ inf |x − y|
>
δ
.
x∈T1
x∈A1
2
d+1
y∈T2
y∈A2
3. Schritt: (Anwendung isoperimetrischer Ungleichung, Lemma 4.49)
Es sei T3 := G \ (T1 ∪ T2 ) und da ρ log-konkav, folgt mit Lemma 4.49 und dem 2. Schritt,
dass
r
2 `
2π
δ
min{πρ (T1 ), πρ (T2 )}.
πρ (T3 ) ≥
D2
d+1
4. Schritt: (Anwendung von Satz 4.36)
Da Mρ,δ reversibel bezüglich πρ und mit den Schritten 1 bis 3 erhalten wir unter Verwendung von Satz 4.36
(
)
r
2π
` −αδ
` δ
ϕ(Mρ,δ ) ≥ e
min 1,
8
2D d+1
und die Behauptung des Satzes ist bewiesen.
Folgerung 4.55
Es sei ρ ∈ Lα mit G = Bd und δ ≤ (d + 1)−1/2 . Dann gilt
r
9δe−αδ
π
√
.
ϕ(Mρ,δ ) ≥
2 1600 · d + 1
Um ϕ(Mρ,δ ) zu maximieren, setzen wir
δ ∗ := min
und erhalten
r
ϕ(Mρ,δ ) ≥
√
1
1
,
d+1 α
π 9e−1
1
√
min
2 1600 d + 1
√
1
1
,
d+1 α
.
Beweis: Die Aussagen sind direkte Folgerungen aus Satz 4.54 und Lemma 4.50.
Nun haben wir alles Beisammen um eine Fehlerabschätzung für die MCMC Methode basierend auf
fρ,δ („lazy version“ von Mρ,δ ) zu formulieren und zu beweisen.
M
99
4.3
Anwendungen und Beispiele
Satz 4.56
Für G = Bd , ρ ∈ Lα und δ ∗ = min{(d + 1)−1/2 , α−1 } sei (Xn )n∈N eine Markovkette mit Übergangsfρ,δ∗ und Startverteilung µB . Mit
kern M
d
n0 := 1189362 α(d + 1) max{d + 1, α2 }
erhalten wir
√
1543 √
13456095
eµBd (Sn,n0 , (f, ρ)) ≤ √
d + 1 max{ d + 1, α} +
(d + 1) max{d + 1, α2 }.
n
n
(f,ρ)∈FLα ,4
sup
Beweis:
1. Schritt (Konvergenzeigenschaft)
Mit Lemma 4.35, Folgerung 4.40 und der Wahl von δ ∗ erhalten wir
fρ,δ∗ ) = ϕ(Mρ,δ∗ ) ≥
ϕ(M
2
r
π 9
1
√
· min
2 3200 e d + 1
1
1
√
,
α
d+1
Mit der Cheeger-Ungleichung (Satz 4.31) erhalten wir dann
81
1
fρ,δ kL0 →L0 ≥ 1 π
min
1 − kM
2
2
2
2
|2 2 3200
{z e } d + 1
1
1
, 2
d+1 α
≈0,841·10−6
2. Schritt (Startverteilung)
Es gilt
Z
µBd (A) =
R
A
ρ(A)
·
ρ(y)dy
Bd
Z Z
ρ(y)dy
=
R
Bd
ρ(y)dy
ρ(x)
dx
λd (Bd )
1
πρ (dx)
ρ(x)λd (Bd )
A Bd
und somit
dµBd
1
(x) =
·
dπρ
ρ(x)
Z
ρ(y)
Bd
100
dy
.
λd (Bd )
.
.
4.3
Mit c :=
R
Bd
Anwendungen und Beispiele
ρ(y)dy folgt
2 Z
dµBd
=
−
1
dπρ
2
Bd
2
c
ρ(x)
−1
ρ(x)
−
1
λd (Bd )
c dx
| {z }
Z =
πρ (dx)
2
2
1
c
c
−2
+1
λd (Bd )2 ρ(x)
λd (Bd )ρ(x)
ρ(x)
dx
c
Bd
Z
c
dx − 1
λd (Bd )2 ρ(x)
Bd
Z Z
dy
ρ(y) dx
≤
ρ(x) λd (Bd ) λd (Bd )
=
(c einsetzen)
Bd Bd
(?)
≤ e2α ,
wobei (?) wegen
ρ(y)
= exp(log(ρ(y)) − log(ρ(x))) ≤ eα|x−y|
ρ(x)
≤
e2α
x,y∈Bd
gilt.
Das heißt
dµBd
α
dπρ − 1 ≤ e .
2
3. Schritt (Fehlerabschätzung, Anwendung von Satz 4.11)
Bestimme n0 : Es gilt
dν
log(k dπ
− 1k)
ρ
1−β
1. und 2. Schritt
≤
α(d + 1) max{d + 1, α2 } · 1189362.
Dann gilt unter Verwendung vom 1. Schritt und Satz 4.11 die Fehlerabschätzung.
101
LITERATUR
Literatur
[1] G. Alsmeyer. Wahrscheinlichkeitstheorie (Vorlesungsskript). http://www.math.uni-muenster.
de/statistik/lehre/SS13/WT/Skript/WT-Buch_8.pdf, 2013. [Online; Zugriff am 31.07.2015].
[2] S. Asmussen and P. Glynn. A new proof of convergence of MCMC via the ergodic theorem.
Statist. Probab. Lett., 81:1482–1485, 2011.
[3] A. Barker. Monte Carlo calculations of the radial distribution functions for a proton-electron
plasma. Aust. J. Phys., 18:119–133, 1965.
[4] R. Dudley. Real analysis and probability, volume 74 of Cambridge Studies in Advanced Mathematics. Cambridge University Press, Cambridge, 2002. Revised reprint of the 1989 original.
[5] O. Forster. Analysis 3: Maß- und Integrationstheorie, Integralsätze im Rn und Anwendungen.
Wiesbaden: Springer Spektrum, 7th revised ed. edition, 2012.
[6] P. Glynn and D. Ormoneit. Hoeffding’s inequality for uniformly ergodic Markov chains. Statist.
Probab. Lett., 56(2):143–146, 2002.
[7] W. Hastings. Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57(1):97–109, 1970.
[8] O. Kallenberg. Foundations of Modern Probability. Probability and its Applications. SpringerVerlag, New York, Second edition, 2002.
[9] A. Klenke. Wahrscheinlichkeitstheorie. Springer, 2006. in German.
[10] K. Łatuszyński and D. Rudolf. Convergence of hybrid slice sampling via spectral gap.
http://arxiv.org/abs/1409.2709, 2014.
[11] G. Lawler and A. Sokal. Bounds on the L2 spectrum for Markov chains and Markov processes:
a generalization of Cheeger’s inequality. Trans. Amer. Math. Soc., 309(2):557–580, 1988.
[12] L. Lovász and S. Vempala. Hit-and-run from a corner. SIAM J. Comput., 35(4):985–1005, 2006.
[13] P. Mathé and E. Novak. Simple Monte Carlo and the Metropolis algorithm. J. Complexity,
23(4-6):673–696, 2007.
[14] N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, and E. Teller. Equation of state calculations by fast computing machines. Journal of Chemical Physics, 21(6):1087–1092, 1953.
[15] S. Meyn and R. Tweedie. Markov chains and stochastic stability. Cambridge University Press,
second edition, 2009.
[16] B. Miasojedow. Hoeffding’s inequalities for geometrically ergodic markov chains on general state
space. Preprint, Available at http://arxiv.org/abs/1201.2265, 2012.
[17] T. Müller-Gronbach, E. Novak, and K. Ritter. Monte Carlo-Algorithmen. Springer, 2012. in
German.
[18] D. Rudolf. Explicit error bounds for Markov chain Monte Carlo. Dissertationes Math., 485:93
pp., 2012.
[19] D. Rudolf and M. Ullrich. Positivity of hit-and-run and related algorithms. Electron. Commun.
Probab., 18:1–8, 2013.
102
LITERATUR
[20] S. Vempala. Geometric random walks: A survey. Combinatorial and computational geometry,
52:573–612, 2005.
[21] D. Werner. Funktionalanalysis. Springer Verlag, 2005.
103
INDEX
Index
L1 -exponentielle Konvergenz, 53
δ ball walk, 21
π-irreduzibel, 55
log-konkav, 89
Übergangsdichte, 34
Übergangskern, 18
Gibbs sampler
Deterministic scan Gibbs sampler, 41
Random scan Gibbs sampler, 40
bedingte Dichte, 16
bedingter Erwartungswert, 10
burn-in, 55
Doeblin, 69
Fehler, 5
Gibbs sampler, 40
gleichmäßig ergodisch, 54
Markovprozess, 18
quasi-konvex, 91
reguläre Version, 13
reguläre Version der bedingten Verteilung, 13,
14
Reversibilität, 27
Slice sampling, 42
Hybride slice sampler, 44
Simple slice sampler, 42
Spektrallücke, 53
Startverteilung, 18
stationäre Verteilung, 26
Support, 89
totaler Variationsabstand, 49
Träger, 89
homogene Markovkette, 18
Update-Funktion, 18
Leitfähigkeit, 71
Levelmenge, 42
local conductance, 93
Version der bedingten Wahrscheinlichkeit, 13
Version des bedingten Erwartungswertes, 9
Verteilungsdichte, 15
Vorschlagskern, 34
Markoveigenschaft, 18
Markovkern, 18
Markovkette, 7, 18
Markovoperator, 46
104
zeitlich diskret, 18
zeitlich stetig, 18
zufällige Vektoren, 15

Download Report