43
Methoden der Statistik
Mathias Trabs
1. September 2015
Inhaltsverzeichnis
1 Grundbegrie der Statistik 21.1 Drei grundlegende Fragestellungen . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.1.1 Schätzprobleme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.1.2 Hypothesentests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61.1.3 Kondenzmengen (Bereichsschätzung) . . . . . . . . . . . . . . . . . . . . . 9
1.2 Minimax- und Bayesansatz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101.3 Ergänzungen: Quantile . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2 Lineares Modell 152.1 Regression und kleinste Quadrate . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.2 Inferenz unter Normalverteilungsannahme . . . . . . . . . . . . . . . . . . . . . . . 212.3 Varianzanalyse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3 Exponentialfamilien and verallgemeinerte lineare Modelle 303.1 Die Informationsungleichung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303.2 Verallgemeinerte Lineare Modelle . . . . . . . . . . . . . . . . . . . . . . . . . . . . 343.3 Ergänzung: Numerische Bestimmung des Maximum-Likelihood-Schätzers . . . . . . 37
4 Klassikation 384.1 Logistische Regression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 384.2 Bayesklassikation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 404.3 Lineare Diskriminanzanalyse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
5 Ausblick 43
Literatur
Georgii, H.-O.: Stochastik. Einführung in die Wahrscheinlichkeitstheorie und Statistik, deGruyter, 2007
James, G., Witten, D., Hastie, T., Tibshirani, R.: An Introduction to Statistical Learning(with Applications in R), Springer, 2013
Lehmann, E.L., Romano, J.P.: Testing Statistical Hypotheses, Springer, 2005
Lehmann, E.L. and G. Casella: Theory of Point Estimation, Springer, 2003
Shao, J: Mathematical Statistics, Springer, 2003
Venables, W.N., Ripley, B.D.: Modern Applied Statistics with S-Plus, Springer, 1997
Wasserman, L.: All of Statistics, Springer, 2003
Witting, H.: Mathematische Statistik I, Teubner, 1985
1
1 Grundbegrie der Statistik
Während die Wahrscheinlichkeitstheorie anhand eines gegebenen Modells die Eigenschaften der(zufälligen) Ereignisse untersucht, ist das Ziel der Statistik genau andersherum: Wie kann manaus den gegebenen Beobachtungen Rückschlüsse auf das Modell ziehen?
Beispiel 1.1 (Werbung). Wir verwenden den Advertising-Datensatz aus James et al. (2013).Für 200 Märkte haben wir die Anzahl der verkauften Produkte Y sowie das jeweilige Budget fürFernsehwerbung XF , für Radiowerbung XR und für Zeitungsannoncen XZ gegeben.
Betrachten wir das Modell
Yi = aXFi + b+ εi, i = 1, . . . , 200,
wobei die zufälligen Störgröÿen εi Marktunsicherheiten, externe Einüsse etc. modellieren. Plau-sible Annahmen an das Modell sind
(i) (εi) sind unabhängig (näherungsweise),
(ii) (εi) sind identisch verteilt,
(iii) E[εi] = 0 (kein systematischer Fehler)
(iv) εi normalverteilt (wegen ZGWS).
Naheliegende Ziele/Fragestellungen:
(i) Es sollen a, b anhand der Daten ermittelt werden. Ein mögliches Schätzverfahren ist derKleinste-Quadrate-Schätzer
(a, b) := arg mina,b
n∑i=1
(Yi − aXi − b)2
(wir minimieren die Summe der quadrierten Residuen). Mit a, b erhalten wir die Regressi-onsgrade
y = axF + b.
(ii) Sind die Modellannahmen erfüllt? Histogramm, Boxplot und QQ-Plot (Quantil-Quantil-Plot) der Residuen.
(iii) Wenn wir die Verteilung von a kennen (Verteilungsannahme an ε nötig!), können wir Inter-valle der Form I = [a− c, a+ c] für c > 0 konstruieren, so dass der tatsächlich Parameter amit vorgegebener Wahrscheinlichkeit in I liegt.
(iv) Wir wollen testen, ob es einen Eekt gibt, d.h. gilt die Hypothese H0 : a = 0 oder kannsie verworfen werden? Beispielsweise kann man die Hypothese verwerfen, falls |a| > c füreinen kritischen Wert c > 0. Um einen sinnvollen Wert zu bestimmen, benötigen wir wiederVerteilungsannahmen an die Fehler (εi).
Wir können das Modell auf polynomielle Regression Yi = a0 + a1XFi + · · · + an(XF
i )n + εi odermultiple Regression Yi = a0 + a1X
Fi + a2X
Ri + a3X
Zi + εi erweitern. Dies führt auf das Problem
der Modellwahl.
Denition 1.2. Ein messbarer Raum (X ,F ) versehen mit einer Familie (Pϑ)ϑ∈Θ von Wahr-scheinlichkeitsmaÿen mit einer beliebigen Parametermenge Θ 6= ∅ heiÿt statistisches Experimentoder statistisches Modell. X heiÿt Stichprobenraum. Jede (F ,S )-messbare Funktion Y : X → Sheiÿt Beobachtung oder Statistik mit Werten in (S,S ) und induziert das statistische Modell(S,S , (PYϑ )ϑ∈Θ). Sind die Beobachtungen Y1, . . . , Yn für jedes Pϑ unabhängig und identisch ver-teilt (iid.), so nennt man Y1, . . . , Yn eine mathematische Stichprobe.
2
Beispiel 1.3 (mathematische Stichprobe). Für n ∈ N seien X1, . . . , Xn iid. verteilte Zufalls-variablen mit Werten in X und Randverteilung X1 ∼ Pϑ mit Parameter ϑ ∈ Θ. Dann ist derStichprobenvektor (X1, . . . , Xn) gemäÿ dem Produktmaÿ Pnϑ(dx) =
∏ni=1 Pϑ(dxi) auf (Xn,F⊗n)
verteilt.
Wir werden uns in dieser Vorlesung weitgehend mit (verallgemeinerten) linearen Modellenbefassen, d.h. die Abhängigkeit der Zufallsvariablen Xi bzw. deren Verteilung vom unbekanntenParameter kann durch eine lineare Abbildung dargestellt werden.
1.1 Drei grundlegende Fragestellungen
Die meisten statistischen Fragestellungen kann man einer der drei Grundprobleme Schätzen, Testenund Kondenzintervalle zuordnen. Diese werden im folgenden kurz umrissen und im Laufe derVorlesung weiter vertieft.
1.1.1 Schätzprobleme
Ziel ist es, aufgrund der vorhandenen Beobachtungen den unbekannten Parameter im statisti-schen Modell (X ,F , (Pϑ)ϑ∈Θ) zu bestimmen, also einen einzelnen (bestmöglichen) Wert anzuge-ben (Punktschätzung). Damit ist ein Schätzer eine Abbildung, die nur von den Beobachtungenabhängt.
Denition 1.4. Sei (X ,F , (Pϑ)ϑ∈Θ) ein statistisches Modell, ρ : Θ → Rd ein (abgeleiteter) d-dimensionaler Parameter, d ∈ N. Ein Schätzer ist eine messbare Abbildung ρ : X → Rd. GiltEϑ[ρ] = ρ(ϑ) so heiÿt ρ unverzerrt oder erwartungstreu (engl.: unbiased).
Beispiel 1.5. Seien X1, . . . , Xn eine Bernoulli-verteilte mathematische Stichprobe mit Parameterp ∈ (0, 1). Betrachte den Schätzer pn := n−1
∑ni=1Xi. Dann gilt Ep[pn] = n−1
∑ni=1 E[Xi] = p.
Also ist pn erwartungstreu. Um die Streuung des Schätzers um den wahren Parameter p zu messen,berechnen wir
Varp(pn) = n−2n∑i=1
Varp(Xi) =p(1− p)
n.
Für gröÿer werdenden Stichprobenumfang konzentriert sich also pn um p.
Wie gut ein Schätzer ist, wird mithilfe einer Verlustfunktion bestimmt. Diese misst den Abstandzwischen geschätztem und wahrem Parameter.
Denition 1.6. Eine Funktion L : Θ × Rd → R+ heiÿt Verlustfunktion, falls L(ϑ, ·) für jedesϑ ∈ Θ messbar ist. Der erwartete Verlust R(ϑ, ρ) := Eϑ[L(ϑ, ρ)] eines Schätzers ρ heiÿt Risiko.Typische Verlustfunktionen sind
(i) der 0-1-Verlust L(ϑ, r) = 1r 6=ρ(ϑ),
(ii) der absolute Verlust L(ϑ, r) = |r − ρ(ϑ)| (euklidischer Abstand im Rp) sowie
(iii) der quadratische Verlust L(ϑ, r) = |r − ρ(ϑ)|2.
Lemma 1.7 (Bias-Varianz-Zerlegung). Sei (X ,F , (Pϑ)ϑ∈Θ) ein statistisches Modell und ρ : X →Rd ein Schätzer des Parameters ρ(ϑ) mit Eϑ[|ρ|2] <∞ für alle ϑ ∈ Θ. Dann gilt für den quadra-tischen Verlust
Eϑ[|ρ− ρ(ϑ)|2
]= Varϑ(ρ) + |Eϑ[ρ]− ρ(ϑ)︸ ︷︷ ︸
Bias
|2 für alle ϑ ∈ Θ.
Beweis. Es gilt
Eϑ[|ρ− ρ(ϑ)|2
]=Eϑ
[|ρ− Eϑ[ρ] + Eϑ[ρ]− ρ(ϑ)|2
]=Eϑ[|ρ− Eϑ[ρ]|2] + 2Eϑ
[(ρ− Eϑ[ρ])>(Eϑ[ρ]− ρ(ϑ))
]+ |Eϑ[ρ]− ρ(ϑ)|2
= Varϑ(ρ) + |Eϑ[ρ]− ρ(ϑ)|2.
3
Beispiel. In der Situation von Beispiel 1.5, betrachten wir den Schätzer pn := (∑ni=1Xi+1)/(n+
2). Dieser hat den Bias
E[pn]− p =1− 2p
n+ 2
und die Varianz
Var(pn) =np(1− p)(n+ 2)2
.
Damit hat pn einen kleineren quadratischen Fehler als pn, wenn |p− 1/2| 6 1/√
8.
Bemerkung 1.8. Ein Schätzproblem, bei dem der interessierende Parameter nur endliche vieleWerte annehmen kann, heiÿt auch Klassikationsproblem und der entsprechende Schätzer heiÿtKlassizierer (mehr dazu in Kapitel 4).
Obwohl wir in dieser Vorlesung keine Asymptotik, d.h. das Verhalten der Schätzer bei Stich-probenumfängen n→∞, behandeln, seien noch zwei weitere wichtige Grundbegrie erwähnt.
Denition 1.9. Sei X1, . . . , Xniid.∼ Pϑ eine mathematische Stichprobe. Dann heiÿt ein Schätzer
ρn vom abgeleiteten Parameter ρ(ϑ) konsistent, falls
ρnPϑ→ ρ(ϑ) für n→∞.
Der Schätzer ρn heiÿt asymptotisch normalverteilt, falls E[|ρn|2] <∞ und
ρn − Eϑ[ρn]√Varϑ(ρn)
L→ N (0, 1) unter Pϑ.
Aufgrund des zentralen Grenzwertsatzes sind viele Schätzer asymptotisch normalverteilt, soauch in Beispiel 1.5. Daher kommt der Untersuchung von statistischen Modellen unter Normal-verteilungsannahme eine besondere Bedeutung zu.
Zwei wichtige Konstruktionsprinzipien von Schätzern sind die Momentenmethode undMaximum-Likelihood-Schätzer:
Methode 1: Momentenmethode. Sei X1, . . . , Xn eine mathematische Stichprobe reellerZufallsvariablen mit E[|X1|d] < ∞. Oensichtlich hängen i.A. die Momente einer Verteilungmk = mk(ϑ) := Eϑ[Xk
1 ], k ∈ N, von ihrem Parameter ϑ ∈ Rd ab. Aufgrund des Gesetzesder groÿen Zahlen ist der kanonische Schätzer von mk gegeben durch das Stichprobenmomentmk := 1
n
∑nj=1X
kj . Der Momentenschätzer ϑ von ϑ ist deniert als die Lösung der d-Gleichungen
m1(ϑ) = m1,
m2(ϑ) = m2,
......
...
md(ϑ) = md.
Beispiel 1.10. Sei X1, . . . , Xniid.∼ N (µ, σ2). Dann ist m1 = Eµ,σ2 [X1] = µ und m2 = Eµ,σ2 [X2
1 ] =Varµ,σ2(X1) + Eµ,σ2 [X1]2 = σ2 + µ2. Folglich müssen wir die Gleichungen
µ =1
n
n∑j=1
Xj und σ2 + µ2 =1
n
n∑j=1
X2j
lösen. Bezeichnen wir das Stichprobenmittel mit Xn := 1n
∑nj=1Xj , erhalten wir die Lösung
µ = Xn, σ2 =1
n
n∑j=1
(Xj −Xn)2.
4
Die Momentenmethode kann auf die Erwartungswerte allgemeinerer Funktionale verallgemei-nert werden (siehe Übung 2). Für die zweite Methode benötigen wir etwas mehr Struktur, die wirauch im weiteren Verlauf der Vorlesung immer wieder aufgreifen.
Denition 1.11. Ein statistisches Modell (X ,F , (Pϑ)ϑ∈Θ) heiÿt dominiert, falls es ein σ-endlichesMaÿ µ gibt, so dass Pϑ absolut stetig bzgl. µ ist (Pϑ µ) für alle ϑ ∈ Θ. Die durch ϑ parametri-sierte Radon-Nikodym-Dichte
L(ϑ, x) :=dPϑdµ
(x), ϑ ∈ Θ, x ∈ X
heiÿt Likelihoodfunktion, wobei diese meist als durch x parametrisierte Funktion in ϑ aufgefasstwird.
Beispiel 1.12.
(i) X = R,F = B(R),Pϑ ist gegeben durch die Lebesguedichte fϑ, beispielsweise Pµ,σ2 =N (µ, σ2) oder Pϑ = U([0, ϑ]). Dann ist L(ϑ, x) = fϑ(x).
(ii) Jedes statistische Modell auf dem Stichprobenraum (N,P(N)) oder allgemeiner auf einemabzählbaren Raum (X ,P(X )) ist vom Zählmaÿ dominiert. Die Likelihoodfunktion ist durchdie Zähldichte gegeben.
(iii) Ist Θ = ϑ1, ϑ2, . . . abzählbar, so ist µ =∑i ciPϑi mit ci > 0 und
∑i ci = 1 ein dominie-
rendes Maÿ.
Methode 2: Maximum-Likelihood-Prinzip. Für ein dominiertes statistisches Modell mitLikelihoodfunktion L(ϑ, x) heiÿt eine Statistik ϑ : X → Θ (Θ trage eine σ-Algebra) Maximum-Likelihood-Schätzer (MLE: maximum likelihood estimator), falls
L(ϑ, x) = supϑ∈Θ
L(ϑ, x) für Pϑ-f.a. x ∈ X und alle ϑ ∈ Θ.
Beispiel 1.13. Betrachten wir wieder eine mathematische StichprobeX1, . . . , Xn normalverteilterZufallsvariablen. Dann ist (Rn,B(Rn),Pnµ,σ2) mit Pµ,σ2 = N (µ, σ2) ein vom Lebesguemaÿ auf Rndominiertes Modell mit Likelihoodfunktion, x = (x1, . . . , xn) ∈ Rn,
L(µ, σ2;x) = (2πσ2)−n/2n∏j=1
exp(− (xj − µ)2
2σ2
).
Um den Maximum-Likelihood-Schätzer zu berechnen, nutzen wir die Monotonie des Logarithmusund betrachten
logL(µ, σ2;x) = −n2
(log(2π) + log σ2)−n∑j=1
(xj − µ)2
2σ2→ max
µ,σ2.
Ableiten nach µ und σ2 führt auf die Gleichungen
0 = σ−2n∑j=1
(xj − µ),n
2σ2=
1
2σ4
n∑j=1
(xj − µ)2.
Umstellen der ersten Gleichung nach µ liefert µ = Xn und Einsetzen in die zweite Gleichung ergibtσ2 = n−1
∑j(Xj−Xn)2. Es ist leicht nachzuprüfen, dass µ und σ2 tatsächlich das Maximierungs-
problem lösen (und messbar sind). In diesem Fall stimmt der Maximum-Likelihood-Schätzer alsomit dem Momentenschätzer überein.
5
Beispiel 1.14. Sei X1, . . . , Xn eine Poisson-verteilte mathematische Stichprobe mit Parameterλ > 0, d.h. X = Zn+,F = P(X ) und Pλ(X1 = k) = λke−k
k! . Dann ist die gemeinsame Verteilunggegeben durch
Pλ(X1 = k1, . . . , Xn = kn) =λ∑i kie−nλ
(k!)n, k1, . . . , kn ∈ Z+.
Ableiten nach λ und null setzen führt auf den Maximum-Likelihood-Schätzer λ = Xn (hinreichendeBedingung prüfen!).
1.1.2 Hypothesentests
Häug interessiert man sich weniger für die gesamte zugrunde liegende Verteilung, als die Frage,ob eine bestimmte Eigenschaft erfüllt ist, oder nicht. Beispielsweise möchte man wissen, ob eineneue Behandlungsmethode I besser ist als die alte bisher genutzte Methode II. Aufgrund einerBeobachtung soll entschieden werden, ob die Hypothese I ist besser als II akzeptiert werdenkann oder verworfen werden sollte.
Um derartige Fragestellungen in einem statistischen Modell (X ,F , (Pϑ)ϑ∈Θ) zu formalisieren,wird die Parametermenge in zwei disjunkte Teilmengen Θ0 und Θ1 zerlegt, d.h. Θ = Θ0 ∪Θ1 und∅ = Θ0 ∩Θ1. Das Testproblem liest sich dann als
H0 : ϑ ∈ Θ0 versus H1 : ϑ ∈ Θ1.
Dabei werden H0, H1 als Hypothesen bezeichnet, genauer heiÿt H0 Nullhypothese und H1 Alterna-tivhypothese oder Alternative. Ein statistischer Test entscheidet nun zwischen H0 und H1 aufgrundeiner Beobachtung x ∈ X .
Denition 1.15. Ein (nicht-randomisierter) statistischer Test ist eine messbare Abbildungϕ : (X ,F ) → (0, 1,P(0, 1)), wobei ϕ(x) = 1 heiÿt, dass die Nullhypothese verworfen/ dieAlternative angenommen wird und ϕ(x) = 0 bedeutet, dass die Nullhypothese nicht verworfenwird/ akzeptiert wird. Die Menge ϕ = 1 = x ∈ X : ϕ(x) = 1 heiÿt Ablehnbereich von ϕ.
Allgemeiner ist ein randomisierter statistischer Test eine messbare Abbildung ϕ : (X ,F ) →([0, 1],B([0, 1])). Im Fall ϕ(x) ∈ (0, 1) entscheidet ein unabhängiges Bernoulli-Zufallsexperimentmit Erfolgswahrscheinlichkeit p = ϕ(x), ob die Hypothese verworfen wird.
Testen beinhaltet mögliche Fehlerentscheidungen:
(i) Fehler 1. Art (α-Fehler, type I error): Entscheidung für H1, obwohl H0 wahr ist,
(ii) Fehler 2. Art (β-Fehler, type II error): Entscheidung für H0, obwohl H1 wahr ist.
Denition 1.16. Sei ϕ ein Test der Hypothese H0 : ϑ ∈ Θ0 gegen die Alternative H1 : ϑ ∈ Θ1
im statistischen Modell (X ,F , (Pϑ)ϑ∈Θ). Die Gütefunktion von ϕ ist deniert als
βϕ : Θ→ R+, ϑ 7→ Eϑ[ϕ]
Ein Test ϕ erfüllt das Signikanzniveau α ∈ [0, 1] (oder ϕ ist Test zum Niveau α), falls βϕ(ϑ) 6 αfür alle ϑ ∈ Θ0. Ein Test ϕ zum Niveau α heiÿt unverfälscht, falls βϕ(ϑ) > α für alle ϑ ∈ Θ1.
Somit hat ein nicht-randomisierten Test das Niveau α ∈ (0, 1), falls
Pϑ(ϕ = 1) 6 α, für alle ϑ ∈ Θ0,
beschränkt also die Wahrscheinlichkeit des Fehlers 1. Art mit der vorgegeben oberen Schranke α.In der Regel ist es nicht möglich, die Wahrscheinlichkeiten für die Fehler 1. und 2. Art gleichzeitigzu minimieren. Daher werden diese typischerweise asymmetrisch betrachtet:
(i) Begrenzung der Fehlerwahrscheinlichkeit 1. Art durch ein vorgegebenes Signikanzniveau α.
6
(ii) Unter der Maÿgabe (i) wird die Wahrscheinlichkeit für Fehler 2. Art minimiert.
Eine zum Niveau α statistisch abgesicherte Entscheidung kann also immer nur zu Gunsten vonH1 getroen werden. Daraus folgt die Merkregel Was nachzuweisen ist, stets als Alternative H1
formulieren.
Beispiel 1.17 (Einseitiger Binomialtest). Von den 13 Todesfällen unter 55- bis 65-jährigen Ar-beitern eines Kernkraftwerkes im Jahr 1995 waren 5 auf einen Tumor zurückzuführen. Die To-desursachenstatistik 1995 weist aus, dass Tumore bei etwa 1/5 aller Todesfälle die Ursache in derbetreenden Altersklasse (in der Gesamtbevölkerung) darstellen. Ist die beobachtete Häufung vontumorbedingten Todesfällen signikant zum Niveau 5%?
Bezeichne X die Anzahl der Tumortoten unter n = 13 Todesfällen. Dann ist das statistischeModell gegeben durch X = 0, . . . n,F = P(X ) und Pp = Bin(13, p) mit Parameter p ∈ [0, 1]und das Testproblem ist gegeben durch
H0 : p 6 1/5 versus H1 : p > 1/5.
Ziel ist ein nicht-randomisierter Test zum Niveau α = 0, 05. Naheliegenderweise konstruieren wirϕ(x) = 1x>c wobei der kritische Wert c > 0 so gewählt wird, dass supp61/5 Pp(X > c) 6 α.Um eine möglichst groÿe Güte zu erreichen, sollte c unter dieser Nebenbedingung möglichst kleingewählt werden. Für k ∈ X gilt
Pp(X 6 k) =
k∑l=0
(13
l
)pl(1− p)13−l.
Da p 7→ Pp(X 6 k) für alle k ∈ X monoton fallend auf [0, 1] ist (ableiten), folgt supp61/5 Pp(X >c) = P1/5(X > c). Wegen
P1/5(X 6 4) ≈ 0, 901 und P1/5(X 6 5) ≈ 0, 970,
wählen wir c = 5. Somit kann die Hypothese zum Niveau 0, 05 nicht verworfen werden. DieGütefunktion von ϕ
βϕ(p) = Pp(X > 5) =
13∑l=6
(13
l
)pl(1− p)13−l, p ∈ [0, 1],
ist monoton wachsend und somit ist ϕ auch unverfälscht.
Dieses Beispiel führt uns auf ein allgemeines Konstruktionsprinzip von Tests einer HypotheseH0 : ϑ ∈ Θ0 vs. H1 : ϑ ∈ Θ1 mit Θ0 6= ∅ und Θ1 = Θ \Θ0.
Methode 3: Teststatistiken. Für Ablehnbereiche (Γα)α∈(0,1) ⊆ B(R) und eine TeststatistikT : (X ,F )→ (R,B(R)) sei ein Test gegeben durch
ϕ(x) = 1T (x)∈Γα, x ∈ X . (1.1)
Oft werden die Ablehnbereiche als Intervalle Γα = (cα,∞) konstruiert für kritische Werte
cα = infc ∈ R : sup
ϑ∈Θ0
Pϑ(T (X) > c) 6 α, α ∈ (0, 1). (1.2)
Ist Θ0 = ϑ0 einelementig, dann sind die kritischen Werte genau das (1−α)-Quantil der Ver-teilung von T unter Pϑ0
. Ein wichtiges Konzept in der Testtheorie, insbesondere in Anwendungen,sind die p-Werte.
Denition 1.18. Sei (X ,F , (Pϑ)ϑ∈Θ) ein statistisches Modell und der Test ϕ der HypotheseH0 : ϑ ∈ Θ0 6= ∅ gegeben durch (1.1). Dann ist der p-Wert einer Realisierung x ∈ X bezüglich ϕdeniert als
pϕ(x) = infα:T (x)∈Γα
supϑ∈Θ0
Pϑ(T (X) ∈ Γα).
7
Statt nur zu prüfen, ob ein Test eine Hypothese akzeptiert oder ablehnt, gibt der p-Wert (dieSignikanzwahrscheinlichkeit) das kleinste Signikanzniveau an, zu dem eine Hypothese abgelehntwürde. Damit gibt der p-Wert Aufschluss darüber wie stark die Daten der Hypothese widerspre-chen.
Satz 1.19. Sei (X ,F , (Pϑ)ϑ∈Θ) ein statistisches Modell und sei ϕ ein Test der Hypothese H0 :ϑ ∈ Θ0 6= ∅ gegeben durch ϕ = 1T>cα für eine Teststatistik T : X → R und kritische Werten(cα)α∈(0,1) aus (1.2). Dann ist der p-Wert einer Realisierung x ∈ X bezüglich ϕ gegeben durch
pϕ(x) = supϑ∈Θ0
Pϑ(T (X) > t∗) mit t∗ := T (x).
Sei α ∈ (0, 1) ein fest vorgegebenes Niveau. Ist die Verteilung PTϑ stetig für alle ϑ ∈ Θ0, gilt
ϕ(x) = 1 ⇐⇒ pϕ(x) < α Pϑ − f.s. für alle ϑ ∈ Θ0.
Ist PTϑ (topologisch) diskret verteilt für alle ϑ ∈ Θ0, gilt
ϕ(x) = 1 ⇐⇒ pϕ(x) 6 α Pϑ − f.s. für alle ϑ ∈ Θ0.
Beweis. Deniere P0 := supϑ∈Θ0Pϑ. Da c 7→ P0(T > c) monoton fallend ist, gilt
pϕ(x) = infα:t∗>cα
P0(T > cα) > P0(T > t∗).
Da cα < t∗ äquivalent zur Existenz eines c < t∗ mit P0(T > c) 6 α ist, folgt aus P0(T > cα) 6 α(Rechtsstetigkeit der Verteilungsfunktion), dass
pϕ(x) 6 infα : cα < t∗ 6 infα : P0
(∩c<t∗ T > c
)6 α
= P0(T > t∗).
Zusammen erhalten wir pϕ(x) = P0(T > t∗).Sei nun α fest und T zunächst stetig verteilt. Aus pϕ(x) = P0(T > t∗) < α und limc↑t∗ P0(T ∈
(c, t∗)) = 0 folgt P0(T > c) 6 α für ein c < t∗. Dann muss aber ϕ(x) = 1 gelten. Andersherum gilt
ϕ(x) = 1 ⇒ ∃c < t∗ : P0(T > t∗) 6 α− P0(T ∈ (c, t∗)).
Dabei gilt Pϑ(T ∈ (c, t∗)) = Pϑ(T ∈ (c, T (x))) > 0 für Pϑ-f.a. x ∈ X und für alle ϑ ∈ Θ0. Ist Tdiskret verteilt, bleibt zu bemerken, dass pϕ(x) = P0(T > t∗) = P0(T > c) für ein c < t∗.
Bemerkung 1.20.
(i) Der Vorteil von p-Werten ist, dass sie unabhängig von einem a priori festgesetzten Si-gnikanzniveau α berechnet werden können. Deshalb werden in allen gängigen Statistik-Softwaresystemen statistische Hypothesentests über die Berechnung von p-Werten imple-mentiert.
(ii) Warnung: Alle Rahmenbedingungen des Experiments, insbesondere also das Signikanzni-veau, müssen vor dessen Durchführung festgelegt werden! Ein Signikanzniveau darf nichta posteriori aufgrund der erzielten p-Werte festgelegt werden. Dies widerspricht richtigerstatistischer Praxis! Insbesondere wäre α eine Zufallsvariable (als Funktion in den Beobach-tungen) und obiger Satz kann nicht angewendet werden.
(iii) Der p-Wert gibt eine Antwort auf die Frage: Wie wahrscheinlich sind die gemessenen Daten,gegeben, dass die Nullhypothese stimmt? (und nicht auf die Frage Wie wahrscheinlich istes, dass die Nullhypothese wahr ist, gegeben den gemessenen Daten?)
Beispiel 1.21. Geburten in Berlin:
8
(i) Hypothese: Es werden genauso viele Jungen wie Mädchen geboren.Sind von n ∈ N Geburten w 6 n Mädchen zur Welt gekommen, ist das statistische Modellgegeben durch den Stichprobenraum X = 0, . . . , n und somit (X ,P(X ), (Pϑ)ϑ∈[0,1]) mitBinomialverteilungen Pϑ = Bin(n, ϑ). Die Hypothese führt auf das zweiseitige Testproblem
H0 : ϑ = 1/2 versus H1 : ϑ 6= 1/2,
wobei w ∈ X beobachtet wird. Wir setzten das Niveau α = 0, 05. Die Teststatistik T (w) =|wn − ϑ| führt auf einen zweiseitigen Binomialtest.
(ii) Hypothese: Höchstens die Hälfte der geborenen Kinder hat nicht verheiratete Eltern.Von n ∈ N geboren Kindern haben v 6 n verheiratete Eltern. Mit (X ,P(X ), (Pϑ)ϑ∈[0,1])wie oben betrachten wir hier das einseitige Testproblem
H0 : ϑ 6 1/2 versus H1 : ϑ > 1/2,
wobei v ∈ X beobachtet wird. Das Niveau α = 0, 05 zusammen mit der Teststatistik T (w) =wn − ϑ führt auf einen einseitigen Binomialtest.
Bemerkung 1.22. Bei groÿen Stichprobenumfängen ist es sinnvoll, einen Gauÿ-Test für geeignetnormalisierter Teststatistik zu verwenden, um Binomialtest zu approximieren: Für ϑ ∈ (0, 1)normalisieren wir die Beobachtung X ∼ Bin(n, ϑ) durch Y := X−nϑ√
nϑ(1−ϑ). Aus dem Zentralen
Grenzwertsatz folgt dann für eine standardnormalverteilte Zufallsvariable Z ∼ N (0, 1), dass
Pϑ(T (X) > cα) = Pϑ( |X − nϑ|√
nϑ(1− ϑ)>
√n
ϑ(1− ϑ)cα
)n→∞−→ P
(|Z| >
√n
ϑ(1− ϑ)cα
)=
1
2
(1− Φ
(√ n
ϑ(1− ϑ)cα
))!= α,
Mit der Verteilungsfunktion Φ(x) = P(Z 6 x). Folglich wählen wir cα =√
ϑ0(1−ϑ0)n q1−2α =√
ϑ0(1−ϑ0)n Φ−1(1− 2α) mit ϑ = ϑ0 unter H0.
1.1.3 Kondenzmengen (Bereichsschätzung)
Während ein (Punkt-)Schätzer einen einzelnen Wert angibt, möglichst in der Nähe des wahren Pa-rameters, um Rückschlüsse auf das zugrunde liegende Modell zu ziehen, geben Kondenzbereicheein Intervall an, in dem der Parameter mit gegebener Wahrscheinlichkeit liegt.
Denition 1.23. Sei (X ,F , (Pϑ)ϑ∈Θ) ein statistisches Modell mit abgeleitetem Parameterρ : Θ → Rd. Eine mengenwertige Abbildung C : X → P(Rd) heiÿt Kondenzmenge zumKondenzniveau 1 − α (oder zum Irrtumsniveau α) für α ∈ (0, 1), falls die Messbarkeitsbedin-gung x ∈ X : ρ(ϑ) ∈ C(x) ∈ F für alle ϑ ∈ Θ erfüllt ist und es gilt
Pϑ(ρ(ϑ) ∈ C) = Pϑ(x ∈ X : ρ(ϑ) ∈ C(x)
)> 1− α für alle ϑ ∈ Θ.
Im Fall d = 1 und falls C(x) für jedes x ∈ X ein Intervall ist, heiÿt C Kondenzintervall.
Beachte, dass ρ(ϑ) x ist, während C zufällig ist. Man muss Kondenzmengen also wie folgtinterpretieren: Werden inm unabhängigen Experimenten für (verschiedene) Parameter Kondenz-mengen zum Niveau 0, 95 konstruiert, dann liegt der unbekannte Parameter in 95% der Fälle imder jeweiligen Kondenzmenge (für m groÿ genug; starkes Gesetz der groÿen Zahlen).
Ein verbreitetes Konstruktionsprinzip für die Kondenzintervalle ist die Verwendung einesSchätzers und dessen Verteilung, wie im nächsten Beispiel illustriert.
9
Beispiel 1.24. Im Bernoulli-Experiment von Beispiel 1.5 gilt für Cn := [pn − εn, pn + εn]
Pp(p ∈ Cn) = Pp(|pn − p| < εn) = Pp(∣∣ n∑
i=1
(Xi − p)∣∣ < nεn
) !> 1− α.
Da∑ni=1Xi ∼ Bin(n, p) können wir εn mithilfe der Quantile der Binomialverteilung bestim-
men. Für groÿe n könnte man wieder eine Normalapproximation verwenden. Das resultierendeKondenzintervall besitzt dann aber nur asymptotisch das Niveau 1− α.
Eine alternative Konstruktion von Kondenzmengen bietet folgender Korrespondenzsatz:
Satz 1.25. Sei (X ,F , (Pϑ)ϑ∈Θ) ein statistisches Modell und α ∈ (0, 1). Dann gilt:
(i) Liegt für jedes ϑ0 ∈ Θ ein Test ϕϑ0 der Hypothese H0 : ϑ = ϑ0 zum Signikanzniveau α vor,so deniert C(x) = ϑ ∈ Θ : ϕϑ(x) = 0 eine Kondenzmenge zum Kondenzniveau 1− α.
(ii) Ist C eine Kondenzmenge zum Niveau 1− α, dann ist ϕϑ0(x) = 1− 1C(x)(ϑ0) ein Niveau-
α-Test der Hypothese H0 : ϑ = ϑ0.
Beweis. Nach Konstruktion erhält man in beiden Fällen,
∀ϑ ∈ Θ : ∀x ∈ X : ϕϑ(x) = 0 ⇐⇒ ϑ ∈ C(x).
Damit ist ϕϑ ein Test zum Niveau α für alle ϑ genau dann, wenn
1− α 6 Pϑ(ϕ = 0) = Pϑ(x : ϑ ∈ C(x))
und somit ist C eine Kondenzmenge zum Niveau α.
Beispiel 1.26. Mit Hilfe des Korrespondenzsatzes können wir ein Kondenzintervall zum Niveau0, 95 für die Geburtswahrscheinlichkeit von Mädchen in Berlin berechnen. Im Modell aus Beispiel1.21(i) ist das Kondenzintervall gegeben durch
C(w) = ϑ ∈ [0, 1] : |wn− ϑ| 6 c0,05 = ϑ ∈ [0, 1], pϕ(w) > 0, 05,
wobei pϕ(w) den zu ϕ gehörigen p-Wert der Realisierung w bezeichnet. Ist C sogar ein Konden-zintervall? (Übung 2).
1.2 Minimax- und Bayesansatz
Wir haben bereits verschiedene Schätzmethoden, wie den Maximum-Likelihood-Schätzer oder dieMomentenmethode kennen gelernt. Natürlich gibt es noch viel mehr Konstruktionen. Wie sollteeine Methode anhand des gegeben Schätzproblems ausgewählt werden? Sei also (X ,F , (Pϑ)ϑ∈Θ)ein statistisches Modell mit abgeleitetem Parameter ρ : Θ→ Rd und Verlustfunktion L. Als mögli-ches Vergleichskriterium käme die Risikofunktion R(ϑ, ρ) = Eϑ[L(ϑ, ρ)] eines Schätzers ρ in Frage.Beachte jedoch folgendes Beispiel:
Beispiel 1.27. Sei X ∼ N (µ, 1), µ ∈ R, und L(µ, µ) = (µ − µ)2. Betrachte die zwei Schätzerµ1 = X und µ2 = 5. Die Risiken sind dann gegeben durch
R(µ, µ1) = Eϑ[(X − µ)2] = 1 und R(µ, µ2) = (5− µ)2.
Damit hat µ1 kleineres Risiko als µ2 genau dann, wenn µ /∈ [4, 6].
Denition 1.28. Im statistischen Modell (X ,F , (Pϑ)ϑ∈Θ) mit abgeleitetem Parameter ρ : Θ →Rd und Verlustfunktion L, heiÿt ein Schätzer ρ minimax, falls
supϑ∈Θ
R(ϑ, ρ) = infρ
supϑ∈Θ
R(ϑ, ρ),
wobei sich das Inmum über alle Schätzer (d.h. messbaren Funktionen) ρ : X → Rd erstreckt.
10
Denition 1.29. Der Parameterraum Θ trage eine σ-Algebra FΘ, die Verlustfunktion L sei pro-duktmessbar und ϑ 7→ Pϑ(B) sei messbar für alle B ∈ F . Die a-priori-Verteilung π des Parametersϑ ist gegeben durch ein Wahrscheinlichkeitsmaÿ auf (Θ,FΘ). Das zu π assoziierte Bayesrisiko einesSchätzers ρ ist
Rπ(ρ) := Eπ[R(ϑ, ρ)] =
ˆΘ
ˆXL(ϑ, ρ(x))Pϑ(dx)π(dϑ).
Der Schätzer ρ heiÿt Bayesschätzer oder Bayes-optimal (bezüglich π), falls
Rπ(ρ) = infρRπ(ρ),
wobei sich das Inmum über alle Schätzer (d.h. messbaren Funktionen) ρ : X → Rd erstreckt.
Während ein Minimaxschätzer den maximal zu erwartenden Verlust minimiert, kann das Bayes-risiko als ein (mittels π) gewichtetes Mittel der zu erwartenden Verluste angesehen werden. Al-ternativ wird π als die subjektive Einschätzung der Verteilung des zugrundeliegenden Parametersinterpretiert.
Beispiel 1.27 (fortgesetzt). Oensichtlich kann µ2 kein Minimaxschätzer sein. Zunächst ist esaber nicht klar, ob es einen besseren Schätzer als µ2 gibt. Tatsächlich werden wir später beweisen,dass µ1 minimax ist. Unter der a-priori-Verteilung µ ∼ π = U([4, 6]) hat jedoch µ2 das kleinereBayesrisiko Rπ(µ2) = 1
3 < 1 = Rπ(µ1).
Das Bayesrisiko kann auch als insgesamt zu erwartender Verlust in folgendem Sinne verstandenwerden: Deniere Ω := X × Θ und die gemeinsame Verteilung von Beobachtung und ParameterP auf (X ×Θ,F ⊗FΘ) gemäÿ P(dx, dϑ) = Pϑ(dx)π(dϑ). Bezeichnen X und T die Koordinaten-projektionen von Ω auf X bzw. Θ, dann gilt Rπ(ρ) = EP[L(T, ρ(X))].
Wiederholung: Auf einem Wahrscheinlichkeitsraum (Ω,F ,P) ist die bedingte Wahrscheinlich-keit eines Ereignisses A ∈ F gegeben B ∈ F mit P(B) > 0 deniert als P(A|B) = P(A∩B)/P(B).Sei Ω =
⋃i∈I Bi eine abzählbare Zerlegung in paarweise disjunkte Ereignisse Bi ∈ F , dann besagt
die Bayesformel für jedes A ∈ F mit P(A) > 0 und alle k ∈ I
P(Bk|A) =P(Bk)P(A|Bk)∑i∈I P(Bi)P(A|Bi)
.
Mittels bedingten Erwartungswerten (Stochastik II) kann diese Formel auf Dichten ausgedehntwerden.
Denition 1.30. Sei (X ,F , (Pϑ)ϑ∈Θ) ein von µ dominiertes statistisches Modell mit DichtenfX|T=ϑ := dPϑ
dµ . Sei π eine a-priori-Verteilung auf (Θ,FΘ) mit Dichte fT bzgl. einem Maÿ ν. IstfX|T=· : X × Θ → R+ (F ⊗ FΘ)-messbar, dann ist die a-posteriori-Verteilung des Parametersgegeben der Beobachtung X = x deniert durch die ν-Dichte
fT |X=x(ϑ) =fX|T=ϑ(x)fT (ϑ)´
ΘfX|T=t(x)fT (t)ν(dt)
, ϑ ∈ Θ (PX -f.ü.). (1.3)
Das a-posteriori-Risiko eines Schätzers ρ gegeben X = x ist deniert durch
Rπ(ρ|x) =
ˆΘ
L(ϑ, ρ(x))fT |X=x(ϑ)ν(dϑ).
Beachte, dass im Nenner in (1.3) die Randdichte fX =´
ΘfX|T=t(·)fT (t)ν(dt) bzgl. µ von X
in (X ×Θ,F ⊗FΘ, P) steht, so dass der Nenner in (1.3) für PX -f.a. x ∈ X gröÿer als null ist.
11
Beispiel 1.31. Setze Θ = 0, 1, L(ϑ, r) = |ϑ − r| (0-1-Verlust) und betrachte eine a-priori-Verteilung π mit π(0) =: π0 und π(1) =: π1 = 1 − π0. Die Wahrscheinlichkeitsmaÿe P0
und P1 mögen Dichten p0 und p1 bzgl. einem Maÿ µ besitzen (z.B. µ = P0 + P1). Dann ist diea-posteriori-Verteilung durch die Zähldichte
fT |X=x(i) =πipi(x)
π0p0(x) + π1p1(x), i = 0, 1 (PX -f.ü.)
gegeben. Damit ist das a-posteriori-Risiko eines Schätzers ϑ : X → 0, 1 gegeben durch
Rπ(ϑ|x) =ϑ(x)π0p0(x) + (1− ϑ(x))π1p1(x)
π0p0(x) + π1p1(x).
Satz 1.32. Es gelten die Bedingungen der vorangegangenen Denition. Für das Bayesrisiko einesSchätzers ρ gilt
Rπ(ρ) =
ˆRπ(ρ|x)fX(x)µ(dx).
Minimiert ρ(x) für PX-f.a. das a-posteriori-Risiko mint∈ran(ρ)Rπ(t|x), dann ist ρ Bayesschätzer.
Beweis. Aus (1.3) folgt fT |x=x(ϑ)fX(x) = fX|T=ϑ(x)fT (ϑ). Der Satz von Fubini ergibt
Rπ(ρ) =
ˆΘ
ˆXL(ϑ, ρ(x))Pϑ(dx)π(dϑ)
=
ˆΘ
ˆXL(ϑ, ρ(x))fT |x=x(ϑ)fX(x)µ(dx)ν(dϑ) =
ˆXRπ(ρ|x)µ(dx).
Korollar 1.33. Unter quadratischem Verlust ist der Bayesschätzer gegeben durch
ρ(x) =
ˆΘ
ρ(ϑ)fT |X=x(ϑ)ν(dϑ) =: E[ρ(ϑ)|X = x].
Der Bayesschätzer bzgl. absolutem Verlust ist gegeben durch den Median der a-posteriori-Verteilung. Für den 0-1-Verlust ist der Bayesschätzer der Modus der a-posteriori-Verteilung.
Beweis. Übung 2.
Methode 4: Bayesschätzer. Durch die Wahl einer Verlustfunktion und einer a-priori-Verteilung im statistischen Modell erhalten wir nach Berechnung der a-posteriori-Verteilung unddurch das vorangegangene Korollar einen expliziten Bayesschätzer.
Beispiel 1.34. Sei X1, . . . , Xn ∼ N (µ, σ2) eine mathematische Stichprobe mit bekanntem σ2 > 0und a-priori-Verteilung µ ∼ N (a, b2). Mittels Bayesformel kann die a-posteriori-Verteilung für eineRealisierung x = (x1, . . . , xn) berechnet werden:
fT |X=x(µ) ∼fX|T=µ(x)fT (µ)
∼ exp(−
n∑i=1
(xi − µ)2
2σ2
)exp
(− (µ− a)2
2b2
)∼ exp
(− µ2 − 2µxn
2σ2/n− µ2 − 2aµ
2b2
)∼ exp
(− (b2 + σ2/n)µ2 − 2µ(b2xn + aσ2/n)
2b2σ2/n
)∼ exp
(− 1
2
( nσ2
+1
b2
)(µ− b2xn
b2 + σ2/n− aσ2/n
b2 + σ2/n
)2).
12
Gegeben der Beobachtung X ist ϑ also a-posteriori verteilt gemäÿ
N( b2
b2 + σ²n
Xn −σ²n
b2 + σ²n
a,( nσ2
+1
b2
)−1).
Der Bayesschätzer bzgl. quadratischem Verlust, gegeben durch den a-posteriori Mittelwert, istdamit
ϑn =b2
b2 + σ²n
Xn −σ²n
b2 + σ²n
a.
Bemerkung 1.35. Erhalten wir bei Wahl einer Klasse von a-priori-Verteilungen für ein statisti-sches Modell dieselbe Klasse (i.A. mit anderen Parametern) als a-posteriori-Verteilung zurück, sonennt man die entsprechenden Verteilungsklassen konjugiert. Im obigen Beispiel haben wir ge-sehen, dass die Normalverteilungen zur den Normalverteilungen konjugiert sind (genauer müssteman sagen, dass für unbekannten Mittelwert in der Normalverteilung a-priori Normalverteilungenkonjugiert sind). Als weiteres Beispiel sind die Beta-Verteilungen zur Binomialverteilung konju-giert sind (siehe Übung 2). In diesen (Einzel-)Fällen ist es besonders einfach, die Bayesschätzer zukonstruieren. Für komplexere Modelle werden häug computer-intensive Methoden wie MCMC(Markov Chain Monte Carlo) verwendet, um die a-posteriori-Verteilung zu berechnen (Problem:i.A. hochdimensionale Integration).
Lemma 1.36. Unter den Bedingungen der vorangegangen Denition gilt für jeden Schätzer ρ
supϑ∈Θ
R(ϑ, ρ) = supπRπ(ρ),
wobei sich das zweite Supremum über alle a-priori-Verteilungen π erstreckt. Insbesondere ist dasRisiko eines Bayesschätzers stets kleiner oder gleich dem Minimaxrisiko.
Beweis. Natürlich gilt Rπ(ρ) =´
ΘR(ϑ, ρ)π(dϑ) 6 supϑ∈ΘR(ϑ, ρ). Durch Betrachtung der a-
priori-Verteilung δϑ folgt daher die Behauptung.
Durch dieses Lemma können wir untere Schranken für das Minimaxrisiko durch das Risiko vonBayesschätzern abschätzen. Mögliche Anwendungen illustriert folgender Satz.
Satz 1.37. Sei X1, . . . , Xn eine N (µ, σ2)-verteilte mathematische Stichprobe mit unbekanntemµ ∈ R und bekanntem σ2 > 0. Bezüglich quadratischem Risiko ist das arithmetische Mittel Xn einMinimaxschätzer von µ.
Beweis. Wir betrachten a-priori-Verteilungen µ ∼ π = N (0, b2). Nach Beispiel 1.34 ist die a-posteriori-Verteilung
N( b2Xn
b2 + σ2
n
,( nσ2
+ b−2)−1)
,
der Bayesschätzer bzgl. quadratischem Risiko ist gegeben durch den a-posteriori-Erwartungswertµn = b2Xn/(b
2 + σ2n−1) und dessen a-posteriori-Risiko ist gegeben durch die Varianz der a-posteriori-Verteilung. Ist fX die Randdichte von X von P, folgt aus Satz 1.32
Rπ(µn) =
ˆRn
VarT |X=x(µ)fX(x)dx
=
ˆRn
(nσ−2 + b−2
)−1fX(x)dx =
(nσ−2 + b−2
)−1.
Somit können wir das Minimaxrisiko nach unten abschätzen:
infµ
supµ∈R
R(µ, µ) = infµ
supπRπ(µ) > inf
µsupb>0
RN (0,b2)(µ)
> supb>0
infµRN (0,b2)(µ) = sup
b>0
(nσ2 + b−2
)−1=σ2
n,
wie behauptet, da R(µ,Xn) = σ2/n.
13
1.3 Ergänzungen: Quantile
Denition. Sei P ein Wahrscheinlichkeitsmaÿ auf (R,B(R)) mit Verteilungsfunktion F (x) =P((−∞, x]). Für α ∈ (0, 1) ist das α-Quantil qα ∈ R von P deniert durch
P((−∞, qα)) 6 α 6 P((−∞, qα]).
Die Quantilfunktion ist deniert als verallgemeinertes Inverses von F :
F−1(p) := infx ∈ R : F (x) > p, p ∈ [0, 1].
α-Quantile sind nicht eindeutig, falls F auf dem Niveau α irgendwo konstant ist. Es gilt aber
Lemma. F−1(α) ist ein α-Qunatil.
Beweis. Aufgrund der Rechtsstetigkeit von F gilt F (F−1(α)) > α. Für alle x < F−1(α) giltF (x) < α und wegen der linken Grenzwerte von F
α > limr↑F−1(α)
F (x) = limr↑F−1(α)
P((−∞, r]) = P((−∞, r)).
Das verallgemeinerte Inverse hat folgende Eigenschaften:
(i) F−1(p) 6 x⇔ p 6 F (x);
(ii) F F−1(p) > p und Gleichheit gilt genau dann, wenn p ∈ ranF. Die Gleichheit kann nurdann nicht gelten, wenn F unstetig bei F−1(p) ist;
(iii) F−1 F (x) 6 x, wobei Gleichheit genau dann nicht gilt wenn x im Inneren oder am rechtenRand einer Ebene (kein Anstieg) von F liegt.
Damit gilt F F−1(p) = p auf (0, 1) genau dann, wenn F stetig ist (d.h. ranF = [0, 1]) undF−1 F (x) = x gilt auf R genau dann, wenn F strikt monoton wachsend ist. Folglich ist F−1 einechtes Inverses genau dann, wenn F stetig und streng monoton wachsend ist.
Satz. Ist U ∼ Uni([0, 1]), dann besitzt die Zufallsvariable F−1(U) die Verteilungsfunktion F(Quantilstransformation). Besitzt X die Verteilungsfunktion F , dann gilt F (X) ∼ Uni([0, 1]) ge-nau dann, wenn F stetig ist.
Beweis. Aus (i) folgt P(F−1(U) 6 x) = P(U 6 F (x)) = F (x) für alle x ∈ R. Anderseits gilt fürp ∈ (0, 1) wegen (i) und (ii)
P(F (X) 6 p) = P(X 6 F−1(p)) = F (F−1(p)) = p ⇐⇒ p ∈ ranF.
Schlieÿlich wollen wir noch den QQ-Plot (Qunatil-Quantil-Plot) verstehen: Die empi-rische Verteilungsfunktion einer mathematischen Stichprobe X1, . . . , Xn ist gegeben durchFn(x) = 1
n
∑ni=1 1Xi6x. Die Verteilungsfunktion der Standardnormalverteilung ist Φ(x) =´ x
−∞(2π)−1/2e−y2/2dy. Für groÿe n approximiert Fn die wahre Verteilungsfunktion F , da nach
dem starken Gesetz der groÿen Zahlen Fn(x)→ E[1X16x] = F (x) P-f.s. für alle x ∈ R gilt (tat-sächlich gilt diese Konvergenz sogar gleichmäÿig auf R nach dem Satz von Borel-Cantelli). FallsXi ∼ N (µ, σ2), so gilt F (x) = Φ(x−mσ ). Für die Quantilfunktion gilt also
F−1(Φ(x)) = Φ−1(Φ(x)) · σ +m = σ · x+m,
d.h. F−1 Φ ist eine Gerade. Im QQ-Plot wird F−1n (die empirischen Quantile) gegen Φ−1 aufge-
tragen und unter einer N (µ, σ2)-Annahme sollten die Werte in etwa auf einer Geraden liegen.
14
2 Lineares Modell
2.1 Regression und kleinste Quadrate
Regression ist eine Methode um den Zusammenhang zwischen einer Zielgröÿe (Response-Variable)Y und einem Vektor von erklärenden Variablen (Kovariablen, Regressoren) X = (x1, . . . , xk) zuanalysieren. Beginnen wir mit dem einfachen linearen Modell
Yi = axi + b+ εi, i = 1, . . . , n,
mit Zufallsvariablen ε1, . . . , εn, die zentriert sind (Ei[εi] = 0) und endliche Varianz Var(εi) = σ2 >0 haben. Die Parameter a, b ∈ R, σ > 0 sind unbekannt. Gesucht ist eine Regressionsgerade derForm y = ax+ b, die die Beobachtungen möglichst gut erklärt. Der Parameter σ ist typischerweisenicht das Ziel der statistischen Inferenz und somit ein Störparameter.
Beispiel 2.1. Yi ist das Wachstum von Deutschlands Bruttoinlandsproduktes im Jahr i. DieKovariable xi ist die Veränderung der Arbeitslosenquote im Vergleich zum Vorjahr. Unter Ver-wendung der Daten von 1992 bis 2012 aus den World Development Indicators der Weltbankerhalten als Regressionsgrade erhalten wir y = −1, 080 ·x+1, 338. Betrachten wir alle sechs Grün-dungsmitglieder der EU im gleichen Zeitraum ergibt ganz ähnlich y = −1, 075 · x + 1, 819 Derlineare Zusammenhang beider Gröÿen ist als Okuns Gesetz bekannt.
Um die Situation weiter zu vereinfachen nehmen wir zunächst an, dass ε1, . . . , εn unabhängigund N (0, σ2)-verteilt sind. Nun können wir den Maximum-Likelihood-Schätzer bestimmen: DerBeobachtungsvektor ist verteilt gemäÿ der Lebesguedichte
L(a, b, σ; y) =
n∏i=1
(2πσ2)−1/2 exp(− (yi − axi − b)2
2σ2
)= (2πσ2)−n/2 exp
(− 1
2σ2
n∑i=1
(yi − axi − b)2), y ∈ Rn.
Somit ist die Loglikelihoodfunktion
l(a, b, σ; y) := logL(a, b, σ; y) = −n2
(log σ2 + log(2π))− 1
2σ2
n∑i=1
(yi − axi − b)2.
Das Maximieren der Likelihood über a, b ist also äquivalent zum Minimieren der Summe derquadrierten Residuen (RSS: residual sum of squares). Auch wenn die Fehler nicht normalverteiltsind, kann diese Methode gute Ergebnisse erzielen.
Methode 5: Methode der kleinesten Quadrate. Im einfachen linearen Modell sind dieKleinste-Quadrate-Schätzer a, b durch Minimierung der Summe quadratischen Abstände
(a, b) := arg mina,b
n∑i=1
(Yi − axi − b)2
gegeben.
Satz 2.2. Im einfachen linearen Modell mit unabhängigen und N (0, σ2)-verteilten Fehlern, ist derMaximum-Likelihood-Schätzer gleich dem Kleinste-Quadrate-Schätzer und es gilt
a =
∑ni=1(xi − xn)(Yi − Y n)∑n
i=1(xi − xn)2und b = Y n − axn,
wobei Y n = 1n
∑ni=1 Yi und xn = 1
n
∑ni=1 xi.
15
Beweis. Es bleibt festzustellen, dass wir durch Dierentiation folgende Normalgleichungen erhal-ten:
0 =
n∑i=1
xi(Yi − axi − b) und 0 =
n∑i=1
(Yi − axi − b),
die leicht gelöst werden können.
Bemerkung 2.3. Bei der Wahl anderer Fehlerverteilungen ergibt das Maximum-Likelihood-Prinzipandere (nicht weniger sinnvolle) Schätzer (Übung 2), die aber im Allgemeinen nicht in ge-schlossener Form darstellbar sind. Populäre nicht gauÿsche Fehlerverteilungen sind Laplace- undExponential-Verteilungen.
Haben wir k > 2 Kovariablen und n Beobachtungen Yi, führt das zur multiplen linearenRegression
Yi = β0 +
k∑j=1
βjxi,j + εi, i = 1, . . . , n,
wobei die Fehlerterme (εi) iid. und zentriert sind mit 0 < Var(εi) =: σ2 <∞. In Vektorschreibweiseerhalten wir
Y = (Y1, . . . , Yn)> ∈ Rn Response-Vektor,
X :=
1 x1,1 · · · x1,k
......
...1 xn,1 · · · xn,k
∈ Rn×(k+1) Design-Matrix,
ε := (ε1, . . . , εn)> ∈ Rn Vektor der Fehlerterme,
β := (β0, . . . , βk)> ∈ Rk+1 Parametervektor,
so dass das multiple Regressionsmodell in der Form
Y = Xβ + ε
geschrieben werden kann. Der kleinste-Quadrate-Schätzer löst folglich das Minimierungsproblem
minb|Xb− Y |2.
Beispiel 2.4. Im crime-Datensatz von Agresti and Finlay (1997, Kap. 9) stehen für die 51Staaten der USA die beiden Responsevariablen
Anzahl der Gewaltverbrechen pro 100.000 Einwohnern (crime),
Morde pro 1.000.000 Einwohner (murder),
und folgende Kovariablen zur Verfügung:
Prozentualer Anteil der Bevölkerung die in Ballungs-/ Groÿstadtgebieten leben (pctmetro),
Prozentualer Anteil der weiÿen Bevölkerung (pctwhite),
Prozentualer Anteil der Bevölkerung mit einem High-School-Abschluss (pcths),
Prozentualer Anteil der Bevölkerung der unter der Armutsgrenze leben (poverty) und
Prozentualer Anteil der Bevölkerung mit alleinerziehenden Eltern (single).
16
Bemerkung 2.5. Wechselwirkungen zwischen zwei Kovariablen xi und xj werden durch Interak-tionsterme xi · xj modelliert. Kategorielle Kovariablen sollten durch eine Menge von sogenann-ten Dummy-Indikatoren kodiert werden, um nicht implizit eine (inadäquate) Metrisierung aufdem diskreten Wertebereich solcher Kovariablen zu induzieren. Eine kategorielle Kovariable mit` möglichen Ausprägungen wird dabei durch (` − 1) Indikatoren (d.h. 0, 1-wertige Variablen)repräsentiert. Der j-te Dummy-Indikator kodiert dabei das Ereignis, dass die Kategorie (j + 1)bei der zugehörigen Kovariablen vorliegt, j = 1, . . . , `− 1. Sind also alle (`− 1) Indikatoren gleichNull, so entspricht dies der (Referenz-) Kategorie 1 der zugehörigen kategoriellen Kovariable (vgl.Varianzanalyse).
Dies führt uns zur allgemeinen Denition des linearen Modells:
Denition 2.6. Ein lineares Modell mit n reellwertigen Beobachtungen Y = (Y1, . . . , Yn)> undk-dimensionalem Parameter β ∈ Rk, k < n, besteht aus einer reellen Matrix X ∈ Rn×k vonvollem Rang k, der Designmatrix, und einem Zufallsvektor ε = (ε1, . . . , εn)>, den Fehler- oderStörgröÿen, mit E[εi] = 0, Cov(εi, εj) = Σi,j für eine Kovarianzmatrix Σ > 0. Beobachtet wirdeine Realisierung von
Y = Xβ + ε.
Der (gewichtete) Kleinste-Quadrate-Schätzer β von β minimiert den gewichteten EuklidischenAbstand zwischen Beobachtungen und Modellvorhersage:
|Σ−1/2(Xβ − Y )|2 = infb∈Rk|Σ−1/2(Xb− Y )|2.
Im gewöhnlichen Fall Σ = σ2En mit Fehlerniveau σ > 0 erhalten wir den gewöhnlichen Kleinste-Quadrate-Schätzer (OLS: ordinary least squares)
|Xβ − Y |2 = infb∈Rk|Xb− Y |2,
der unabhängig von der Kenntniss von σ2 ist.
Bemerkung 2.7. Wir schreiben Σ > 0, falls Σ eine symmetrische, strikt positiv-denite Matrixist. Dann ist Σ diagonalisierbar mit Σ = TDT>, D = diag(λ1, . . . , λn) Diagonalmatrix und T
orthogonale Matrix, und wir setzen Σ−1/2 := TD−1/2T> mit D1/2 := diag(λ−1/21 , . . . , λ
−1/2n ).Wie
erwartet, gilt (Σ−1/2)2 = Σ−1 und somit |Σ−1/2v|2 = 〈Σ−1v, v〉.Zusätzlich zur einfachen und multiplen Regression umfasst das lineare Modell weitere Beispiele.
Beispiel 2.8 (Polynomiale Regression). Wir beobachten
Yi = a0 + a1xi + a2x2i + · · ·+ ak−1x
k−1i + εi, i = 1, . . . , n.
Damit ergibt sich als Parameter β = (a0, . . . , ak−1)> und eine Designmatrix vom Vandermonde-Typ
X =
1 x1 x21 · · · xk−1
1...
......
...1 xn x2
n · · · xk−1n
.
Die Matrix hat vollen Rang, sofern k der Designpunkte (xi) verschieden sind.
Lemma 2.9. Setze XΣ := Σ−1/2X. Mit ΠXΣwerde die Orthogonalprojektion von Rn auf den
Bildraum ran(XΣ) bezeichnet. Dann gilt
ΠΣ = XΣ(X>ΣXΣ)−1X>Σ
und für den Kleinste-Quadrate-Schätzer
β = (X>Σ−1X)−1X>Σ−1Y.
Insbesondere existiert der Kleinste-Quadrate-Schätzer, ist eindeutig und erwartungstreu.
17
Beweis. Zunächst beachte, dass X>ΣXΣ = X>Σ−1X invertierbar ist wegen der Invertierbarkeitvon Σ und der Rangbedingung an X:
X>Σ−1Xv = 0⇒ v>X>Σ−1Xv = 0⇒ |Σ−1/2Xv| = 0⇒ |Xv| = 0⇒ v = 0.
Setze PXΣ:= XΣ(X>ΣXΣ)−1X>Σ und w = PXΣ
v für ein v ∈ Rn. Dann folgt w ∈ ran(XΣ) und imFall v = XΣu durch Einsetzen w = PXΣXΣu = v, so dass PXΣ eine Projektion auf ran(XΣ) ist.Da PXΣ selbstadjungiert (symmertrisch) ist, handelt es sich um die Orthogonalprojektion ΠXΣ :
∀u ∈ Rn,∀w ∈ ranXΣ : 〈u− PXΣu,w〉 = 〈u,w〉 − 〈u, PXΣw〉 = 0.
Aus der Eigenschaft β = arg minb |Σ−1/2(Y −Xb)|2 folgt, dass β die beste Approximation vonΣ−1/2Y durch XΣb liefert. Diese ist durch die Orthogonalprojektionseigenschaft ΠXΣΣ−1/2Y =
XΣβ bestimmt. Es folgt
X>Σ ΠXΣΣ−1/2Y = (X>ΣXΣ)β ⇒ (X>ΣXΣ)−1X>Σ−1Y = β.
Schlieÿlich folgt aus der Linearität des Erwartungswertes und E[ε] = 0:
E[β] = E[(X>ΣXΣ)−1X>Σ−1(Xβ + ε)] = β + 0 = β.
Bemerkung 2.10.
Im gewöhnlichen linearen Modell bzw. der multiplen linearen Regression gilt β =(X>X)−1X>Y und ist somit unabhängig vom unbekannten Parameter σ > 0.
X†Σ := (X>ΣXΣ)−1X>Σ heiÿt auch Moore-Penrose-(Pseudo-)Inverse von XΣ, so dass β =
X†ΣΣ−1/2Y bzw. β = X†Y im gewöhnlichen linearen Modell gilt.
Wir kommen zum zentralen Satz in der Regressionsanalyse:
Satz 2.11 (Gauÿ-Markov). Ist der Parameter ρ = 〈β, v〉 für ein v ∈ Rk im linearen Modell
zu schätzen, so ist ρ = 〈β, v〉 ein (in den Daten Y ) linearer erwartungstreuer Schätzer, derunter allen linearen erwartungstreuen Schätzern minimale Varianz besitzt, nämlich Var(ρ) =|XΣ(X>ΣXΣ)−1v|2.
Beweis. Die Linearität ist klar und aus dem vorangegangen Lemma folgt, dass ρ erwartungstreuist. Sei nun ρ = 〈Y,w〉 ein beliebiger linearer erwartungstreuer Schätzer von ρ. Dies impliziert füralle β ∈ Rk
E[〈Y,w〉] = ρ⇒ 〈Xβ,w〉 = 〈β, v〉 ⇒ 〈X>w − v, β〉 = 0
und somit v = X>w = X>Σ Σ1/2w. Nach Pythagoras erhalten wir
Var(ρ) = E[〈ε, w〉2] = E[w>εε>w]
= w>Σw = |Σ1/2w|2 = |ΠXΣ(Σ1/2w)|2 + |(En −ΠΣ)(Σ1/2w)|2.
Damit gilt Var(ρ) > |ΠXΣ(Σ1/2w)|2 = |XΣ(X>ΣXΣ)−1X>w| = |XΣ(X>ΣXΣ)−1v| = Var(ρ).
Bemerkung 2.12. Man sagt, dass der Schätzer ρ im Satz von Gauÿ-Markov bester linearererwartungstreuer Schätzer (blue: best linear unbiased estimator) ist. Eingeschränkt auf lineareSchätzer ist der Kleinste-Quadrate-Schätzer damit minimax. Ob es einen besseren nichtlinearenSchätzer geben kann, werden wir in Kapitel 3 beantworten.
Im gewöhnlichen linearen Modell ist die optimale Varianz insbesondere σ2|X(X>X)−1v|2. Indiesem Spezialfall ist es auch von Interesse das Rauschniveau σ2 zu schätzen. Dies ermöglicht esinsbesondere Tests und Kondenzbereiche zu konstruieren.
18
Lemma 2.13. Im gewöhnlichen linearen Modell mit σ > 0 und Kleinste-Quadrate-Schätzer β giltXβ = ΠXY und R := Y − Xβ bezeichne den Vektor der Residuen. Die geeignet normalisierteStichprobenvarianz
σ2 :=|R|2
n− k=|Y −Xβ|2
n− kist erwartungstreuer Schätzer von σ2.
Beweis. Xβ = ΠXY folgt aus Lemma 2.9. Einsetzen zeigt E[|Y − Xβ|2] = E[|Y − ΠXY |2] =E[|(En − ΠX)ε|2]. Ist nun e1, . . . , en−k eine Orthnormalbasis vom (n − k)-dimensionalen Bildran(En −ΠX) ⊆ Rn, so folgt
E[|(En −ΠX)ε|2] =
n−k∑i=1
E[〈ε, ei〉2] = σ2(n− k),
was die Behauptung impliziert.
Beachte, dass der Maximum-Likelihood-Schätzer von σ2 gegeben ist durch σ2ML = n−1|R|2 6=
σ2 (Übung 2). Der erwartungstreue Schätzer σ2 wird in der Praxis bevorzugt, hat jedoch gröÿereVarianz als σ2
ML.Bevor wir uns mit statistischer Inferenz, also der Konstruktion von Tests und Kondenzinter-
vallen, im linearen Modell beschäftigen, soll der Baysianische Ansatz auf das Regressionsproblemangewendet werden.
Satz 2.14. Im gewöhnlichen linearen Modell Y = Xβ + ε mit ε ∼ N (0, σ2En) und bekanntemσ > 0 genüge β ∈ Rk der a-priori-Verteilung
β ∼ N (m,σ2M)
mit Parametern m ∈ Rk und symmetrisch positiv deniter Matrix M ∈ Rk×k. Dann ist die a-posteriori-Verteilung von β gegeben einer Realisierung y ∈ Rn gegeben durch
β|Y = y ∼ N (µy,Σy) mit Σy = σ2(X>X +M−1
)−1, µy = Σy(σ−2X>y + σ−2M−1m).
Insbesondere ist der Bayesschätzer bzgl. quadratischem Verlust gegeben durch βBayes =(X>X +
M−1)−1
(X>Y +M−1m).
Beweis. Für die a-posteriori-Dichte an der Stelle t ∈ Rk gilt
fβ|Y=y(t) ∼ exp(− 1
2σ2(y −Xt)>(y −Xt)
)exp
(− 1
2σ2(t−m)>M−1(t−m)
)∼ exp
( 1
σ2t>X>y − 1
2σ2t>X>Xt− 1
2σ2t>M−1t+
1
σ2t>M−1m
)= exp
( 1
σ2t>(X>y +M−1m
)− 1
2σ2t>(X>X +M−1
)t).
Daher ist β gegeben Y = y normalverteilt mit Kovarianzmatrix Σy = (σ−2X>X + σ−2M−1)−1
und Mittelwert µy = Σy(X>y +M−1m)/σ2.
Es ist erneut bemerkenswert, dass der Bayesschätzer βBayes nicht von σ2 abhängt.
Bemerkung 2.15. Indem wir auch den Parameter σ2 mit einer a-priori-Verteilung versehen, er-halten wir ein (mehrstuges) Bayesmodell. Da wir besonders an konjugierten Verteilungsklasseninteressiert sind, wird hierzu oft die inverse Gamma-Verteilung verwendet: Ist Z ∼ Γ(a, b) so ist1/Z ∼ IG(a, b) invers Gamma-verteilt mit Parametern a, b > 0 und Lebesguedichte
fa,b(x) =ba
Γ(a)x−(a−1)e−a/x1(0,∞)(x), x ∈ R.
19
Das Bayesmodell ist also gegeben durch
Y |β, σ2 ∼ N (Xβ, σ2En), β|σ2 ∼ N (m,σ2M), σ ∼ IG(a, b).
Die gemeinsame Verteilung von (β, σ2) ∼ NIG(m,M, a, b) wird Normal-inverse Gammaverteilunggenannt und besitzt die Dichte
f(β, σ2) =1
(2πσ2)k/2|M |1/2exp
(− 1
2σ2(β −m)>M−1(β −m)
) ba
Γ(a)(σ2)a+1e−a/σ
2
∼ 1
(σ2)k/2+a+1exp
( 1
2σ2
((β −m)>M−1(β −m) + b
)), β ∈ Rk, σ2 > 0.
In diesem Modell ist die a-posteriori-Verteilung von σ2 gegeben β und Y gegeben durch σ2|β, Y ∼IG(a′, b′) mit a′ = a+ n
2 + k2 und
b′ = b+1
2(Y −Xβ)>(Y −Xβ) +
1
2(β −m)>M−1(β −m).
Die a-posteriori-Verteilung von (β, σ2) gegeben Y ist (β, σ2)|Y ∼ NIG(m, M , a, b) mit Parametern
M = (X>X +M−1)−1, m = M(M−1m+X>y),
a = a+n
2, b = b+
1
2
(Y >Y +m>M−1m− m>M−1m
),
siehe Fahrmeir et al. (2009, Kap. 3.5).
Korollar 2.16. Unter den Voraussetzungen des vorangegangenen Satzes mit m = 0 und M =τ2Ek, τ > 0, gilt für den Bayesschätzer unter quadratischem Verlust
βBayes = arg minβ∈Rk
|Y −Xβ|2 +1
τ2|β|2.
Beweis. Im Spezialfall m = 0 und M = τ2Ek folgt aus obigem Satz βBayes =(X>X +
τ−2Ek)−1
X>y. Andererseits gilt
arg minβ
((Y > − β>X>)(Y −Xβ) +
1
τ2β>β
)= arg min
β
(− 2Y >Xβ + β>(X>X +
1
τ2Ek)β
).
Null setzen des Dierenzials der Funktion β 7→ −2Y >Xβ + β>(X>X + 1τ2Ek)β liefert 0 =
−2Y >X+2β>(X>X+ 1τ2Ek), so dass aus der positiv Denitheit und Symmetrie von X>X+ 1
τ2Ekdie Behauptung folgt.
Der Bayesansatz führt uns also zu einer neuen Schätzmethode im linearen Modell:
Methode 6: Ridge-Regression. Im linearen Modell Y = Xβ + ε ist derRidge-Regressionsschätzer oder Schrumpfungsschätzer (engl.: Shrinkage) mit Schrumpfungs-koezient λ > 0 deniert als
βridge = arg minβ∈Rk
|Y −Xβ|2 + λ|β|2.
Durch Einführung des Strafterms (engl.: penalty) λ|β|2 wird die Varianz auf Kosten eines Biasverringert. Dies ist insbesondere sinnvoll, wenn einige (wenige) Koezienten von β groÿ sind unddie übrigen klein und liefert in diesen Fällen gute Schätzergebnisse auch wenn die Parameterdi-mension in einer ähnlichen Gröÿenordnung liegt wie die Anzahl der Beobachtungen (n ∼ p). Dieswird im nächsten Beispiel illustriert. Die richtige Wahl des Strumpfungsparameters λ ist allerdingsein schwieriges Problem.
20
Beispiel 2.17. Betrachten wir das Modell Yi = x>i β + εi mit Kovariablenvektor xi ∈ Rp, Pa-rameter β ∈ Rp und εi
iid.∼ N (0, 1) mit i = 1, . . . , n. Wir wählen n = 50 und p = 30 wobei 10Koezienten groÿ sind (zwischen 0,5 und 1) und 20 klein (zwischen 0 und 0,3) und bestimmenden mittleren Quadratischen Fehler aus 200 Simulationen für verschiedene Werte von λ ∈ [0, 20](Übung 2).
2.2 Inferenz unter Normalverteilungsannahme
Im Folgenden werden wir das gewöhnliche lineare Modell unter der Normalverteilungsannahme(εi) ∼ N (0, σ2En) betrachten.
Beispiel 2.18. Sind die Messfehler (εi) ∼ N (0, σ2En) gemeinsam normalverteilt und ρ = 〈v, β〉für v ∈ Rk, so gilt
β ∼ N (β, σ2(X>X)−1) und ρ = 〈v, β〉 ∼ N (γ, σ2v>(X>X)−1v).
Ist σ > 0 bekannt, so ist ein Kondenzintervall zum Niveau 95% für ρ gegeben durch
I0,95(ρ) :=[ρ− 1, 96σ
√v>(X>X)>v, ρ+ 1, 96σ
√v>(X>X)>v
].
Dabei ist der Wert 1, 96 gerade das 0,975-Qunatil bzw. 0,025 Fraktil der Standardnormalverteilung.Analog (Korrespondenzsatz) wird der zweiseitige Gauÿ-Test der Hypothese H0 : ρ = ρ0 gegenH1 : ρ 6= ρ0 zum Niveau α ∈ (0, 1) konstruiert: Wähle die Teststatistik |ρ− ρ0| und den kritischenWert q1−α/2σ
√v>(X>X)>v mit dem (1− α/2)-Quantil von N (0, 1).
Ist σ unbekannt, so ist eine Idee, einfach σ durch den Schätzer σ in obiger Formel zu ersetzen.Allerdings wird dann das vorgegebene Niveau nur noch asymptotisch erreicht für einen konsistentenSchätzer (Slutsky-Lemma). Im vorliegenden Fall können wir aber sogar die Verteilung für endlicheStichprobenumfänge exakt bestimmen.
Denition 2.19. Die t-Verteilung t(n) (oder Student-t-Verteilung) mit n ∈ N Freiheitsgradenauf (R,B(R)) ist gegeben durch die Lebesguedichte
tn(x) =Γ(n+1
2 )
Γ(n2 )√πn
(1 +
x2
n
)−(n+1)/2
, x ∈ R.
Die F-Verteilung F (m,n) (oder Fisher-Verteilung) mit (m,n) ∈ N2 Freiheitsgraden auf (R,B(R))ist gegeben durch die Lebesguedichte
fm,n(x) =mm/2nn/2
B(m2 ,n2 )
xm/2−1
(mx+ n)(m+n)/21R+(x), x ∈ R.
Dabei bezeichnet Γ(p) =´∞
0tp−1e−tdt die Gamma-Funktion und B(p, q) = Γ(p)Γ(q)
Γ(p+q) die Beta-Funktion.
Erinnerung: Für X1, . . . , Xm ∼ N (0, 1) ist X :=∑mi=1X
2i ∼ χ2(n) verteilt mit Lebesguedichte
fX(x) = (2m/2Γ(m2 ))−1xm/2−1e−x/21R+(x).
Lemma 2.20. Es seien X1, . . . , Xm, Y1, . . . , Yn unabhängige N (0, 1)-verteilte Zufallsvariablen.Dann gilt
Tn :=X1√
1n
∑nj=1 Y
2j
∼ t(n) und Fm,n :=1m
∑mi=1X
2i
1n
∑nj=1 Y
2j
∼ F (m,n).
Beweis. Es gilt T 2n = F1,n, so dass mittels Dichtetransformation f|Tn|(x) = fF1,n
(x2)2x, x > 0,gilt. Da Tn symmetrisch (wie −Tn) verteilt ist, folgt fTn = FF1,n
(x2)|x|, x ∈ R, und Einsetzenzeigt die Behauptung für Tn, sofern F1,n F (1, n)-verteilt ist.
21
Um die Behauptung für Fm,n nachzuweisen, benutze, dass X :=∑mi=1X
2i χ
2(m)-verteilt undY :=
∑nj=1 Y
2j χ2(n)-verteilt sind. Wegen Unabhängigkeit von X und Y gilt für z > 0 (setze
w = x/y)
P(X/Y 6 z) =
ˆ ˆ1x/y6zfx(x)fY (y)dxdy
=
ˆ1w6z
(ˆfX(wy)fY (y)ydy
)dw,
so dass sich die Dichte wie folgt ergibt (setze w = (z + 1)y)
fX/Y (z) =
ˆfX(zy)fY (y)ydy
=2−(m+n)/2
Γ(m2 )Γ(n2 )
ˆ ∞0
(zy)m/2−1yn/2e−(zy+y)/2dy
=2−(m+n)/2
Γ(m2 )Γ(n2 )
ˆ ∞0
(zw/(z + 1))m/2−1(w/(z + 1))n/2e−w/2(z + 1)−1dw
=Γ(m+n
2 )
Γ(m2 )Γ(n2 )zm/2−1(z + 1)−(m+n)/2, z > 0.
Dichtetransformation ergibt damit für Fm,n = nmXY die Dichte m
n fX/Y (mn x) = fm,n(x).
Bemerkung 2.21. Es gilt T 2n = F1,n. Für n = 1 ist die t(n)-Verteilung gerade die Cauchy-Verteilung
und für n → ∞ konvergiert sie schwach gegen die Standardnormalverteilung. Für jedes n ∈ Nbesitzt t(n) nur Momente bis zur Ordnung p < n (sie ist heavy-tailed). Ähnliches gilt für dieF-Verteilung, insbesondere konvergiert die Verteilung von mFm,n für n → ∞ gegen die χ2(m)-Verteilung.
Aus diesem Lemma ergeben sich die Standardtests für die Parameter der Normalverteilung,siehe Witting (1985, S. 200-204).
Bevor wir zur Konstruktion von Tests und Kondenzbändern im linearen Modell kommen nochein weiteres nützliches Hilfsresultat zur Verteilung quadratischer Formen:
Lemma 2.22. Seien X ∼ N (0, En) und R eine symmetrische, idempotente (n× n)-Matrix (d.h.R = R> und R2 = R) mit rank(R) = r 6 n. Dann gilt
(i) X>RX ∼ χ2(r),
(ii) X>RX ist unabhängig von BX für jede Matrix B ∈ Rp×n mit p 6 n und BR = 0,
(iii) für jede weitere symmetrische, idempotente Matrix S ∈ Rn×n mit rank(S) = s 6 n undRS = 0 sind X>RX und X>SX unabhängig und
s
r
X>RX
X>SX∼ F (r, s).
Beweis. (i) Da R symmetrisch und idempotent ist, existiert eine Orthogonalmatrix P mit R =
PDrP>, wobei Dr =
(Er 00 0
). Da P orthogonal ist und X standardnormalverteilt, folgt W :=
P>X ∼ N (0, En). Wegen
X>RX =X>R2X = (RX)>(RX) = (PDrW )>(PDrW ) = W>DrW =
r∑i=1
W 2i
ist X>RX χ2(r)-verteilt.
22
(ii) Wir setzen Y := BX ∼ N (0, B>B) und Z := RX ∼ N (0, R). Dann gilt
Cov(Y,Z) = BVar(X)R> = BR = 0.
Da (X,Y ) als Lineartransformation von X gemeinsam normalverteilt ist, folgt aus der Unkorre-liertheit bereits die Unabhängigkeit.
(iii) Genau wie in (ii) folgt die Unabhängigkeit von Y := SX und Z := RX und somitauch die Unabhängigkeit von Y >Y = X>SX und Z>Z = X>RX. Zusammen mit (i) und demvorangegangenen Lemma folgt die Behauptung.
Als Korollar erhalten wir Kondenzbereiche für die Schätzung von β und linearen Funktionalenim gewöhnlichen linearen Modell unter der Normalverteilungsannahme.
Satz 2.23. Im gewöhnlichen linearen Modell unter der Normalverteilungsannahme (εi) ∼N (0, σ2En) für σ > 0 gelten folgende Kondenzaussagen für gegebenes Niveau α ∈ (0, 1) :
(i) Ist qF (k,n−k);1−α das (1− α)-Quantil der F (k, n− k)-Verteilung, so ist
C :=β ∈ Rk
∣∣|X(β − β)|2 < kσ2qF (k,n−k);1−α
ein Kondenzellipsoid zum Kondenzniveau 1− α für β.
(ii) Ist qt(n−k);1−α/2 das (1− α2 )-Quantil der t(n− k)-Verteilung, so ist
I :=[ρ− σ
√v>(X>X)−1vqt(n−k);1−α/2, ρ+ σ
√v>(X>X)−1vqt(n−k);1−α/2
]ein Kondenzintervall zum Kondenzniveau 1− α für ρ = 〈v, β〉.
Beweis. (i) Nach Konstruktion gilt
Xβ = XX†Y = ΠXY = Xβ + ΠXε, σ2 =|(En −ΠX)ε|2
(n− k).
Da ΠX und (En −ΠX) symmetrische, idempotente Matrizen mit Rang k bzw. (n− k) sind (Pro-jektionen auf ranX bzw. (ranX)⊥) und es gilt (En + ΠX)ΠX = 0, folgt aus Lemma 2.22:
|X(β − β)|2
kσ2=
(n− k)
k
ε>(En −ΠX)ε
ε>ΠXε∼ F (k, n− k).
Durch die Wahl des Quaniles folgt die Kondenzaussage Pβ(β ∈ C) = 1− α.(ii) Wegen ρ ∼ N (ρ, σ2v>(X>X)−1v) nach dem Satz von Gauÿ-Markov, ist
ρ− ρσ√v>(X>X)−1v
∼ N (0, 1).
Andererseits sind ρ und σ2 unabhängig und es gilt σ2 = σ2Z/(n − k) für eine ZufallsvariableZ ∼ χ2(n− k). Damit ist
ρ− ρ√σ2v>(X>X)−1v
∼ t(n− k).
Bemerkung 2.24. Ebenso kann man ein Kondenzintervall für die Varianz konstruieren (Übung2).
Zusammen mit dem Korrespondenzsatz liefert dieses Resultat:
23
Methode 7: t-Test und F-Test. Im gewöhnlichen linearen Modell unter Normalverteilungs-annahme (εi) ∼ N (0, σ2En) ist der (zweiseitige) t-Test der Hypothese H0 : ρ = ρ0 gegen dieAlternative H1 : ρ 6= ρ0 für ρ0 = 〈v, β0〉 zum Niveau α ∈ (0, 1) gegeben durch
ϕρ0(Y ) = 1|Tn−k(Y )|>qt(n−k);1−α/2 mit Tn−k(Y ) :=ρ0 − ρ
σ√v>(X>X)−1v
.
Der F-Test der Hypothese H0 : β = β0 vs. H1 : β 6= β0 zum Niveau α ∈ (0, 1) ist gegeben durch
ϕβ0(Y ) = 1Fk,n−k(Y )>qF (k,n−k);1−α mit Fk,n−k(Y ) :=|X(β0 − β)|2
kσ2.
Schlieÿlich wollen wir Hypothesentests noch für den allgemeineren Fall von linearen (bzw.anen) Hypothesen konstruieren.
Denition 2.25. Im gewöhnlichen linearen Modell ist ein (zweiseitiges) lineares Testproblemgegeben durch
H0 : Kβ = d versus H1 : Kβ 6= d
für eine (deterministische) Matrix K ∈ Rr×k mit vollem Rang rank(K) = r 6 k und einem Vektord ∈ Rr. K wird Kontrastmatrix genannt. Unter der Hypothese H0 sind also insgesamt r 6 k linearunabhängige Bedingungen an die Parameter des linearen Modells gestellt.
Beispiel 2.26. Test auf Gleichheit zweier Regressionskoezienten: Für 2 6 j < l 6 k ist dasTestproblem gegeben durch
H0 : βj = βl versus H1 : βj 6= βl.
Damit ist die Kontrastmatrix K = (a1,i) ∈ R1×k gegeben durch a1,i = 1i=j − 1i=l und d = 0.Weitere Beispiele sind der Globaltest (Übung 2):
H0 : ∀j ∈ 1, . . . , k : βj = 0 versus H1 : ∃j ∈ 1, . . . , k : βj 6= 0
sowie der Test eines Subvektors β∗ = (β∗1 , . . . , β∗r )> mit r 6 k (Übung 2):
H0 : ∀j ∈ 1, . . . , r : βj = β∗j versus H1 : ∃j ∈ 1, . . . , r : βj 6= β∗j .
Die Grundidee für das Testen linearer Hypothesen ist, die Residuen RSS = |Y − Xβ|2 desKleinste-Quardrate-Schätzers mit den Residuen des auf H0 : Kβ = d eingeschränkten Kleinste-Quardate-Schätzers βH0 , d.h.
RSSH0 := |Y −XβH0 |2 mit |Y −XβH0 |2 = minβ∈Rk:Kβ=d
|Y −Xβ|2,
zu vergleichen. Ist die Abweichung (relativ zu RSS) zu groÿ, spricht dies gegen die Hypothese.
Satz 2.27. Im gewöhnlichen linearen Modell unter Normalverteilungsannahme (εj) ∼ N (0, σ2En)ist die lineare Hypothese
H0 : Kβ = d versus H1 : Kβ 6= d
mit Kontrastmatrix K ∈ Rr×k und d ∈ Rr zu testen. Es gilt
(i) βH0= β − (X>X)−1K>(K(X>X)−1K>)−1(Kβ − d),
(ii) RSSH0−RSS = (Kβ − d)>(K(X>X)−1K>)−1(Kβ − d) und (RSSH0
−RSS)/σ2 ∼ χ2(r)unter H0
(iii) die Fisher-Statistik F := n−kr
RSSH0−RSS
RSS ist unter H0 gemäÿ F (r, n− k) verteilt.
24
Beweis. (i) Für jeden Vektor γ ∈ Rk, der die Nebenbedingung Kγ=d erfüllt, gilt
|Y −Xγ|2 = |Y −Xβ +X(β − γ)|2 = |Y −Xβ|2 + |X(β − γ)|2
nach Pythagoras, da Y −Xβ = (En −ΠX)Y ⊥ ran(X). Auÿerdem ist
|X(β − γ)|2 = |X(β − βH0)|2 + |X(βH0
− γ)|2 + 2〈X(β − βH0), X(βH0
− γ)〉.
Die Wahl von βH0 impliziert jedoch
〈X(β − βH0), X(βH0
− γ)〉 =((X>X)−1K>(K(X>X)−1K>)−1(Kβ − d)
)>X>X(βH0
− γ)
=(Kβ − d)>(K(X>X)−1K>)−1(KβH0−Kγ) = 0,
denn βH0 erfüllt die Nebenbedingung:
KβH0= Kβ −K(X>X)−1K>(K(X>X)−1K>)−1(Kβ − d) = d.
Insgesamt erhalten wir also
|Y −Xγ|2 = |Y −Xβ|2 + |X(β − βH0)|2 + |X(βH0 − γ)|2, (2.1)
was oensichtlich für γ = βH0minimal ist.
(ii) Aus (2.1) mit γ = βH0folgt durch Einsetzen von βH0
RSSH0−RSS =|Y −XβH0
|2 − |Y −Xβ|2 = |X(β − βH0)|2
=(β − βH0)>X>X(β − βH0
)
=(Kβ − d)>(K(X>X)−1K>)−1(Kβ − d).
Unter H0 gilt für die Zufallsvariable Z := Kβ, dass E[Z] = d und Var(Z) = σ2K(X>X)−1K>.Aus der Normalverteilung von β folgt daher (RSSH0 −RSS)/σ2 ∼ χ2(r).
(iii) Da RSSH0− RSS eine Funktion von β ist und somit unabhängig von RSS ist
(Lemma 2.22), folgt die Verteilungsaussage für F aus der Charakterisierung der F (r, n − p)-Verteilung.
Bemerkung 2.28. W := rF heiÿt auchWald-Statistik. Im Fall d = 0 ist L := Xβ|β ∈ Rk,Kβ = 0ein linearer Unterraum von ranX und XβH0
= ΠLY die Orthogonalprojektion der BeobachtungenY auf L. In diesem Fall gilt nach Pythagoras
RSSH0= |Y −ΠLY |2 = |Y −ΠXY + (ΠX −ΠL)Y |2 = |Y −ΠXY |2 + |Xβ −XβH0
|2,
so dass die Fisher-Statistik auch als
F =|Xβ −XβH0
|2
rσ2
geschrieben werden kann.
Beispiel 2.26 (fortgesetzt). Einsetzen von K und d liefert
F =n− kRSS
(βj − βl)2
K(X>X)−1K>.
Wegen Var(βj − βl) = Var(Kβ) = σ2K(X>X)−1K> ist Var(βj-βl) = K(X>X)−1K>σ2 mitσ2 = RSS/(n− k) der natürliche (plug-in) Varianzschätzer. Damit können wir die Test-StatistikF als
F =(βj − βl)2
Var(βj − βl)H0∼ F (1, n− k)
25
schreiben. Dieser F-Test ist äquivalent zum (zweiseitigen) t-Test mit der Teststatistik
T =βj − βl
(Var(βj − βl))1/2∼ t(n− k).
Beispiel 2.29 (Klimaentwicklung). Wir folgen Beispiel 12.24 von Georgii (2007) undbetrachten die mittleren Augusttemperaturen von 1799 bis 2008 in Karlsruhe (Quelle:http://www.klimadiagramme.de/Europa/special01.htm). Für die Jahre 1854 und 1945 liegenkeine Daten vor, so dass wir n = 208 Beobachtungen haben. Eine polynomielle Regression in derZeit t (in Jahrhunderten beginnend bei 1799) mit Graden d = 1, . . . , 4 liefert
p1(t) = 18, 7 + 0, 1t,
p2(t) = 20, 0− 3, 5t+ 1, 7t2,
p3(t) = 19, 5− 0, 6t− 1, 7t2 + 1, 1t3,
p4(t) = 19, 4 + 0, 5t− 4, 1t2 + 2, 9t3 − 0, 4t4.
Zunächst ist es plausibel, dass die zufälligen Schwankungen unabhängig von einander sind und alsnäherungsweise normalverteilt angenommen werden können (QQ-Plot). Um statistisch verwert-bare Aussagen zu treen, setzen wir noch das Niveau α = 0, 05 fest. Der Parametervektor istβ = (β0, . . . , βd)
>. Welcher Grad des Regressionspolynoms ist sinnvoll?Frage 1 : Ist der positive Trend von p1 signikant? H0 : β1 6 0 vs. H1 : β1 > 0. Die zugehörige
t-Statistik T = β1
σ√v>(X>X)−1v
≈ 0, 62 liegt deutlich unter dem kritischen Wert qt(n−2),1−α ≈ 1, 65
(einseitiger T-Test), so dass die Hypothese nicht verworfen werden kann.Frage 2 : Liegt den Beobachtungen ein linearer Zusammenhang zugrunde (im Modell mit d =
4)? H0 : β2 = β3 = β4 = 0. Mittels Bemerkung 2.28 berechnen wir die Fisher-Statistik
F =
∑nk=1(p4(tk)− p1(tk))2
3σ2≈ 13, 68 > 2, 65 ≈ qF (3,n−5),1−α.
Folglich kann die Hypothese abgelehnt werden und wir schlussfolgern, dass eine Regressionsgeradeunzureichend ist.
Frage 3 : Benötigen wir ein Polynom vierten Grades?H0 : β4 = 0. Die zugehörige t-Statistik hatden Wert −0, 41 dessen Absolutbetrag kleiner als das Quantil qt(n−5),0.975 ≈ 1, 97 ist (zweiseitigert-Test). Diese Nullhypothese kann also akzeptiert werden.
Frage 4: Benötigen wir ein Polynom dritten Grades ? H0 : β3 = 0 (im Modell mit d =3). Die zugehörige t-Statistik hat den Wert 2, 05 dessen Absolutbetrag gröÿer als das Quantilqt(n−4),0.975 ≈ 1, 97 ist. Die Hypothese kann also abgelehnt werden und der kubische Anteil imRegressionspolynom ist signikant, d.h. p3 ist signikant besser geeignet die Beobachtungen zubeschreiben als p2.
p3 zeigt einen deutlichen Anstieg der Temperaturen im 19. Jahrhundert. Es sei bemerkt, dasswir hier nur eine Zeitreihe betrachtet haben und somit nicht auf einen allgemeinen Zusammenhangschlieÿen können (Aufgabe der Klimatologen).
2.3 Varianzanalyse
Beispiel 2.30. Um den Einuss von k ∈ N verschiedenen Düngemitteln auf den Ernteertrag zuvergleichen wird jedes Düngemittel i ∈ 1, . . . , k auf ni verschiedenen Agrarächen ausgebracht.Der durch Witterungseinüsse etc. zufällige Ernteertrag kann mittels Yij = µi+εij für j = 1, . . . , niund i = 1, . . . , k modelliert werden, wobei µi der mittlere Ernteertrag von Düngemittel i ist undεij unabhängige, zentrierte Störgröÿen sind. Wir fragen uns also ob µ1 = · · · = µk gilt oder nicht.
Denition 2.31. Das Modell der einfaktoriellen Varianzanalyse (ANOVA1: (one-way) analysisof variance) ist gegeben durch Beobachtungen
Yij = µi + εij , i = 1, . . . , k, j = 1, . . . , ni,
26
mit iid.-verteilten Störgröÿen εij ∼ N (0, σ2). Wir bezeichnen die erste Dimension als den Faktorund den Wert i = 1, . . . , k als die Faktorstufe. Folglich geben (ni)i=1,...,k die Anzahl der unabhän-gigen Versuchswiederholungen pro Faktor an und n :=
∑ki=1 ni ist der Gesamtstichprobenumfang.
Gilt n1 = · · · = nk, so sprechen wir von balanciertem Design.
Damit ist das ANOVA1-Modell ein Spezialfall des gewöhnlichen linearen Modells der Form
Rn 3 Y :=
Y11
...Y1n1
...Yk1
...Yknk
=
1 0 · · · 0...
......
1 0 · · · 0...
......
0 0 · · · 1...
...0 0 · · · 1
︸ ︷︷ ︸
=:X∈Rn×k
·
µ1
...µk
︸ ︷︷ ︸=:µ∈Rk
+
ε11
...ε1n1
...εk1
...εknk
.
Beachte, dass rankX = k. Die klassische Fragestellung der Varianzanalyse lautet: ExistierenUnterschiede in den Faktorstufen-spezischen Mittelwerten µi? oder anders formuliert Hat derFaktor einen Einuss auf die Response oder nicht?. Dies führt auf das Testproblem
H0 : µ1 = · · · = µk versus H1 : ∃i, l ∈ 1, . . . , k : µi 6= µl.
Satz 2.32 (Streuungszerlegung). Im ANOVA1-Modell denieren wir das i-te Gruppenmittel, i =1, . . . , k, bzw. das Gesamtmittel als
Y i• :=1
ni
ni∑j=1
Yij bzw. Y •• :=1
n
k∑i=1
ni∑j=1
Yij
sowie
SSB :=
k∑i=1
ni(Y i• − Y ••)2 und SSW :=
k∑i=1
ni∑j=1
(Yij − Y i•)2
(SSB: sum of squares between groups; SSW: sum of squares within groups ). Dann gilt
SST :=
k∑i=1
ni∑j=1
(Yij − Y ••)2 = SSB + SSW.
Beweis. Es gilt
SST =∑i
∑j
(Yij − Y ••)2 =∑i
∑j
(Yij − Y i• + Y i• − Y ••)2
=∑i
∑j
((Yij − Y i•)2 + 2(Yij − Y i•)(Y i• − Y ••) + (Y i• − Y ••)2
),
wobei ∑i
∑j
(Yij − Y i•)(Y i• − Y ••) =∑i
(Y i• − Y ••)∑j
(Yij − Y i•)
=∑i
(Y i• − Y ••)(niY i• − niY i•) = 0.
Oenbar spricht es gegen die Nullhypothese, wenn die Streuung zwischen den Gruppen gröÿerist als die Streuung innerhalb der Gruppen. Dies motiviert sowohl den Namen ANOVA als auchfolgende Methode:
27
Fg Quadratsummen Quadratmittel F-Statistikzwischen k − 1 SSB =
∑ki=1 ni(Y i• − Y ••)2 SSB/(k − 1)
n− kk − 1
SSB
SSWinnerhalb n− k SSW =
∑ki=1
∑nij=1(Yij − Y i•)2 SSW/(n− k)
total n− 1 SST =∑ki=1
∑nij=1(Yij − Y ••)2 SST/(n− 1)
Tabelle 1: ANOVA-Tafel
Methode 8: Einfaktorielle Varianzanalyse (ANOVA1). Im Modell der einfaktoriellen Va-rianzanalyse testen wir
H0 : µ1 = · · · = µk versus H1 : ∃i, l ∈ 1, . . . , k : µi 6= µl
zum Niveau α ∈ (0, 1) durch den F-Test
ϕµ(Y ) = 1F (Y )>qF (k−1,n−k);1−α mit F (Y ) :=n− kk − 1
SSB
SSW,
wobei qF (k−1,n−k);1−α das (1− α)-Quantil der F (k − 1, n− k)-Verteilung ist.
Satz 2.33. Im einfaktoriellen Varianzanalysemodell gilt:
(i) Der Kleinste-Quadrate-Schätzer von µ = (µ1, . . . , µk)> ist gegeben durch µ =(Y 1•, . . . , Y k,•)
>.
(ii) SSW/σ2 ∼ χ2(n− k) und unter H0 gilt SSB/σ2 ∼ χ2(k − 1)
(iii) SSW und SSB sind unabhängig und somit F := n−kk−1
SSBSSW
H0∼ F (k − 1, n− k).
Beweis. (i) Nachrechnen zeigt
µ = (X>X)−1X>Y =
1/n1 0. . .
0 1/nk
∑n1
j=1 Y1j
...∑nkj=1 Ykj
=
Y 1•...
Y k•
.
(ii)+(iii) Wegen RSS = |Y −Xµ|2 = SSW folgt SSW/σ2 ∼ χ2(n−k) und die Unabhängigkeitvon SSW und µ aus Lemma 2.22. Nach dem vorangegangen Satz gilt weiterhin SSB = SST −SSW . Somit folgt die Behauptung aus Satz 2.27, falls SST = RSSH0 . Nun gilt
RSSH0= min
µ∈R
∣∣∣Y −Xµ...µ
︸ ︷︷ ︸∈Rk
∣∣∣2 = minµ∈R
∣∣∣Y −1...1
︸ ︷︷ ︸
=:X0∈Rn•×1
µ∣∣∣2.
Dieses Minimierungsproblem wird gelöst durch µH0= (X>0 X0)−1X>0 Y = n−1
∑i,j Yij = Y ••.
Damit folgt RRSH0= SST .
Bemerkung 2.34. In der Eektdarstellung wird das einfaktorielle Varianzanalysemodell als
Yij = µ0 + αi + εij , i = 1, . . . , k, j = 1, . . . , ni,
geschrieben mit Intercept µ0 := 1n
∑ki=1 niµi = E[Y••] und αi := µi − µ0, den Eekt der
Faktorstufe i = 1, . . . , k. Insbesondere muss in dieser Darstellung die Nebenbedingung 0 =∑ki=1 niαi oder äquivalent nkαk = −
∑k−1i=1 niαi beachtet werden, damit die Designmatrix weiter
vollen Rang hat. Der Parametervektor ist also gegeben durch (µ0, α1, . . . , αk−1)>. Die F-Statistikum die Globalhypothese H0 : α1 = · · · = αk−1 = 0 zu überprüfen, ist identisch zur Statistik ausSatz 2.33.
28
Beispiel 2.35 (Zweistichproben t-Test). Soll die Gleichwertigkeit von bspw. zwei Düngemittelngetestet werden, ist k = 2 und das Testproblem H0 : µ1 = µ2 versus H1 : µ1 6= µ2. WegennY •• = n1Y 1• + n2Y 2• gilt
SSB = n1(Y 1• − Y ••)2 + n2(Y 2• − Y ••)2
= n1Y2
1• + n2Y2
2• + nY2
•• − 2(n1Y 1• + n2Y 2•)Y ••
= n1Y2
1• + n2Y2
2• −1
n
(n1Y
2
1• + n2Y2
2•)2
=n1n2
n
(Y 1• − Y 2•
)2.
Somit ist
ϕ = 1|T |>qt(n−2),1−α/2 mit T :=Y 1• − Y 2•√
( 1n1
+ 1n2
)SSW/(n− 2)
mit dem (1− α/2)-Quantil der t(n− 2)-Verteilung qt(n−2),1−α/2 ein Test der Hypothese H0 zumNiveau α ∈ (0, 1).
Denition 2.36. Das Modell der zweifaktoriellen Varianzanalyse mit balanciertem Design (ANO-VA2) ist gegeben durch Beobachtungen
Yijk =µij + εijk, i = 1, . . . , I, j = 1, . . . , J, k = 1, . . . ,K
=µ0 + αi + βj + γij + εijk
mit I, J,K > 2, iid.-verteilten Störgröÿen εijk ∼ N (0, σ2) und Nebenbedingungen (der Eektdar-stellung)
I∑i=1
αi =
J∑j=1
βj =
I∑i=1
γij =
J∑j=1
γij = 0.
Wir haben also zwei Faktoren mit Faktorstufen i = 1, . . . , I und j = 1, . . . , J . (αi) bzw.(βj) heiÿen Haupteekte des ersten bzw. zweiten Faktors. (γij) heiÿen Interaktions- bzw.Wechselwirkungseekte.
Das ANOVA2-Modell ist also ein lineares Modell mit zwei kategoriellen Kovariablen. Die Ge-samtanzahl an Beobachtungen ist gegeben durch n = I · J ·K. Die typische Testprobleme sind
H0 : ∀i : αi = 0 versus H1 : ∃i ∈ 1, . . . , I : αi 6= 0, (2.2)
H0 : ∀j : βj = 0 versus H1 : ∃j ∈ 1, . . . , J : βj 6= 0, (2.3)
H0 : ∀i, j : γij = 0 versus H1 : ∃i ∈ 1, . . . , I, j ∈ 1, . . . , J : γij 6= 0. (2.4)
Satz 2.37. Im zweifaktoriellen Varianzanalysemodell mit balanciertem Design gilt:
(i) Die Kleinsten-Quadrate-Schätzer für µ0, αi, βj und γij, i = 1, . . . , I − 1, j = 1, . . . , J − 1,sind gegeben durch (• heiÿt, dass über die jeweilige Koordinate gemittelt wird)
µ0 = Y •••, αi = Y i•• − Y •••, βj = Y •j• − Y •••,
γij = (Y ij• − Y •••)− αi − βj = Y ij• − Y i•• − Y •j• + Y •••.
(ii) Denieren wir
SSW :=
I∑i=1
J∑j=1
K∑k=1
(Yijk − Y ij•)2,
SSB1 := JK
I∑i=1
(Y i•• − Y •••)2, SSB2 := IK
J∑j=1
(Y •j• − Y •••)2,
SSB12 := K
I∑i=1
J∑j=1
(Y ij• − Y i•• − Y •j• + Y •••)2,
29
dann können die Hypothesen (2.2), (2.3) bzw. (2.4) mit den F-Statistiken
IJ(K − 1)
I − 1
SSB1
SSW∼ F (I − 1, IJ(K − 1)),
IJ(K − 1)
J − 1
SSB2
SSW∼ F (J − 1, IJ(K − 1)) bzw.
IJ(K − 1)
(I − 1)(J − 1)
SSB12
SSW∼ F
((I − 1)(J − 1), IJ(K − 1)
)getestet werden.
Beweis. Übung 2.
Bemerkung 2.38. Selbstverständlich erhält man analoge Resultate, wenn wir für jede Zelle (i, j) ∈1, . . . , I × 1, . . . , J verschiedene Stichprobenumfänge nij > 2 beobachten.
Beispiel 2.39. Ein Bauer möchte wissen ob die Gröÿe seiner geernteten Kohlköpf sich für zweiverschiedene Kultursorten unterscheidet. Auch der Panztag könnte eine Rolle spielen.
3 Exponentialfamilien and verallgemeinerte lineare Modelle
3.1 Die Informationsungleichung
Der Satz von Gauÿ-Markov hat uns bereits ein Optimalitätsresultat geliefert, dass allerdings auflineare Schätzer im linearen Modell eingeschränkt ist. Wir suchen nun allgemeiner nach unverzerr-ten Schätzern deren Schätzwerte möglichst wenig um den korrekten Wert streuen.
Denition 3.1. Sei (X ,F , (Pϑ)ϑ∈Θ) ein statistisches Modell. Ein erwartungstreuer Schätzer Teines abgeleiteten Parameters ρ(ϑ) heiÿt varianzminimierend bzw. (gleichmäÿig) bester Schätzer(UMVUE: uniformly minimum variance unbiased estimator), wenn für jeden weiteren erwartungs-treuen Schätzer S gilt:
Varϑ(T ) 6 Varϑ(S) für alle ϑ ∈ Θ.
Wir werden zunächst eine untere Schranke für die Varianz beweisen und anschlieÿend unter-suchen, für welche Schätzer diese erreicht wird.
Denition 3.2. Ein vom Maÿ µ dominiertes, statistisches Modell (X ,F , (Pϑ)ϑ∈Θ) heiÿt regulär,wenn die folgenden Eigenschaften erfüllt sind:
(i) Θ ist eine oene Menge in Rd, d > 1.
(ii) Die Likelihood-Funktion L(ϑ, x) ist auf Θ×X strikt positiv und nach ϑ stetig dierenzierbar.Bezeichnen wir den Gradienten in ϑ mit ∇ϑ = ( ∂
∂ϑ1, . . . , ∂
∂ϑd)>, existiert insbesondere die
Scorefunktion
Uϑ(x) := ∇ϑ logL(ϑ, x) =∇ϑL(ϑ, x)
L(ϑ, x).
(iii) Für jedes ϑ ∈ Θ existiert die Fisher-Information
I(ϑ) := Eϑ[Uϑ(X)Uϑ(X)>
]und ist positiv denit.
(iv) Es gilt die Vertauschungsrelationˆh(x)∇ϑL(ϑ, x)µ(dx) = ∇ϑ
ˆh(x)L(ϑ, x)µ(dx) (3.1)
für h(x) = 1.
30
Ein Schätzer T : X → R heiÿt regulär, falls E[|T (X)|2] <∞ und (3.1) auch für h(x) = T (x) gilt.
Bemerkung 3.3.
(i) Der Satz von Lebesgue liefert eine hinreichende Bedingung für die Vertauschungsrelation(3.1): Sie gilt falls für jedes ϑ0 ∈ Θ eine Umgebung Vϑ0 ⊆ Θ existiert, so dass
ˆX
supϑ∈Vϑ0
∣∣∣∇ϑL(ϑ, x)∣∣∣µ(dx) <∞.
Auÿerdem kann man (3.1) für jedes gegebene Modell (und jeden Schätzer) explizit nachprü-fen.
(ii) Als Konsequenz von (3.1) ergibt sich
Eϑ[Uϑ] =
ˆ∇ϑL(ϑ, x)µ(dx) = ∇ϑ
ˆL(ϑ, x)µ(dx) = ∇ϑ1 = 0
und damit Varϑ(Uϑ) = I(ϑ).
(iii) Ist L(ϑ, x) in ϑ zweimal stetig dierenzierbar und gilt (3.1) mit h(x) = 1 und L ersetzt mit∂L∂ϑi
für alle i ∈ 1, . . . , d, dann gilt I(ϑ) = −Eϑ[HUϑ(X)(ϑ)] für die Hesse-Matrix HU·(x) der
Scorefunktion ϑ 7→ Uϑ(x) (Übung 2).
(iv) Warum heiÿt I(ϑ) Information? Erstens: I(ϑ) = 0 gilt auf einer Umgebung Θ0 ⊆ Θ genaudann, wenn Uϑ(x) = 0 für alle ϑ ∈ Θ0 und µ-f.a. x ∈ X , also wenn L(ϑ, x) für µ-f.s. konstantist und somit keine Beobachtung die Parameter in Θ0 unterscheiden kann (dieser Fall istdaher in der Denition ausgeschlossen). Zweitens, verhält sich die Fisher-Information beiunabhängigen Beobachtungen additiv: Ist (X ,F , (Pϑ)ϑ∈Θ) ein reguläres Modell mit Fisher-Information I, so hat das Produktmodell (Xn,F⊗n, (P⊗nϑ )ϑ∈Θ) die Fisher-Information I⊗n =nI (Beweis als Übung 2).
Satz 3.4 (Cramér-Rao-Ungleichung, Informationsschranke). Gegeben seien ein reguläres statisti-sches Modell (X ,F , (Pϑ)ϑ∈Θ), eine zu schätzende stetig dierenzierbare Funktion ρ : Θ→ R, undein regulärer erwartungstreuer Schätzer T von ρ. Dann gilt
Varϑ(T ) >(∇ρ(ϑ)
)>I(ϑ)−1∇ρ(ϑ) für alle ϑ ∈ Θ. (3.2)
Beweis. Aus der Zentriertheit von Uϑ und der Regularität und Erwartungstreue von T erhaltenwir
Covϑ(Uϑ, T ) = Eϑ[TUϑ] =
ˆXT (x)∇ϑL(ϑ, x)µ(dx)
= ∇ˆXT (x)L(ϑ, x)µ(dx) = ∇Eϑ[T ] = ∇ρ
für alle ϑ ∈ Θ. Für jeden Vektor e ∈ Rd ergibt die Cauchy-Schwarz-Ungleichung somit
(〈e,∇ρ〉)2 = Covϑ(〈e, Uϑ〉, T )2 6 Varϑ(〈e, Uϑ〉) Varϑ(T ) = 〈I(ϑ)e, e〉Varϑ(T ),
also
Varϑ(T ) >(〈∇ρ, e〉)2
〈I(ϑ)e, e〉.
Maximieren über e ∈ Rd ergibt mit e = I(ϑ)−1∇ρ(ϑ) die Behauptung.
Denition 3.5. Ein regulärer erwartungstreuer Schätzer für den Gleichheit in (3.2) gilt, heiÿtCramér-Rao-ezient.
31
Im Folgenden beschränken wir uns auf einparametrische (d = 1) Modelle.
Satz 3.6. Unter den Bedingungen von Satz 3.4 mit Θ ⊆ R erreicht der Schätzer T die untereSchranke für alle ϑ ∈ Θ genau dann, wenn µ-f.ü. gilt
T − ρ(ϑ) = ρ′(ϑ)I(ϑ)−1Uϑ für alle ϑ ∈ Θ.
Falls ρ′ 6= 0 ist dies äquivalent zu
L(ϑ, x) = exp(η(ϑ)T (x)− ζ(ϑ)
)c(x),
wobei η : Θ → R eine Stammfunktion von I/ρ′, c : X → (0,∞) messbar und ζ(ϑ) =log´c(x) exp(η(ϑ)T (x))µ(dx) eine Normierungsfunktion sind.
Beweis. Denieren wir v(ϑ) := ρ′(ϑ)I−1(ϑ) (konstant in x) erhalten wir wegen Covϑ(Uϑ, T ) =ρ′(ϑ)
0 6Varϑ(T − v(ϑ)Uϑ
)= Varϑ(T ) + v(ϑ)2 Varϑ(Uϑ)− 2v(ϑ) Covϑ(Uϑ, T ) = Varϑ(T )− ρ′(ϑ)2I−1(ϑ),
also wieder die Informationsungleichung. Gleichheit gilt genau dann, wenn T − v(ϑ)Uϑ Pϑ-f.s.konstant also gleich seinem Erwartungswert ρ(ϑ) ist. Da Pϑ eine strikt positive µ-Dichte hat giltµ(T − ρ(ϑ) 6= v(ϑ)Uϑ) = 0. Wenn dies nun für alle ϑ ∈ Θ gilt, so folgt sogar
µ(T − ρ(ϑ) 6= v(ϑ)Uϑ für ein ϑ ∈ Θ
)= 0,
denn aus Stetigkeitsgründen kann man sich auf rationale ϑ beschränken und die abzählbare Ver-einigung von Nullmengen ist wieder eine Nullmenge. Die explizite Form der Likelihood-Funktionfolgt durch unbestimmte Integration bzgl. ϑ.
Dieser Satz führt uns in natürlicher Weise auf eine wichtige Klasse von statistischen Modellen:
Denition 3.7. Es sei (X ,F , (Pϑ)ϑ∈Θ) ein von µ dominiertes statistisches Modell mit Θ ⊆ Roen. Dann heiÿt (Pϑ)ϑ∈Θ (einparametrische) Exponentialfamilie in η(ϑ) und T , wenn messbareFunktionen η : Θ→ R, T : X → R und c : X → (0,∞) existieren, so dass
dPϑdµ
(x) = c(x) exp(η(ϑ)T (x)− ζ(ϑ)
), x ∈ X , ϑ ∈ Θ,
wobei ζ(ϑ) := log´c(x) exp(η(ϑ)T (x))µ(dx). Dabei wird angenommen, dass T nicht µ-f.s. kon-
stant ist. η(ϑ) heiÿt natürlicher Parameter der Exponentialfamilie und
Ξ :=η ∈ R :
ˆXc(x)eηT (x)µ(dx) ∈ (0,∞)
heiÿt natürlicher Parameterraum. Ist die Exponentialfamilie durch η ∈ Ξ parametrisiert, dannwird sie als natürliche Exponentialfamilie bezeichnet.
Bemerkung 3.8.
(i) Die Darstellung ist nicht eindeutig, mit a 6= 0 erhält man beispielsweise eine Exponen-tialfamilie in η(ϑ) = aη(ϑ) und T (x) = T (x)/a. Auÿerdem kann die Funktion c in dasdominierenden Maÿ absorbiert werden: µ(dx) := c(x)µ(dx).
(ii) Die Identizierbarkeitsfordung Pϑ 6= Pϑ′ für alle ϑ 6= ϑ′ ist äquivalent zur Injektivität von η.
Beispiel 3.9.
32
(i) (N (µ, σ2))µ∈R mit σ > 0 bekannt ist eine Exponentialfamilie in η(µ) = µ/σ2 und T (x) = x:
L(ϑ, x) =1√
2πσ2e−(x−µ)2/(2σ2) =
1√2πσ2
e−(x2−2µx+µ2)/(2σ2), x ∈ R.
(ii) Die Familie der Poissonverteilungen (Poiss(λ))λ>0 mit Intensitätsparameter λ bildet eineExponentialfamilie mit natürlichem Parameter η(λ) = log λ und T (x) = x:
L(λ, x) =λxe−λ
x!=
1
x!ex log λ−λ, x ∈ Z+.
Lemma 3.10. Ist ein statistisches Modell durch eine Exponentialfamilie in η : Θ→ R und T : X →R mit dierenzierbarem η gegeben, so ist dieses regulär. Ferner gilt
(i) Jede Statistik S : X → R mit existierendem Erwartungswert ist regulär. ρ(ϑ) := Eϑ[T ] iststetig dierenzierbar mit ρ′(ϑ) = η′(ϑ) Varϑ(T ) 6= 0, ϑ ∈ Θ.
(ii) Die Normierungsfunktion ζ ist auf Θ ⊆ R stetig dierenzierbar mit ζ ′(ϑ) = η′(ϑ)Eϑ[T ] fürϑ ∈ Θ. Die Scorefunktion ist Uϑ = η′(ϑ)T − ζ ′(ϑ).
(iii) Für die Fisher-Information gilt I(ϑ) = η′(ϑ)ζ ′′(ϑ) = η′(ϑ)ρ′(ϑ) für alle ϑ ∈ Θ.
Beweis. O.B.d.A. ist η(ϑ) = ϑ und somit η′ = 1 für alle ϑ ∈ Θ. Der allgemeine Fall ergibt sichdurch Reparametrisierung und Anwendung der Kettenregel.
Schritt 1: Sei S eine beliebige reelle Statistik mit S ∈ L1(Pϑ) für alle ϑ ∈ Θ. Dann ist dieFunktion
uS(ϑ) := eζ(ϑ)Eϑ[S] =
ˆXS(x)eϑT (x)c(x)µ(dx)
auf Θ wohl deniert. Wir zeigen nun, dass us beliebig oft dierenzierbar ist.Ist ϑ ∈ Θ und t ∈ R so klein, dass auch ϑ± t ∈ Θ, so gilt mittels monotoner Konvergenz∑
k>0
|t|k
k!
ˆX|S(x)||T (x)|keϑT (x)c(x)µ(dx) =
ˆX|S(x)|eϑT (x)+|tT (x)|c(x)dx
6ˆX|S(x)|
(e(ϑ+t)T (x) + e(ϑ−t)T (x)
)c(x)dx <∞.
Also ist ST k ∈ L1(Pϑ) für alle ϑ ∈ Θ und insbesondere T ∈ L2(Pϑ) für alle ϑ. Ferner ist die Reihe∑k>0
tk
k!
ˆXS(x)T (x)keϑT (x)c(x)µ(dx)
absolut konvergent und Summation und Integration können vertauscht werden. Die Reihe nimmtdaher den Wert uS(ϑ+ t) an. Damit ist uS sogar analytisch.
Schritt 2: Es folgt u′S(ϑ) = eζ(ϑ)Eϑ[ST ] und insbesondere u′1(ϑ) = u1(ϑ)Eϑ[T ] sowie u′′1(ϑ) =u1(ϑ)Eϑ[T 2]. Für ζ(ϑ) = log u1(ϑ) bekommen wir also ζ ′(ϑ) = Eϑ[T ] =: ρ(ϑ) und
ρ′(ϑ) = ζ ′′(ϑ) = u′′1(ϑ)/u1(ϑ)− (u′1(ϑ)/u1(ϑ))2 = Varϑ(T ).
Aus der Dierenzierbarkeit von ζ folgt
Uϑ =∂
∂ϑlogL(ϑ, x) = T − ζ ′(ϑ), ϑ ∈ Θ
und somit I(ϑ) = Varϑ(Uϑ) = Varϑ(T ) > 0. Weiter können wir schreiben
d
dϑEϑ[S] =
(uS(ϑ)e−ζ(ϑ)
)′=(u′S(ϑ)− uS(ϑ)ζ ′(ϑ)
)e−ζ(ϑ)
=Eϑ[ST ]− Eϑ[S]ζ ′(ϑ) = Eϑ[SUϑ]
=
ˆXS(x)
∂
∂ϑL(ϑ, x)µ(dx).
Daher gilt einerseits (3.1) für alle h ∈ L1(Pϑ) und anderseits folgt die Regularität des Modells.
33
Korollar 3.11 (Existenz von besten Schätzern). Für jedes statistische Modell gegeben durch eineExponentialfamlie mit dierenzierbarem η und η′ 6= 0 ist die zugrunde liegende Statistik T einbester und Cramér-Rao-ezienter Schätzer für ρ(ϑ) := Eϑ[T ] = ζ ′(ϑ)/η′(ϑ). In dem Fall gilt
Varϑ(T ) = ρ′(ϑ)/η′(ϑ) und I(ϑ) = η′(ϑ)ρ′(ϑ) für alle ϑ ∈ Θ.
Für natürliche Exponentialfamilien gilt also insbesondere Varη(T ) = I(η).
Beweis. Folgt unmittelbar aus Satz 3.4 und Lemma 3.10. Für natürliche Exponentialfamiliengilt also Varη(T ) = ρ′(η) = I(η) und die Informationsschranke ist gegeben durch ρ′(η)2/I(η) =I(η).
Beispiel 3.9 (fortgesetzt).
(i) (N (µ, σ2))µ∈R und bekanntem σ > 0 ist wie oben gesehen eine Exponentialfamilie in η(µ) =µ/σ2, T (x) = x und mit ζ(µ) = µ2/(2σ2). Somit ist ρ(µ) = Eµ[T ] = µ und Varµ(T ) = σ2.Da T nicht von σ > 0 abhängt, ist T sogar bester Schätzer für den Erwartungswert für alleNormalverteilungen.
(ii) Für die Exponentialfamilie (Poiss(λ))λ>0 in η(λ) = log λ und T (x) = x gilt ζ(λ) = λ. Wegenρ(λ) = Eλ[T ] = λ und Varλ(T ) = λ ist T bester Schätzer für λ.
Lemma 3.12. Ist (Pϑ)ϑ∈Θ auf (X ,F ) eine Exponentialfamilie in η : Θ → R und T : X → Rso ist (P⊗nϑ )ϑ∈Θ eine Exponentialfamilie auf (Xn,F⊗n) mit zugrundeliegender Statistik Tn =1n
∑ni=1 T Xi. Ist η dierenzierbar mit η′ 6= 0, folgt insbesondere, dass Tn ein bester Schätzer für
ρ(ϑ) = Eϑ[T ] ist.
Beweis. Übung 2.
Abschlieÿend klären wir noch die Frage was das Maximum-Likelihood-Prinzip für natürlicheExponentialfamilien ergibt.
Lemma 3.13. Ist (Pϑ)ϑ∈Θ auf (X ,F ) eine natürliche Exponentialfamilie in η ∈ Ξ und T : X →R, dann ist T auf dem Ereignis T (X) ∈ ran(ζ ′) der eindeutige Maximum-Likelihood-Schätzerdes Parameters ρ(η) := Eη[T ]. Ferner ist ζ ′ : Θ → R invertierbar und der eindeutige Maximum-Likelihood-Schätzer des natürlichen Parameters η ist gegeben durch
η = (ζ ′)−1(T (X)).
Beweis. Um die Maximalstelle der Likelihood-Funktion zu nden, setzen wir die Scorefunktiongleich null. Auf T (X) ∈ ran(ζ ′) gilt
∂η logL(η, x) = Uη(x) = 0 ⇔ T (x) = ζ ′(η).
Da ∂2η logL(η, x) = −ζ ′′(η) = −Varη(T ) < 0, ist η 7→ − logL(η, x) konvex und somit T der eindeu-
tige Maximum-Likelihood-Schätzer des Parameters ρ(η) = ζ ′(η). Aus ζ ′′ > 0 folgt auÿerdem, dassζ ′ invertierbar ist, so dass der Maximum-Likelihood-Schätzer des natürlichen Parameters gegebenist durch (ζ ′)−1 T .
3.2 Verallgemeinerte Lineare Modelle
Mit Hilfe von Exponentialfamilien wollen wir nun lineare Modelle verallgemeinern. Wie in Bei-spiel 3.9 gesehen bildet (N (µ, σ2))µ∈R eine Exponentialfamilie mit natürlichem Parameter η(µ) =µ/σ2 und Statistik T (x) = x, die ein ezienter Schätzer des Parameters ρ(µ) = Eµ[T ] = µ ist. Imgewöhnlichen linearen Modell sind nun die Beobachtungen gegeben durch
Rn 3 Y = Xβ + ε,
34
mit Parametervektor β ∈ Rk, Designmatrix X ∈ Rn×k und iid. Fehlervariablen εi ∼ N (0, σ2) mitVarianz σ > 0. Schreiben wir die Designmatrix als
X =
x1
...xn
mit Zeilenvektoren x1, . . . , xn ∈ Rk,
ist Beobachtung Yi gemäÿ N (xiβ, σ2) verteilt, folgt also einer Exponentialfamilie mit ηi(β) =
xiβ/σ2 und ρi(β) = xiβ, i = 1, . . . , n. Lassen wir nun andere Exponentialfamilien zu, können wir
sowohl Situationen modellieren in den der Zusammenhang zwischen E[Yi] und den Kovariablen(codiert in der Designmatrix X) nichtlinear ist als auch diskrete Beobachtungen Yi zulassen.
Denition 3.14. Auf einem Produktmodell (Xn,F⊗n) liegt ein verallgemeinertes lineares Modell(GLM: generalized linear model) mit n unabhängigen Beobachtungen Y1, . . . , Yn vor, falls dieRandverteilungen von Yi durch natürliche Exponentialfamilien gegeben sind mit Dichten
dPYiηidµ
(yi) = exp(ηiyi − ζ(ηi)
ϕ
)c(yi, ϕ), i = 1, . . . , n,
bzgl. einem dominierenden Maÿ µ, mit unbekanntem Dispersionsparameter ϕ > 0,
ηi ∈ Ξ =η ∈ R :
ˆXeηy/ϕc(y, ϕ)µ(dy) ∈ (0,∞)
⊆ R
für alle i und bekannten Funktionen ζ : Ξ → R und c : X → R+ mit ζ ′′(η) > 0 für alle innerenPunkte η ∈ Ξ. Setze ρ(ηi) := Eη[Yi]. Für einen unbekannten Parametervektor β ∈ Rk, eineDesignmatrix X ∈ Rn×k und eine bijektive, stetig dierenzierbare Funktion g : R→ R gelte weiterg(ρ(η1))
...g(ρ(ηn))
= Xβ.
g heiÿt Linkfunktion. Falls ρ = g−1, gilt (η1, . . . , ηn)> = Xβ und g heiÿt kanonische Linkfunktion(oder kanonischer Link).
Während β der interessierende Parameter ist, wird ϕ als Störparameter angesehen. Für -xiertes ϕ ist Yi also gemäÿ einer natürlichen Exponentialfamilie in T (y) = y/ϕ verteilt. Aus denEigenschaften natürlicher Exponentialfamilien folgt
Eβ,ϕ[Yi] = ζ ′(ηi) und Varβ,ϕ(Yi) = ϕζ ′′(ηi), i = 1, . . . , n.
Beispiel 3.15. Das gewöhnliche lineare Modell ist ein GLM mit kanonischer Linkfunktiong(x) = x, ζ(η) = η2/2 und Dispersionsparameter ϕ = σ2. Lassen wir allgemeinere Linkfunk-tionen zu erhalten wir nicht-lineare Regressionsmodelle (mit normalverteilten Fehlern) gegebendurch Beobachtungen Yi ∼ N (g−1((Xβ)i), ϕ).
Der Dispersionsparameter wird dazu verwendet eine Unterschätzung der (empirisch beobach-teten) Varianz durch das Modell auszugleichen (siehe Übung 2).
Um den unbekannten Parametervektor β in einem verallgemeinerten linearen Modell zu schät-zen, verwenden wir den Maximum-Likelihood-Ansatz. Da ζ ′ streng monoton wachsend und dieLinkfunktion g invertierbar sind, existiert die Funktion ψ := (g ρ)−1. Ist xi ∈ Rk wieder die i-teZeile von X, kann Loglikelihood-Funktion geschrieben werden als
logL(β, ϕ; y) =
n∑i=1
(ψ(xiβ)yi − ζ(ψ(xiβ))
ϕ+ log(c(yi, ϕ))
).
Als notwendige Bedingung an einen Maximum-Likelihood-Schätzer β erhalten wir durch Ableiten
∇β logL(β, ϕ; y) =1
ϕ
n∑i=1
(yi − ρ(ψ(xiβ))
)ψ′(xiβ)x>i = 0. (3.3)
35
Lemma 3.16. In einem verallgemeinerten linearen Modell mit kanonischer Linkfunktion ist dieFisher-Information gegeben durch
I(β) =1
ϕ
n∑i=1
ζ ′′(xiβ)x>i xi ∈ Rk×k.
Ist I(β) positiv denit für alle β und existiert eine Lösung β von (3.3), so ist β der eindeutigeMaximum-Likelihood-Schätzer von β.
Beweis. Aus Lemma 3.10 folgt, dass die Fisher-Information im natürlichen Parameter(η1, . . . , ηn)> gegeben ist durch 1
ϕ
∑ni=1 ζ
′′(ηi). Die Reparametrisierung ηi = xiβ zusammen mitder Kettenregel ergibt die Darstellung von I(β).
Der kanonische Link ist gegeben durch g = ρ−1, so dass ψ in (3.3) die Identität ist. Wegenρ = ζ ′, gilt also
∂2 logL(β, ϕ; y)
∂β∂β>= − 1
ϕ
n∑i=1
ζ ′′(xiβ)x>i xi = −I(β).
Da I(β) > 0, ist β 7→ − logL(β, ϕ; y) streng konvex und somit β der eindeutige Maximum-Likelihood-Schätzer.
Bemerkung 3.17.
(i) Typischerweise besitzt β keine geschlossene Form mehr und muss durch numerische Verfahrenbestimmt werden. Fishers Scoring-Methode verwendet hierfür das iterative Verfahren
β(t+1) = β(t) + I(β)−1∇β logL(β(t), ϕ; y), t = 0, 1, . . .
(Beachte, dass sich der unbekannte Dispersionsparameter ϕ gerade rauskürzt). Für den ka-nonischen Link ist dieses Verfahren äquivalent zur Newton-Raphson-Methode.
(ii) Ist g nicht der kanonische Link ist eine Lösung von (3.3) nicht notwendigerweise einMaximum-Likelihood-Schätzer.
Zwei wichtige Beispielklassen für verallgemeinerte lineare Modelle sind die Poisson-Regressionund die logistische Regression, die abschlieÿend eingeführt werden.
Die Poisson-Regression modelliert unabhängige Poisson-verteilte Beobachtungen, deren Inten-sitätsparameter von Kovariablen abhängen. Sie eignet sich also für Beobachtungen die Zähldaten-struktur haben. Wir hatten bereits gesehen dass die Familie (Poiss(λ))λ>0 eine Exponentialfamiliein η(λ) = log λ und T (x) = x ist: Bezüglich des Zählmaÿes ist die Likelihood-Funktion gegebendurch
L(λ, x) =λxe−λ
x!=
1
x!ex log λ−λ, x ∈ Z+,
und es gilt ρ(λ) = Eλ[T ] = λ.
Denition 3.18. Ein verallgemeinertes lineares Modell auf (Zn+,P(Zn+)) heiÿtPoisson-Regression, falls die unabhängigen Beobachtungen Yi Poiss(λi)-verteilt sind, wobeiλi = eηi mit natürlichen Parametern ηi ∈ R, i = 1, . . . , n, und
η = (η1, . . . , ηn)> = Xβ
mit unbekanntem β ∈ Rk und Designmatrix X ∈ Rn×k.
Bemerkung 3.19. Wir verwenden hier also den kanonischen Link g(λ) = log λ. In der Praxis wirdoft das erweiterte Modell Yi ∼ Poiss(λi · si) verwendet für einen so genannten Zählrahmen si > 0,i = 1, . . . , n. Dann gilt Eβ [Yi] = exp(xiβ + log(si)) mit den Zeilen xi von X. Der Term log(si)wird als Oset bezeichnet, da er jeder Beobachtung einen individuellen Intercept zuweist.
36
Beispiel 3.20. In einem groÿen Krankenhaus wird die Anzahl der Beschwerden über n = 44Notfallärzten untersucht (Daten aus Le (2003)). Der Zählrahmen pro Arzt ist die Anzahl anPatientenbesuchen, die vier zu berücksichtigenden Kovariablen sind Vergütung (in $/h), Erfahrung(in h), Geschlecht und Facharztausbildung (ja/nein).
Lemma 3.21. Die Familie der Bernoulli-Verteilungen (Bernoulli(p))p∈(0,1) bildet eine Exponen-tialfamilie in η(p) = log p
1−p und T (x) = x.
Beweis. Die vom Zählmaÿ dominierte Familie besitzt die Likelihood-Funktion
L(p, x) = px(1− p)1−x = (1− p)( p
1− p
)x= exp
(x log
( p
1− p
)+ log(1− p)
), x ∈ 0, 1.
Denition 3.22. Ein verallgemeinertes lineares Modell auf (0, 1n,P(0, 1n)) heiÿt logistischeRegression, falls die unabhängigen Beobachtungen Yi Bernoulli(pi)-verteilt sind, i = 1, . . . , n, mitnatürlichem Parameterraum R, der kanonischen Link-Funktion g : (0, 1)→ R, g(p) = log(p/(1−p))und
η = (g(p1), . . . , g(pn))> = Xβ
mit unbekanntem β ∈ Rk und Designmatrix X ∈ Rn×k. Die Funktion g heiÿt Logit-Funktion undihre Umkehrfunktion g−1 : R→ (0, 1), g−1(x) = (1 + e−x)−1 heiÿt logistische Funktion.
Bemerkung 3.23. Es gilt also E[Yi] = g−1(ηi) = eηi/(1 + eηi), wobei die Funktion g−1 gerade dieVerteilungsfunktion der standardisierten logistischen Verteilung ist (welche im Allgmeinen einenMittelwerts- und einen Streuungsparameter besitzt). Das motiviert ein populäres Beispiel für einenicht kanonische Linkfunktion: die Probit-Funktion g(λ) = Φ−1(λ) mit der Verteilungsfunktionder Standardnormalverteilung Φ.
Da wir hier ein Modell gefunden haben um 0, 1-wertige Zufallsvariablen durch Kovariablen zuerklären, werden wir die logistische Regression im nächsten Kapitel zur Klassikation verwenden.
3.3 Ergänzung: Numerische Bestimmung des Maximum-Likelihood-
Schätzers
Das vermutlich grundlegendste numerische Verfahren zur Bestimmung von Nullstellen ist dasNewton-Verfahren oder Newton-Raphson-Verfahren:
Ziel: Finde x∗ ∈ R: f(x∗) = 0 für eine Funktion f : R→ R.
Verfahren:
(i) Wähle einen Startpunkt x0 ∈ R (der möglichst nahe an x∗ liegen sollte).
(ii) Approximiere x∗ mit der rekursiven Vorschrift
xn+1 := xn −f(xn)
f ′(xn)falls f ′(xn) 6= 0
Abbruchkriterien: |f(xn)| < ε oder |xn+1 − nn| < ε für ein ε > 0.
Geometrisch ist xn+1 genau die Nullstelle der Tangente y = f(xn)+f ′(xn)(x−xn) an f im Punkt(xn, f(xn)). Im allgemeineren Fall f : Rk → Rk erhalten wir die Rekursionsvorschrift
Jf (xn)(xn+1 − xn) = −f(xn) ⇐⇒ xn+1 = xn − Jf (xn)−1f(xn)
mit der Jacobi-Matrix Jf (x) = ( ∂fi∂xj)i,j=1,...,k ∈ Rk×k falls diese positiv denit ist.
Das Newton-Verfahren soll nun verwendet werden um den Maximum-Likelihood-Schätzer βin einem verallgemeinerten linearen Modell (Xn,F⊗n, (P⊗nβ,ϕ)β∈Rk,ϕ>0) mit kanonischem Link zubestimmen. Zur Erinnerung ist in diesem Fall die Likelihood-Funktion gegeben durch
L(β, ϕ; y) =
n∏i=1
exp( (xiβ)yi − ζ(xiβ)
ϕ
)c(yi, ϕ)
37
mit n Zeilenvektoren xi ∈ Rk. Setzen wir also
f(β) = ∇β logL(β, ϕ; y) =1
ϕ
n∑i=1
(yi − ζ ′(xiβ)
)x>i ,
dann ist die Jacobi-Matrix gleich der Hesse-Matrix der Loglikelihood-Funktion
Jf (β) =(∂βl logL(β, ϕ; y)
∂βj
)l,j=1,...,k
= − 1
ϕ
( n∑i=1
ζ ′′(xiβ)xi,lxi,j
)l,j=1,...,k
= − 1
ϕ
n∑i=1
ζ ′′(xiβ)x>i xi.
Da diese nicht mehr von y abhängt, erhalten wir
Jf (β) = Eβ[HlogL(·,ϕ,Y )(β)
]= −I(β).
Einsetzen in obige Iterationsvorschrift liefert Fishers Scoring-Methode:
β(t+1) = β(t) + I(β)−1∇β logL(β(t), ϕ; y)
= β(t) −( n∑i=1
ζ ′′(xiβ(t))x>i xi
)−1 n∑i=1
(Yi − ζ ′(xiβ(t))
)x>i , t ∈ N,
wobei wir β(0) = 0 setzen.
4 Klassikation
Während im linearen Modell die Zielvariable quantitativ ist, gibt es viele Situationen in denendie Daten qualitativ bzw. kategoriell sind. Das Grundprinzip der Klassikation ist anhand einersogenannten Trainingsmenge (x1, Y1), . . . , (xn, Yn) zu lernen, die Klassen zu unterscheiden, umanschlieÿend vorherzusagen, zu welcher Klasse Beobachtungen zu neuen xn+1, . . . , xn+m gehören(statistisches Lernen). Anders ausgedrückt, sollen xn+1, . . . , xn+m klassiziert werden.
Beispiel 4.1. Auf Grundlage vom monatlichen Kontostand der Kreditkarte und dem Jahresein-kommen soll vorhergesagt werden ob jemand zahlungsunfähig wird oder nicht. Als Trainingsda-tensatz haben wir Daten (xi,1, xi,2, Yi) für i = 1, . . . , n Personen gegeben, wobei xi,1 bzw. xi,2der monatliche Kontostand der Kreditkarte bzw. das Jahreseinkommen von Person i sind und Yidie Frage Zahlungsunfähig? mit Ja oder Nein beantwortet (simulierter default-Datensatzaus James et al. (2013)). Etwa 3% der Personen sind zahlungsunfähig. Beachte, dass in realenAnwendungen die Beziehung zwischen Ko- und Zielvariablen typischerweise nicht so eindeutigsind.
4.1 Logistische Regression
Stammen die Zielvariablen nur aus zwei verschiedenen Klassen (die mit 0 und 1 codiert werden),bietet sich die logistische Regression aus Kapitel 3 als Modell an. Zur Erinnerung heiÿt ein stati-sches Experiment (0, 1n,P(0, 1n), (Pβ⊗n)β∈Rk) multiple logistische Regression mit Kovaria-blen xi = (1, xi,1, . . . , xi,k−1) ∈ Rk (Zeilenvektor mit Absolutglied), falls Yi Bernoulli(pi)-verteiltist, wobei pi = p(xi, β) gegeben ist durch
log( pi
1− pi
)= β0 +
k−1∑j=1
βjxi,j für i = 1, . . . , n.
Äquivalent gilt
p(xi, β) =exiβ
1 + exiβ.
Die Wahrscheinlichkeit, dass Yi zur Klasse 1 gehört, wird also durch die k−1 Kovariablen erklärt.
38
Methode 9: Klassikation mittels logistischer Regression. Nach Schätzung des Parame-tervektors β auf der Trainingsmenge (xi, Yi)i=1,...,n können wir für eine jede neue Kovariablenrea-lisierung xn+1 = (1, xn+1,1, . . . , xn+1,k−1) (Zeilenvektor) einen zugehörigen Wert
pn+1 = p(xn+1, β) =exn+1β
1 + exn+1β
vorhersagen und xn+1 der Klasse 1 zuordnen, falls pn+1 > τ für einen Schwellenwert τ ∈ [0, 1].Andernfalls klassizieren wir xn+1 mit 0.
Der Maximum-Likelihood-Ansatz führt auf die Maximierung der Loglikelihood-Funktion
`(β, y) := logL(β, y) =
n∑i=1
(yi log p(xi, β) + (1− yi)(1− p(xi, β)
)=
n∑i=1
(yi(xiβ)− log
(1 + exiβ
)).
Nullsetzen des Gradienten führt auf k Gleichungen, die nichtlinear in β sind. Um den Maximum-Likelihood-Schätzer numerisch zu bestimmen, verwenden wir wieder das Newton-Verfahren. Diesesführt uns auf die iterativ neugewichteten Kleinste-Quadrate-Methode (IRLS: iteratively reweightedleast squares):
Lemma 4.2. In der logistischen Regression mit Designmatrix X von vollem Rang ist der (t+1)steIterationsschritt von Fishers Scoring-Methode gegeben durch
β(t+1) = (X>Wβ(t)X)−1X>Wβ(t)Zβ(t) = arg minb
∣∣W 1/2
β(t)(Zβ(t) −Xb)
∣∣2mit adjustiertem Responsevektor Zβ = Xβ +W−1
β (Y − pβ), wobei
pβ =(p(x1, β), . . . , p(xn, β)
)> ∈ Rn und
Wβ = diag(p(x1, β)(1− p(x1, β), . . . , p(xn, β)(1− p(xn, β))
)∈ Rn×n.
Beweis. Für die logistische Funktion g(x) = ex/(1 + ex) gilt g′(x) = g(x)(1 − g(x)). Aus derexpliziten Form der Loglikelihood-Funktion `(β) folgt damit, dass Scorefunction und Hesse-Matrixgegeben sind durch
∇β`(β, y) = X>(y − pβ) bzw. H`(·,y)(β) = −X>WβX.
Somit ist der Iterationsschritt von Fishers Scoring-Methode
β(t+1) = β(t) +(X>Wβ(t)X
)−1X>(Y − pβ(t))
=(X>Wβ(t)X
)−1X>Wβ(t)
(Xβ(t) +W−1(Y − pβ(t))
)=(X>Wβ(t)X
)−1X>Wβ(t)Zβ(t) .
Wie in Kapitel 2 gesehen, ist dies gerade die Lösung des gewichteten Kleinste-Quadrate-Problems.
Bemerkung 4.3.
(i) Der Maximum-Likelihood-Schätzer β erfüllt ∇β`(β, y) =∑ni=1 x
>i (y − p(xi, β)) = 0. Da der
erste Koezient von xi gleich 1 ist folgt∑i yi =
∑i p(xi, β), d.h. die erwartete Anzahl der
Beobachtungen in Klasse eins stimmt mit der beobachteten Anzahl überein.
39
(ii) Insbesondere zeigt dieses Lemma, dass der Maximum-Likelihood-Schätzer β die Lösung einesgewichteten Kleinste-Quadrate-Problems mit Responsevektor Zβ = Xβ+W−1
β(Y −pβ) und
Gewichten wi = pi(1 − pi), wobei beides wieder von β abhängt. Die gewichteten Quadrat-summen der Residuen sind dann
n∑i=1
(Yi − pi)2
pi(1− pi)
und messen die Abweichung der Daten von der Modellvorhersage.
Beispiel 4.4. Wir betrachten wieder den Datensatz aus Beispiel 4.1, wobei wir eine zusätzlicheKovariable Student mit Werten Ja oder Nein zur Verfügung haben. Eine Logistische Re-gression, die nur Student und einen Intercept verwendet führt zu einem positiven Koezientender Dummy-Variable, d.h. die Zahlungsunfähigkeitswahrscheinlichkeit ist für Studenten höher alsfür Nicht-Studenten. Verwenden wir alle drei Kovariablen erhalten wir jedoch einen negativenZusammenhang! Wie kann man diesen so genannten Konfundierungseekt erklären?
Die logistische Regression kann auch auf mehr als zwei Klassen ausgeweitet werden, indem wirstatt der Bernoulli-Verteilung die Multinomialverteilung verwenden. Häug wird jeden die Me-thode des nächsten Abschnittes dieser Variante vorgezogen. Inferenz für die logistische Regressionberuht auf asymptotischen Überlegungen auf die wir in dieser Vorlesung nicht weiter eingehenwerden.
4.2 Bayesklassikation1
Die logistische Regression modelliert die Wahrscheinlichkeit P(Y = 1) unter Verwendung des Re-gressorvektors x für zwei Klassen 0 und 1. Verstehen wir den Kovariablenvektor als ZufallsvariableX, wird also die bedingte Wahrscheinlichkeit P(Y = 1|X = x) der Klasse 1 gegeben einer Kova-riablenrealisierung X = x modelliert. Stattdessen wird nun ein Bayesianischer Ansatz verfolgt.
Gegeben sei das zufällige Paar(X,Y
), welches Werte in Rd ×
1, . . . ,K
annimmt. Hierbei
bezeichnet Y dieKlassizierung vonX. Das heiÿt, dass die Verteilung PX vonX durch die bedingteVerteilung PX|Y festgelegt wird. Dieser Zusammenhang wird später durch die Bayesformel genauererläutert. Wie in der logistischen Klassizierung möchten wir einen Klassikator C konstruieren,der einer Realisierung X = x eine Klasse C(x) ∈
1, . . . ,K
zuordnet. Formal haben wir es also
mit einer deterministischen Abbildung
C : Rd 7→
1, . . . ,K
zu tun. In der Praxis wird C normalerweise über ein Trainingsample Xn =
(x1, y1), . . . , (xn, yn)
konstruiert, worauf wir später genauer eingehen werden. Ein wesentliche Frage bezüglich einesKlassikators ist seine Qualität hinsichtlich einer korrekten Klassizierung. Diese können wir überdie Wahrscheinlichkeit einer fehlerhaften Klassizierung P
(C(X) 6= Y
)beschreiben, und bezeich-
nen sie allgemein mit
R(C)
= P(C(X) 6= Y
).
R(C) entspricht also dem 0-1-Risiko. Ein optimaler Klassizierer Copt wäre daher (gegeben Mess-barkeit)
Copt = arg minC
R(C).
1Vielen Dank an Moritz Jirak für die Ausarbeitung dieses Abschnitts
40
Wie können wir diesen konstruieren? Zunächst ist es günstig das Risiko R(C)umzuformen:
R(C)
= P(C(X) 6= Y
)=
ˆP(C(x) 6= Y
∣∣X = x)PX(dx)
=
ˆ (1− P(C(x) = Y
∣∣X = x))PX(dx)
= 1−ˆ
P(C(x) = Y
∣∣X = x)PX(dx).
Wir sehen also, dass R(C)klein ist, wenn die bedingte Wahrscheinlichkeit P
(C(x) = Y
∣∣X = x)
möglichst groÿ ist. Bedingen auf Y liefert weiter
R(C)
= 1−ˆ ( K∑
k=1
P(C(x) = k
∣∣Y = k, X = x)P(Y = k
∣∣X = x))
PX(dx).
Nun benutzen wir die Tatsache, dass C deterministisch ist (eine leichte Verallgemeinerung istUnabhängigkeit von Y ). Dadurch erhalten wir
P(C(x) = k
∣∣Y = k, X = x)
= P(C(x) = k
)= 1C(x)=k,
und somit
R(C)
= 1−ˆ ( K∑
k=1
1C(x)=kP(Y = k
∣∣X = x))
PX(dx).
Wir haben es nun mit einer überraschend einfachen Optimierung zu tun. Um R(C)zu minimieren,
genügt es den Ausdruck
A(x) :=
K∑k=1
α(k, x)P(Y = k
∣∣X = x)
mit α(k, x) = 1C(x)=k,
für jedes x ∈ Rd zu maximieren. Es gilt nun allerdings α(k, x) ∈ 0, 1. Die Gröÿe A(x) ist folglichgenau dann maximal, wenn wir das meiste Gewicht (und somit α(k, x) = 1) auf
max1≤k≤K
P(Y = k
∣∣X = x)
legen (dies kann leicht bewiesen werden). A(x) ist also genau dann maximal, wenn wir für jedesx ∈ Rd
α(k, x) =
1, für k = k∗,
0, sonst,wobei k∗ = arg maxk=1,...,K P
(Y = k
∣∣X = x).
Dieses k∗ liefert uns automatisch die optimale Klassikation:
Copt(x)
= arg maxk=1,...,K P(Y = k
∣∣X = x). (4.1)
Theorem 4.5. Der deterministische Klassikator, welcher das Risiko R(C)minimiert, ist durch
die Klassikation in (4.1) gegeben und wird Bayesklassikator genannt.
Bemerkung 4.6.
(i) Die Optimalität der Bayesklassikation hängt essentiell mit der Denition des Risikos R(C)
zusammen.
(ii) Falls K = 2, erhalten wir eine sehr einfache Klassizierungsregel: Wenn P(Y = 1
∣∣X = x)≥
1/2, dann wählen wir Klasse k = 1, ansonsten Klasse k = 2.
Viele Klassikationsalgorithmen versuchen die Bayesklassikation zu imitieren. Ein allgemeinerZugang ist die bedingte Wahrscheinlichkeit P
(Y = k
∣∣X = x)anhand eines Trainingsamples zu
schätzen. Ein bekannter Repräsentant dieser Art ist das KNN-Verfahren (K-nearest neighbour,Übung 2). Ein anderes ist die lineare Diskrimimanzanalyse, die im nächsten Abschnitt behandeltwird.
41
4.3 Lineare Diskriminanzanalyse2
Das Problem der Bayesklassikation ist, dass es nicht so einfach ist, gute Schätzer für die bedingtenWahrscheinlichkeiten P
(Y = k
∣∣X = x)zu konstruieren. Allerdings kann die Bayesformel hier
helfen. Wir modellieren die Verteilung von X für jede Klasse k ∈ 1, . . . ,K mit K > 2 (alsogegeben Y ) durch eine Dichte
fk(x) = P(X = dx|Y = k)
und wählen a-priori-Wahrscheinlichkeiten der Klassen πk = P(Y = k) ∈ [0, 1] für k = 1, . . . ,K mit∑k πk = 1. Die Bayesformel liefert dann die a-posteriori-Zähldichte von Y
pk(x) = P(Y = k
∣∣X = dx)
=P(X = dx
∣∣Y = k)P(Y = k
)P(X = x
) =πkfk(x)∑Kl=1 πlfl(x)
.
Die Idee der linearen Diskrimimanzanalyse (LDA) ist nun, fk(x) als Gauÿdichte (univariat) zumodellieren, also
fk(x) =1√
2πσkexp
(− (x− µk)2
2σ2k
),
wobei µk und σk der Mittelwert und die Varianz der k-ten Klasse sind. Der Einfachheit halber seiσ2
1 = . . . = σ2K in der folgenden Diskussion. Dann erhalten wir
pk(x) =πk
1√2πσk
exp(− (x−µk)2
2σ2k
)∑Kl=1 πl
1√2πσl
exp(− (x−µl)2
2σ2l
) . (4.2)
Durch umformen erhalten wir, dass pk(x) genau dann maximal ist, wenn δk(x) maximal ist, gege-ben durch
δk(x) = xµkσ2− µ2
k
2σ2+ log
(πk). (4.3)
Tatsächlich sind die Werte πk, µk und σ2 aber unbekannt, und müssen geschätzt werden.
Methode 10: Lineare Diskriminanzanalyse. Wir denieren
πk =nkn, µk =
1
|nk|∑
j: yj=k
xj und σ2 =1
n−K
K∑k=1
∑j:yj=k
(xj − µk
)2,
wobei n die Gesamtanzahl des Trainingssamples Xn und nk die Anzahl des Trainingssamples inder k-ten Klasse sind. Dann ist der Klassizierer gegeben durch
C(x) = arg maxk=1,...,K δk(x) mit δk(x) = xµkσ2− µ2
k
2σ2+ log
(πk).
Im multivariaten Fall erhalten wir analog die Klassizierungsregel
δk(x) = x>Σ−1µk −1
2µ>k Σ−1µk + log
(πk), (4.4)
wobei Σ die d-dimensionale Kovarianzmatrix von X ∈ Rd ist, und µk ∈ Rd der Vektor der kompo-nentenweisen Erwartungswerte. Dabei können Σ, µ1, . . . , µK sowie π1, . . . , πK wieder über Plug-inbzw. relative Häugkeiten geschätzt werden.
Bemerkung 4.7. Eine weitere Verallgemeinerung stellt die quadratische Diskriminanzanalyse(QDA) dar, wo jede Klasse k eigene, im allgemeinen unterschiedliche Kovarianzmatrizen Σk be-sitzen. Dies führt zu einer quadratischen Klassizierungsregel.
2Vielen Dank an Moritz Jirak für die Ausarbeitung dieses Abschnitts
42
Obwohl die Motivation für die logistische Klassikation und LDA unterschiedlich ist, gibt eseinen engen Zusammenhang. Betrachten wir den Fall K = 2. Dann gilt p2(x) = 1−p1(x) und einekurze Rechnung ergibt für die LDA
log
(p1(x)
1− p1(x)
)= log
(p1(x)
p2(x)
)= c1 + c2x,
wobei die Konstanten c1, c2 von µ1, µ2 und σ2 abhängen. Im Fall der logistischen Klassizierunghaben wir:
log
(p1(x)
1− p1(x)
)= log
(p1(x)
p2(x)
)= β1 + β2x.
Der Unterschied liegt also nur in der Art und Weise, wie die Konstanten geschätzt werden! In derPraxis führt das oft zu sehr ähnlichen Ergebnissen, aber nicht immer.
5 Ausblick
Im letzten Teil der Vorlesung werden (voraussichtlich) noch folgende Themen behandelt:
(i) Modellwahl und statistisches Lernen
(a) Variablenselektion (Cp, AIC, BIC, R2)
(b) Lasso
(c) Dimensionsreduktion / Hauptkompontenanalyse
(ii) Resampling
(a) Bootstrap
(b) Kreuzvalidierung
Literatur
Agresti, A. and Finlay, B. (1997). Statistical Methods for Social Sciences. Prentice Hall.
Fahrmeir, L., Kneib, T., and Lang, S. (2009). Regression: Modelle, Methoden und Anwendungen.Springer, Berlin.
Georgii, H.-O. (2007). Stochastik. Einführung in die Wahrscheinlichkeitstheorie und Statistik. deGruyter, Berlin.
James, G., Witten, D., Hastie, T., and Tibshirani, R. (2013). An Introduction to Statistical Lear-ning (with Applications in R). Springer, New York.
Le, C. T. (2003). Introductory biostatistics. John Wiley & Sons.
Witting, H. (1985). Mathematische Statistik I. Teubner.
43
