skript zur vorlesung - meteo.uni-bonn.de · meteorologisches institut der universit at bonn skript...
TRANSCRIPT
Meteorologisches Institut der Universitat Bonn
Skript zur Vorlesung
Methoden der multivariaten Statistik
Sommersemester 2002
Andreas Hense
Version: April 2002 (mit Korrekturen Nov. 2005, PF)
1
Inhaltsverzeichnis
1 Einfuhrung 1
1.1 Matrix / Vektor Algebra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2 Multivariat normalverteilte Grundgesamtheiten 10
2.1 Multivariate Zufallsvariablen . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.1.1 Momente multivariater ZVA . . . . . . . . . . . . . . . . . . . . . . . 11
2.2 Multivariate Normalverteilung . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.3 Eigenschaften der Kovarianzmatrix . . . . . . . . . . . . . . . . . . . . . . . 13
2.4 Schatzungen der Kovarianzmatrix . . . . . . . . . . . . . . . . . . . . . . . . 16
2.5 Informationskomprimierung . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.6 Beispiele fur Guessmuster . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.7 Bedingte Wahrscheinlichkeitsdichten . . . . . . . . . . . . . . . . . . . . . . 25
3 Bivariate Grundgesamtheiten 28
3.1 Eigenschaften der Grundgesamtheit . . . . . . . . . . . . . . . . . . . . . . . 28
3.2 Schatzer der Eigenschaften der bivariaten GG . . . . . . . . . . . . . . . . . 31
3.3 Konfidenzintervalle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.4 Hypothesentest . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4 Multivariate Hypothesentests 35
4.1 Einfuhrung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
4.2 Die Hotelling T 2 Variable . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
4.3 Vergleich zweier multivariater Stichproben . . . . . . . . . . . . . . . . . . . 40
4.4 Die Macht des T 2 Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.5 Informationskomprimierung und der T 2 Test . . . . . . . . . . . . . . . . . . 42
4.6 Beispiel zum Hotelling T 2 Test . . . . . . . . . . . . . . . . . . . . . . . . . . 45
5 Diskriminanzanalyse 52
5.1 Das Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.2 Klassifikation bei multivariaten Normalverteilungen . . . . . . . . . . . . . . 54
5.3 Klassifikation bei unbekannter pdf . . . . . . . . . . . . . . . . . . . . . . . . 56
5.4 Beispiel zur Anwendung der Diskriminanzanalyse . . . . . . . . . . . . . . . 58
6 Empirische Orthogonalfunktionen 62
6.1 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
6.2 Schatzung der Prinzipalen Vektoren . . . . . . . . . . . . . . . . . . . . . . . 63
6.3 EOF und die Singularwert-Zerlegung SVD . . . . . . . . . . . . . . . . . . . 66
6.4 Schatzung der Konfidenzintervalle . . . . . . . . . . . . . . . . . . . . . . . . 67
6.5 Beweis des Satzes uber prinzipale Vektoren . . . . . . . . . . . . . . . . . . . 70
6.6 Rotation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
6.7 Singular Spectral Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
6.8 Beispiel fur eine EOF Analyse . . . . . . . . . . . . . . . . . . . . . . . . . . 75
7 EOF’s und ihre physikalische Interpretation 76
8 Kanonische Korrelationsanalyse 80
8.1 Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
8.1.1 Bedeutung der Metrikmatrizen . . . . . . . . . . . . . . . . . . . . . 85
8.2 Schatzung der kanonischen Korrelation . . . . . . . . . . . . . . . . . . . . . 87
8.3 Signifikanzschatzungen der kanonischen Korrelationen . . . . . . . . . . . . . 87
8.4 Beispiel fur eine CCA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
1 Einfuhrung
Abbildung 1 Feld der 2m Temperaturanomalie (Abweichung vom Mittelwert 1951-1980), Mittel-
wert fur Januar 1983, Nordhemisphare zwischen 17.5◦N und 77.5◦N, Einheit K
1 Einfuhrung
In der Vorlesung ”Statistik” sind im wesentlichen Methoden der univariaten Statistik be-
schrieben worden. Ein nicht unbedeutender Teil der Statistik in der Klimatologie befasst sich
aber mit Problemen wie etwa die Statistik des Bodendruckfeldes uber dem Nordatlantik, das
Feld der 2m Temperaturanomalie oder des 500 hPa Geopotentials uber der Nordhemisphare
oder das Windfeld in 200 hPa in den Tropen, um mal eine nicht sehr reprasentative Auswahl
zu liefern. Die Daten liegen typischerweise an Gitterpunkten vor, z.B. 72*13 = 936 fur die
Nordhemisphare, wenn man/frau einen Gitterpunktsabstand von 5◦ in N-S Richtung und
eine 5◦ Abstand in E - W Richtung wahlt und sich nur auf die Daten zwischen 17.5◦N und
77.5◦N konzentriert (vergl. Abb. (1)).
Deshalb konnen wir das diskrete Feld ”Temperatur” als eine Realisierung einer multi-
variaten ZVA mit der Dimension q = 936 auffassen oder allgemeiner bei der statistischen
Betrachtung von meteorologischen Daten in Feldform die Notwendigkeit einer multivariaten
Statistik ausmachen. Die Verallgemeinerung diskreter Felddaten auf kontinuierliche Feld-
daten ist ebenfalls moglich und fuhrt auf die Betrachtung von Zufallsfunktionen F (~r), die
wir aber nicht weiter betrachten wollen. Zwischen den Werten an den Gitterpunkten gibt
1
1 Einfuhrung
es Abhangigkeiten, die eine univariate Statistik nicht erfassen wurde. Diese Abhangigkeiten
entstehen aus der Dynamik: so sind z.B.
• in hyperbolischen Systemen Wellenausbreitungen (Rossbywellen, seismische Wellen)
gegeben, die verschiedene Punkte in Raum und Zeit verknupfen.
• in elliptischen Systemen (stationare Diffusion, Warmeleitung) die Randwerte, die die
Losung und damit die raumlich-zeitliche Struktur bestimmen
die entscheidenden Prozesse, die bei einer statistischen Analyse ein multivariates Vorgehen
erzwingen. Generell kann man sagen, dass man Statistik und Dynamik getrennt behandelt,
wenn man nur univariate Statistik betreibt. Will man aber eine mit der Dynamik konsis-
tente Statitik (statistisch-dynamische Analyse) betreiben, muss man zwangslaufig auf die
multivariate Statistik ausweichen. Das Paradebeispiel hierfur ist die numerische Analyse, in
der die multivariate Betrachtung der Beobachtungen an verschiedenen Raum-Zeitpunkten
eine mit der Dynamik konsistente Analyse ergeben soll. Schliesslich ist multivariate Statistik
gefragt, wenn generell Muster untersucht werden sollen. Das Paradebeispiel fur diese Art
von Anwendung ist der Nachweis und die Zuordnung von Klimaanderungssignalen: Ist
eine simuliertes Klimaanderungsmuster (z.B. ein Temperaturtrend) in den beobachteten
Temperaturtrends (vergl. Abb.(1)) wiederzufinden?
Wie wir in der Vorlesung ”Statistik” gesehen haben, sind wesentliche statistische Metho-
den formuliert worden fur univariate und normalverteilte Grundgesamtheiten. Es liegt daher
nahe, analog zu den univariaten Methoden, sich zunachst auf die multivariate Statistik nor-
malverteilter ZVA zu beschranken. Relevante Bucher in diesem Zusammenhang sind
• Morrison, D.F., Multivariate Statistical Methods, McGraw Hill Series in Probability
and Statistics ([9])
• Anderson, T.W., An Introduction to Multivariate Statistical Analysis, 2nd Edition, J.
Wiley & Sons, umfangreich und sehr mathematisch, sehr gut, ([10])
• Hartung, J. und Elpelt, B. Multivariate Statistik, Lehr- und Handbuch der angewand-
ten Statistik: Wie der Titel verrat, ein Buch zur Anwendung ohne umfangreiche ma-
thematische Herleitungen aber mit viel Anwendungsbeispielen ([16])
2
1 Einfuhrung
Abbildung 2 Feld des Trends der 2m Temperaturanomalie Winter DJF 1949-1998 zwischen 32.5◦S
und 77.5◦N, Einheit K/10a
• Hans von Storch und Francis Zwiers, Statistical Analysis in Climate Research, Cam-
bidge University Press: Das Buch, das viele statistische Verfahren, die in der Klimafor-
schung Verwendung finden, vorstellt und mit Beispielen erlautert. Mit dabei sind auch
die multivariaten Untersuchungen. [15]
Multivariate ZVA werden i.A. mit Vektoren im IRq identifizert. Deshalb werden wir zu
Beginn der Vorlesung die Matrix / Vektorrechenregeln zusammenstellen, auf die wir im
Rahmen der Vorlesung immer wieder zuruckgreifen werden.
1.1 Matrix / Vektor Algebra
Alle Vektoren in der Vorlesung werden als Spaltenvektoren aufgefasst, sofern es nicht explizit
anders erwahnt wird. D.h. ein q dimensionaler Vektor wird geschrieben als
~x =
x1
x2
...
xq
(1.1)
3
1 Einfuhrung
Entsprechende ist der transponierte Vektor ~xT ein Zeilenvektor.
~xT = (x1, x2, . . . , xq) (1.2)
Eine q × p dimensionale Matrix A ist eine Anordnung von p Spaltenvektoren der Dimension
q zu einem zweidimensionalen Array Ai,j
A =
A11 A12 . . . A1p
A21 A22 . . . A2p
......
...
Aq1 Aq2 . . . Aqp
(1.3)
Dabei bezeichnet der erste Index die Zeilennummer (1, . . . , q) und der zweite den Spalten-
index (1, . . . , p). Die Transponierte einer Matrix erhalten wir durch Vertauschen der Zeilen-
und Spaltenindices.
(AT )ij = Aji (1.4)
Eine Matrix nennt man/frau symmetrisch, wenn gilt
AT = A (1.5)
Die Spur einer Matrix ist definiert als die Summe uber die Diagonalelemente Aii
Sp(A) =
min(p,q)∑
i=1
Aii (1.6)
Als Diagonalmatrix wird eine Matrix bezeichnet, die nur auf der Diagonalen Eintrage hat,
die von Null verschieden sind
diag(γi, i = 1, q) =
γ1 0 . . . 0
0 γ2 0...
... 0 γq
(1.7)
Die Einheitsmatrix I ist eine Diagonalmatrix mit 1 als Diagonaleintragen.
I = diag(1, i = 1, q) (1.8)
Matrizen werden addiert, indem die Matrixelemente addiert werden
(A + B)ij = Aij + Bij (1.9)
4
1 Einfuhrung
Eine Matrix wird skalar mit c multipliziert, indem jedes Element mit dem Skalar c multipli-
ziert wird
(cA)ij = cAij (1.10)
Zwei Matrizen A, B der Dimension (p × r) und r × q werden multipliziert durch folgende
Rechenregel
(AB)ij =r∑
k=1
AikBkj (1.11)
Diese Multiplikation ist distributiv und assoziativ d.h. es gilt fur die drei Matrizen A(p ×r), B, C(r × q) bzw. C(q × s)
A(B + C) = AB + AC
A(BC) = (AB)C (1.12)
Dagegen ist die Matrixmultiplikation i.A. (d.h. auch fur quadratische Matrizen q = p = r = s
nicht kommutativ
AB 6= BA (1.13)
Lineare Gleichungssysteme fur die q Unbekannten xj der Form
q∑
j=1
Ai,jxj = bi (1.14)
konnen kompakt in Matrix - Vektorschreibweise geschrieben werden als:
A~x = ~b (1.15)
Zwei Vektoren konnen auf verschieden Weise multipliziert werden. Das Skalarprodukt zweier
Vektoren ~x, ~y der Dimension q ist definiert als
~xT ~y =
q∑
i=1
xiyi (1.16)
Da das Ergebnis ein Skalar ist, ist die skalare Multiplikation kommutativ
~xT~y = ~yT~x (1.17)
Das aussere Produkt zweier Vektoren ist eine Matrix A mit folgender Beziehung
A = ~x~yT
Aij = yixj (1.18)
5
1 Einfuhrung
Skalarprodukt und ausseres Produkt hangen uber folgende Beziehung zusammen:
Sp(~y~xT ) = ~xT~y (1.19)
Die Determinanten einer Matrix wollen wir mit det(A) bzw. |A| bezeichen.
Die Inverse einer (quadratischen) Matrix A wollen wir mit A−1 bezeichnen. Dies ist eine
Matrix B, fur die gilt
AB = BA = AA−1 = A−1A = I ⇒ B = A−1 (1.20)
Die Losung eines linearen Gleichungssystems A~x = ~b erhalten wir dann formal zu
~x = A−1~b (1.21)
Damit A−1 existiert, muss A nichtsingular sein bzw. den vollen Rang haben. Dies ist
gegeben, wenn gilt det(A) 6= 0. Existiert A−1 nicht, so heisst A singular und det(A) = 0.
Die Inverse der Transponierten einer Matrix ist die Transponierte der Inversen, sofern die
Matrix nichtsingular ist.
(A−1)T = (AT )−1 (1.22)
Die Inverse eines Matrixproduktes zweier nichtsingularer Matrizen A, B findet man/frau als
(AB)−1 = B−1A−1 (1.23)
Als Rang einer Matrix A mit den Dimensionen q × p bezeichnet man/frau die Anzahl der
linear unabhangigen Zeilen / Spaltenvektoren
rg(A) ≤ min(p, q) (1.24)
Ist bei einer quadratischen Matrix p = q der Rang gleich der Dimension rg(A) = q, so sagt
man/frau, die Matrix A habe den vollen Rang. In diesem Fall existiert eine Inverse A−1. Fur
die Rangberechnung von Matrizen gelten folgende Regeln
rg(AT ) = rg(A)
rg(AT A) = rg(A) = rg(AAT )
rg(A) = a rg(B) = b ⇒ rg(AB) ≤ min(a, b) (1.25)
Sei B eine quadratische Matrix mit vollem Rang p und A eine Matrix der Dimensionen q×p
mit rg(A) = a. Dann gilt
rg(AB) = a (1.26)
6
1 Einfuhrung
Die Multiplikation eines Vektors ~x und einer quadratischen, symmetrischen Matrix A von
der Form
~xT A~x (1.27)
nennt man/frau eine quadratische Form. Die Matrix A heisst positiv (semi-) definit, wenn
fur alle ~x 6= 0 gilt
~xT A~x > 0 positiv − definit
~xT A~x ≥ 0 positive − semidefinit (1.28)
Ist die quadratische Form kleiner oder kleiner gleich Null, so nennt man/frau A negati-
ve (semi-) definit. Ist kein Vorzeichen der quadratischen Form ausgezeichnet, so heisst A
indefinit.
~x ist ein Eigenvektor der quadratischen Matrix A, wenn sich durch die Multiplikation A~x
der Vektor bis auf einen Langenfaktor λ reproduziert.
A~x = λ~x (1.29)
Dies ist ein lineares, homogenes Gleichungssystem
(A − λI)~x = 0, (1.30)
das nur dann nichttriviale Losungen ~x 6= 0 haben kann, wenn gilt
det(A − λI) = 0 (1.31)
Diese Gleichung bestimmt die Nullstellen λ eines Polynoms (charakteristisches Polynom) der
Ordnung q
λ1, λ2 . . . λq (1.32)
wobei allerdings einige λi identisch seien konnen (”entartet”). Zu jedem dieser Eigenwerte
λi existiert ein Eigenvektor ~ei, so dass gilt
A~ei = λi~ei (1.33)
Ordnen wir die Eigenvektoren als Spaltenvektoren zu einer Matrix E an und schreiben die
Eigenwerte als Elemente einer Diagonalmatrix, erhalten wir
E = (~e1, ~e2, . . . , ~eq)
Λ = diag(λ1, . . . , λq)
AE = EΛ (1.34)
7
1 Einfuhrung
Ist A eine symmetrische Matrix, lehrt uns die lineare Algebra, dass gilt
λi ∈ IR (1.35)
Ist A ausserdem positiv (semi-) definit, so gilt
λi ≥ 0 (1.36)
fur alle i ∈ [1, q]. Die Eigenvektoren von symmetrischen Matrizen bilden ein vollstandiges,
orthogonales Basissystem im IRq, d.h. es gilt:
~eTi ~ej = δij
ET E = EET = I (1.37)
Betrachten wir zum Schluss noch Ableitungen von skalaren Funktionen nach Vektoren
(”Gradienten”). Sei f(x1, x2, . . . , xq) = f(~x) eine Abbildung aus dem IRq in IR. Dann konnen
wir nach den bekannten Rechenregeln die partiellen 1. und 2. Ableitungen bilden
∂f
∂xi
∂2f
∂xi∂xj
(1.38)
Sei f nun eine verallgemeinerte quadratische Form
f(~x) = (~a − C~x)T K(~a − C~x) (1.39)
wobei ~a ein fester q dimensionaler Vektor ist, K eine symmetrische Matrix und C eine
beliebige Matrix. Gesucht ist der Gradient von f bzgl. des Vektors ~x:
~∇xf =∂f
∂~x(1.40)
Ublich ware es, die verallgemeinerte quadratische Form als Doppelsumme auszuschreiben
und dann den Gradienten durch partielle Ableitung nach den Komponenten xi auszurechnen.
Man/frau kann sich aber viel Arbeit ersparen, wenn man/frau nach folgenden drei Regeln
verfahrt:
(1) Berechne die Ableitung ∂f
∂~xwie im skalaren Fall, d.h. unter Benutzung der ublichen
Rechenregeln zur Differentiation.
8
1 Einfuhrung
(2) Erhalte in jedem Fall die Reihenfolge der Matrix - Vektormultiplikationen, da diese
nicht kommutativ sind.
(3) Soll die quadratische Form nach dem Komponenten abgeleitet werden, die aus dem
transponierten Vektoranteil ~xT stammen, so verfahre nach den Regeln (1) und (2) ohne
Berucksichtigung der Transponierung und transponiere das Ergebnis anschliessend.
Mit Hilfe dieser Regeln erhalten wir dann
∂f
∂~x=
∂
∂~x(~aT K~a − ~aT KC~x − ~xT CT K~a + ~xT CT KC~x)
= −~aT KC − (CT K~a)T + (CTKC~x)T + ~xT CT KC
= −2~aT KC + 2~xT CT KC (1.41)
Auch die Matrix der zweiten Ableitungen (Hesse - Matrix) erhalten wir durch Anwendung
der Rechenregeln auf den Gradienten
H =∂2f
∂~x∂~xT= 2CT KC (1.42)
9
2 Multivariat normalverteilte Grundgesamtheiten
2 Multivariat normalverteilte Grundgesamtheiten
2.1 Multivariate Zufallsvariablen
Wir definieren eine q-dimensionale multivariate ZV ~X als Vektor
~X = (X1, X2, . . . , Xq)T ∈ � q.
Die Elemente des Vektors ~X sind kontinuierliche univariate ZV mit den Randdichten
f1(x1), f2(x2), . . . , fq(xq) und den Randverteilungen F1(x1), F2(x2), . . . , Fq(xq).
Die Verteilungsfunktion der multivariaten ZVA
~X = {( ~X, F (~x)), ~X ∈ � q}
ist definiert uber die Wahrscheinlichkeit
F (x1, x2, . . . , xq) = P (X1 ≤ x1, X2 ≤ x2, . . . , Xq ≤ xq).
Diese heisst auf die gemeinsame Verteilung oder ’joint distribution’.
Sei F (x1, x2, . . . , xq) kontinuierlich und differenzierbar, so konnen wir die gemeinsame
Wahrscheinlichkeitsdichtefunktion definieren als f(x1, . . . , xq) = F ′(x1, . . . , xq), mit
F (x1, . . . , xq) =
∫ x1
−∞
. . .
∫ xq
−∞
f(u1, . . . , uq)du1 . . . duq.
Unabhangigkeit: Wenn die Elemente der multivariaten ZVA ~X statistisch unabhangig
voneinander sind, dann gilt
f(x1, . . . , xq) = f1(x1) · . . . · fq(xq)
F (x1, . . . , xq) = F1(x1) · . . . · Fq(xq).
In diesem Fall konnten wir die univariate Statistik verwenden und mit den Randverteilungen,
bzw. -dichten arbeiten. Aber dies ist bei meteorologischen Feldern nicht der Fall!
Die Randverteilung eines Elementes eine multivariabten ZVA definiert sich als
fi(xi) =
∫ ∞
−∞
. . .
∫ ∞
−∞
f(x1, . . . , xq)dx1 . . . dxi−1dxi+1dxq. (2.1)
Durch die Integration der gemeinsamen Verteilung uber einzelne Elemente der multivariaten
ZVA wird die Verteilung unabhangig gemacht von diesen Elementen – dieser Prozess heisst
Marginalisierung.
10
2 Multivariat normalverteilte Grundgesamtheiten
Multivariate bedingte Verteilungen bzw. Dichten: Gegeben sei eine MV ZVA mit
f(x1, . . . , xq). Seien nun die Elemente xp + 1, . . . , xq z.B. durch ein Experiment oder Be-
obachtung festgelegt, dann lautet die auf diese Elemente bedingte Wahrscheinlichkeitsdichte
f(x1, . . . , xp|xp+1, . . . , xq) =f(x1, . . . , xq)
f(xp+1, . . . , xq). (2.2)
2.1.1 Momente multivariater ZVA
Der Erwartungswert E( ~X) ist ein Vektor bestehend aus den Erwartungswerten der einzelnen
Komponenten von ~X
E( ~X) = (E(X1), . . . , E(Xq))T = ~µ.
Um auf die allgemeine Form der multivariaten Varianz zu kommen, definieren wir zuerst die
Kovarianz
Cov(Xi, Xj) = E ((Xi − µi) (Xj − µj))
=
∫ ∞
−∞
. . .
∫ ∞
−∞
(xi − µi)(xj − µj)f(x1, . . . , xq)dx1 . . . dxq
=
∫ ∞
−∞
∫ ∞
−∞
(xi − µi)(xj − µj)f(xi, xj)dxidxj
= σ2ij (2.3)
Wenn i = j ist, dann ist σ2ii = σ2
i die Varianz der i-ten Komponente. Bei einer q-
dimensionalen ZVA ~X existieren also q×q Kovarianzen σ2ij, wobei diese naturlich symmetrisch
sind. Wir konnen diese also zusammenfassen zu der Kovarianzmatrix Σ
E(( ~X − E( ~X))( ~X − E( ~X))T
)= Σ. (2.4)
Die Kovarianzmatrix ist also eine symmetrische positiv definite Matrix.
2.2 Multivariate Normalverteilung
Wir bereits erwahnt wurde, werden wir uns bei der multivariaten Statistik auf ZVA be-
schranken, die multivariat normalverteilt sind. Das hat zwei Grunde:
1. Ein Zufallsvektor, der entsteht aus der Summe eine grosseren Anzahl unabhangig und
identisch verteilter Zufallsvektoren ist nach dem zentrale Grenzwertsatz der Statistik asym-
ptotisch multivariat normalverteilt (dies gilt analog zu univariaten auch fur multivariate
11
2 Multivariat normalverteilte Grundgesamtheiten
ZVA). Eine solche Annahme kann in vielen Fallen gemacht werden, sollte aber immer be-
grundet oder getestet werden.
2. Ohne die Annahme der Normalverteilung waren statistische Tests sehr komplex und even-
tuell nur uber nicht parametrische Methoden zu losen.
Sei
~X = {(~x, f(~x)), ~x ∈ � q} (2.5)
eine q-dimensionale ZVA. ~x heißt multivariat NV, wenn f(~x) die Form
f(~x) =1
Zexp(−1
2(~x − ~µ)tB(~x − ~µ)) (2.6)
hat, wobei B eine symmetrische, positiv-definite Matrix ist (d.h. alle Eigenwerte sind positiv)
und Z der Normierungsfaktor. Bedenke, daß eigentlich f(~x) = f(~x, ~µ,B)! Diese multivariat
NV ZVA ist symmetrisch um ~µ, d.h.
∫ ∞
−∞
. . .
∫ ∞
−∞
(~x − ~µ)f(~x, ~µ,B) dx1dx2...dxq = ~0 (2.7)
Damit ist aber
E(~x − ~µ) = ~0
⇒ E(~x) = ~µ (2.8)
In der Bestimmungsgleichung von f(~x, ~µ,B) war B noch unbestimmt. Daher bildet man nun
~∇µ
∫ ∞
−∞
. . .
∫ ∞
−∞
(~x − ~µ) f(~x, ~µ,B) dx1dx2...dxq = ~0 (2.9)
Ausrechnen der Ableitung fuhrt auf (I ist die Einheitsmatrix, B = Bt)
∫ ∞
−∞
. . .
∫ ∞
−∞
(I − (~x − ~µ)(~x − ~µ)tB) f(~x, ~µ,B) dx1dx2...dxq = ~0 (2.10)
Damit wiederum gilt auch (O ist die Nullmatrix)
E(I − (~x − ~µ)(~x − ~µ)tB) = O (2.11)
und daraus folgend
E((~x − ~µ)(~x − ~µ)tB) = E(I) = I (2.12)
d.h. die Matrix B ist die Inverse der Kovarianzmatrix Σ. Damit ist im Fall der multivariaten
NV die gesamte Verteilung durch die Parameter ~µ und Σ vollstandig beschrieben.
12
2 Multivariat normalverteilte Grundgesamtheiten
2.3 Eigenschaften der Kovarianzmatrix
In Kapitel 4.4 der Vorlesung Statistik war die pdf (”probability density function”) einer
multivariaten, normalverteilten ZVA ~X mit der Dimension q definiert worden als:
f(~x) =1
Zexp(−1
2(~x − ~µ)T Σ−1(~x − ~µ)) (2.13)
wobei Z =√
2πq det(Σ) der Normierungsfaktor ist.
Die Kovarianzmatrix Σ ist eine positiv definite, symmetrische Matrix. Die Eigenvektoren
~ej bilden deshalb eine vollstandige orthonormale Basis im IRq und die Eigenwerte λj sind
reell und positiv definit. Die Eigenwerte und die Eigenvektoren seien so angeordnet, daß gilt:
λ1 ≥ λ2 ≥ . . . λq (2.14)
Dann lassen sich die Eigenwertgleichungen zusammenfassen zu der Matrixgleichung
ΣE = EΛ (2.15)
wobei die Spalten der Matrix E die Eigenvektoren ~ej sind und die Matrix Λ eine Diago-
nalmatrix ist mit den Eigenwerten λj als Diagonalelemente. Die Orthonormalitatsrelation
lautet dann
ET E = I (2.16)
wobei I die Einheitsmatrix und ET die Transponierte von E ist. Wegen der Vollstandigkeit
kann jeder Vektor ~x nach den Eigenvektoren entwickelt werden. Faßt man/frau die Entwick-
lungskoeffizienten zum Koeffizientenvektor ~a zusammen, so kann man/frau schreiben
~x =
q∑
i=1
axi ~ei = E~ax ~µ =
q∑
i=1
aµ~ei = E~aµ
(~x − ~µ) = E(~ax − ~aµ) (2.17)
Der Mahalonobisabstand D2
D2 = (~x − ~µ)T Σ−1(~x − ~µ) (2.18)
kann man/frau dann auch schreiben als:
D2 = (~ax − ~aµ)T ET Σ−1E(~ax − ~aµ) (2.19)
13
2 Multivariat normalverteilte Grundgesamtheiten
Wenn ~ej ein Eigenvektor von Σ zum Eigenwert λj ist, ist ~ej ebenfalls auch Eigenvektor von
Σ−1 zum Eigenwert λ−1j . Dann folgt wegen der Orthonormalitatsrelation
D2 = (~ax − ~aµ)TΛ−1(~ax − ~aµ)
=
q∑
j=1
(axj − aµ
j )2
λj
(2.20)
Fuhren wir diese Koordinatentransformation fur das Integral zur Berechnung der Normie-
rungsfaktors Z in Gl. (2.13) ein, so folgt fur das Volumenelement dqx die Relation
dqx = det(E)dqax = dqax (2.21)
oder
Z =
∫ +∞
−∞
. . .
∫ +∞
−∞
exp(−1
2
q∑
j=1
(axj − aµ
j )2
λj
)dqax
=
q∏
j=1
∫ +∞
−∞
exp(−1
2
(axj − aµ
j )2
λj
)daj (2.22)
Die Integrale unter dem Produkt haben jeweils den Wert√
2πλj. Dann erhalt man/frau
fur den Normierungsfaktor
Z = (2π)q
2
√√√√q∏
j=1
λj = (2π)q
2
√det Σ (2.23)
Wie man/frau sieht, sind die Eigenvektoren der Kovarianzmatrix offenbar ausgezeichnete
Richtungen im Vektorraum IRq, namlich genau die, fur die die Kovarianzen cov(aiaj), i 6= j
verschwinden (siehe Darstellung der Mahalanobisdistanz durch die Entwicklungskoeffizienten
a). Geometrisch ist der Fall q = 2 besonders einfach zu interpretieren. Sei die Kovarianzma-
trix Σ gegeben als
Σ =
σ2
1 c
c σ22
(2.24)
Setzt man/frau
c = σ1σ2ρ (2.25)
wobei ρ der Korrelationskoeffizient zwischen ersten und zweiten Komponente der ZVA ~X =
(X1, X2) und σ1, σ2 entsprechenden Standardabweichungen sind, so erhalt man/frau fur die
Mahalanobisdistanz folgenden Ausdruck:
D2 =(x1 − µ1)
2
σ22(1 − ρ2)
− 2ρ(x1 − µ1)(x2 − µ2)
σ1σ2(1 − ρ2)+
(x2 − µ2)2
σ21(1 − ρ2)
(2.26)
14
2 Multivariat normalverteilte Grundgesamtheiten
-3 -2 -1 0 1 2 3-3
-2
-1
0
1
2
32-dim Gauss pdf, Eigenvektoren der Kovarianzmatrix
Abbildung 3 Illustration zum zweidimensionalen normalverteilten Problem
Eine Kontur konstanter Wahrscheinlichkeitsdichte f0 ist dann durch die Gleichung
D2(x1, x2) + 2 ln(f0Z) = 0 (2.27)
gegeben. Dies ist die Gleichung einer Ellipse (weil σ21σ
22(1 − ρ2) > 0 ) mit dem Mittelpunkt
( µ1, µ2). Die große Halbachse bildet einen Winkel α mit der x-Achse mit
tan 2α = −2ρσ1σ2
(σ21 + σ2
2)(2.28)
Die Drehung der Koordinatenachsen (x1, x2) in Richtung der Eigenvektoren von Σ dreht
diese Ellipse in Normalform:
D2(a1, a2) + 2 ln(f0Z) = 0 (2.29)
mit
D2(a1, a2) =(a
(x)1 − a
(µ)1 )2
λ1
+(a
(x)2 − a
(µ)2 )2
λ2
(2.30)
und
λ1,2 =σ2
1 + σ22
2±√
(σ21 − σ2
2)2
4+ σ2
1σ22ρ
2 (2.31)
Die Ellipse liegt nun mit den Hauptachsen parallel zu den Koordinatenachsen (”Haupt-
achsentransformation”), die Langen der Hauptachsen sind proportional zu den Eigenwerten
15
2 Multivariat normalverteilte Grundgesamtheiten
(vergl. Abb. (3)). Im Fall q > 2 betrachtet man/frau Flachen (Hyperflachen der Dimension
q−1) konstanter Wahrscheinlichkeitsdichte. Fur Normalverteilungen sind diese Hyperflachen
Ellipsoide. Dann spannen die Eigenvektoren der Kovarianzmatrix dieses Ellipsoid auf. Die
Lange der q Ellipsoidhalbachsen ist proportional zu den Eigenwerten. Wir werden in ei-
nem spatern Abschnitt sehen, daß die Eigenvektoren der Kovarianzmatrix noch eine weitere,
bemerkenswerte Eigenschaft haben, die der effektivsten Darstellung.
2.4 Schatzungen der Kovarianzmatrix
Analog zur Stichprobe einer univariaten ZVA kann auch die Stichproben ZVA einer multiva-
riaten Grundgesamtheit eingefuhrt werden. Die Stichproben ZVA einer univariaten ZVA bei
einer Stichprobenlange m war eine Vektorwertige ZVA ~X(m) der Dimension m. Im Fall einer
multivariaten Grundgesamtheit mit der Dimension q > 1 ist die Stichproben ZVA eine q×m
dimensionale Matrix X , deren Spalten jeweils ein Stichprobenelement eines q-dimensionalen
Vektors enthalten und deren Zeilen der Lange m univariate Stichprobenvektoren der i-ten
Komponente der Vektor ZVA enthalten. Sei die Grundgesamtheit, der die Stichprobe entnom-
men wurde, multivariat normalverteilt mit den Parametern ~µ und Σ. Es gilt nun, Schatzer
fur Erwartungswert und Kovarianzmatrix anzugeben. In der Vorlesung Statistik hatten wir
die Maximum Likelihood Methode eingefuhrt als eine allgemeine Vorschrift zur Konstruktion
von Schatzern. Dies wollen wir jetzt anwenden.
Mit Hilfe der Definition (2.13) ist die Wahrscheinlichkeit, die m Stichprobenelemente ~xk
zu beobachten, gegeben durch
f(~x1, . . . , ~xm) =1
((2π)qm(detΣ)m)12
exp
(−1
2
m∑
k=1
(~xk − ~µ)T Σ−1(~xk − ~µ)
)(2.32)
Bilden wir die log-likelihood Funktion l, erhalten wir
l = ln f = −qm
2ln 2π − m
2ln detΣ − 1
2
m∑
k=1
(~xk − ~µ)T Σ−1(~xk − ~µ) (2.33)
Fuhren wir nun das arithmetische Mittel ~x der Stichprobe ein,
~x =1
m
m∑
k=1
~xk, (2.34)
16
2 Multivariat normalverteilte Grundgesamtheiten
so konnen wir die log-likelihood Funktion umschreiben zu
l = −qm
2ln 2π − m
2ln detΣ − 1
2
m∑
k=1
(~xk − ~x)T Σ−1(~xk − ~x)
︸ ︷︷ ︸I
+m
2(~x − ~µ)T Σ−1(~x − ~µ)
︸ ︷︷ ︸II
(2.35)
Mit Hilfe der Rechenregeln fur Vektormultiplikationen konne wir den mit I bezeichneten
Term umschreiben zu
m∑
k=1
(~xk − ~x)T Σ−1(~xk − ~x)
=
m∑
k=1
(Σ−1(~xk − ~x))T (~xk − ~x) = Sp(
m∑
k=1
(~xk − ~x)(~xk − ~x)T Σ−1) = Sp(AΣ−1) (2.36)
wobei A die sogenannte Matrix der Produktsummen der Abweichungen vom Mittelwert ist.
Dann lautet die log-likelihood Funktion
l = −qm
2ln 2π − m
2ln detΣ − 1
2Sp(AΣ−1) − m
2(~x − ~µ)T Σ−1(~x − ~µ) (2.37)
Ableiten der log - likelihood Funktion (siehe Rechenregeln oben) nach dem Vektor ~µ und
Nullsetzen (maximum likelihood) fuhrt auf den Schatzer ~µ fur den Erwartungswert
m(~x − ~µ)Σ−1 = 0
~xT Σ−1 = ~µT Σ−1
=⇒ ~µ = ~x (2.38)
d.h. der arithmetische Mittelwert ist der ML Schatzer fur den Erwartungswert. Setzen wir
diesen wieder in die log - likelihood Funktion ein, so ist die quadratische Form II in Gl. (2.35)
identisch Null und wir mussen, um den ML Schatzer fur die Kovarianzmatrix zu bestimmen,
nur noch die Funktion
l = −qm
2ln 2π − m
2ln detΣ − 1
2Sp(AΣ−1) (2.39)
betrachten. Mit Hilfe der Eigenvektorzerlegung fur die (positiv definite, symmetrische) Ko-
varianzmatrix Σ ist es moglich zu zeigen, dass gilt
ln detΣ = Sp(ln Σ) (2.40)
Dann erhalten wir
l = Sp(−m
2lnΣ − 1
2AΣ−1) (2.41)
17
2 Multivariat normalverteilte Grundgesamtheiten
l ist damit eine Funktion der unbekannten Matrix Σ. Spurbildung und Ableitung nach Σ
sind vertauschbare Operationen, so dass wir als Bestimmungsgleichung fur den Schatzer der
Kovarianzmatrix erhalten:
−m
2Σ−1 +
1
2AΣ−2 = 0 (2.42)
Bezeichnen wir mit S den Schatzer der Kovarianzmatrix, so erhalten wir
−mS + A = 0
S =1
mA (2.43)
als ML Schatzer der Kovarianzmatrix. Aufgelost nach den Komponenten ergibt sich
Sij =1
m
m∑
k=1
(xik − xi)(xjk − xj)
xi =1
m
m∑
k=1
xik (2.44)
Definieren wir
di,k = xi,k − xi (2.45)
als ein Matrixelement fur die sogenannte Datenmatrix D der zentrierten Abweichungen, so
ist der ML Schatzer der Kovarianzmatrix gegeben als
S =1
mDDT (2.46)
Wie im univariaten Fall ist jedoch Σ nur asymptotisch unverzerrt, ein unverzerrter Schatzer
ist gegeben durch
S∗ =1
m − 1DDT (2.47)
Die Schatzer fur Erwartungswert und Kovarianzmatrix sind selbstverstandlich wiederum
ZVA. So ist der Schatzer des E-Wertes ~x eine multivariat normalverteilte ZVA mit Erwar-
tungswert ~µ und Kovarianzmatrix Σ/m. Im univariaten Fall war die pdf der aus dem Schatzer
der Varianz σ2 abgeleiteten Große (m− 1)σ2/σ2 die χ2 Verteilung mit (m− 1) Freiheitsgra-
den. Die multivariate Verallgemeinerung fur den unverzerrten Schatzer der Kovarianzmatrix
ist die sogenannte Wishart - Verteilung: eine pdf fur Matrizen!
Sei D eine Stichproben- (Daten-) matrix – Stichprobenlange m und Vektordimension q –
einer multivariat normalverteilten Grundgesamtheit mit Erwartungswert null und Kovari-
anzmatrix Σ (≡ N(~0, Σ)). Dann ist die pdf der positiv definiten symmetrische Matrix
A = (m − 1)S∗ = DDT (2.48)
18
2 Multivariat normalverteilte Grundgesamtheiten
die Wishartverteilung W (A, Σ, m − 1) mit m − 1 Freiheitsgraden (Γ(n) ist die im Kapitel 5
der Statistik I definierte Γ- Funktion):
W (A, Σ, m − 1) =(det A)
m−q−22 exp(−1
2Sp(Σ−1A))
2(m−1)q
2 πq(q−1)
4 (det Σ)m−1
2
∏q
i=1 Γ(m−i2
)(2.49)
Im Fall q = 1 ist die Wishart Verteilung identisch zur χ2 Verteilung.
Zur kompletten Beschreibung der multivariaten Normalverteilung ist aber nicht eine
Schatzung von Σ gefragt, sondern eine von Σ−1. Es ergibt sich also die Frage, ob aus unserer
ML Schatzung fur Σ eine Inverse berechnet werden kann. Die Matrix S (oder auch S∗) kann
invertiert werden, wenn rg(S) = rg(Σ) = q ist. Aus der Gleichung fur den Schatzer und der
Ungleichung fur den Rang von Matrixprodukten folgt
rg(S) = rg(DDT ) ≤ rg(D) ≤ min(m − 1, q) (2.50)
Ist also m ≥ q, so ist die geschatzte Kovarianzmatrix invertierbar, ist dagegen m < q, ist
die geschatzte Kovarianzmatrix singular und damit auch keine Schatzung der pdf moglich.
In diesem Fall ist auch die Kovarianzmatrix nicht mehr positiv definit, sondern nur noch
semidefinit d.h. mindestens ein Eigenwert ist identisch Null. Diese Matrix ist dann auch
nicht mehr Wishart verteilt! Dieses Problem nennt man den Fluch der Dimensionen (Curse
of Dimension). Es lieget darin begrundet, dass die Dichte an Beobachtungspunkten in einem
q dimensionalen Raum mit ∆x−q abnimmt, wenn ∆x ein typischer eindimensionaler Abstand
zwischen den Datenpunkten ist, z.B. der mittlere Zweipunkteabstand
∆x ∼ 1
m2
m∑
k,k′=1
√(~xk − ~xk′)2 (2.51)
Auf der anderen Seite ist der Fall m < q in der Klimatologie eigentlich der ubliche
Fall. Denn wenn man die Kovarianzmatrix der Anomalien der 2m Temperaturs zwischen
32,5◦S und 77.5◦N fur einen bestimmten Monat (z.B. Januar) auf einem 5 × 5◦ Gitter
(q ' 1656) schatzen will, so liegen nur 118 Januarfelder vor (1881 - 1998) (das ist schon sehr
viel im Vergleich mit anderen Daten). Eine komplette Darstellung der (angenommenen)
multivariaten Normalverteilung der ZVA Temperaturanomalien im Januar im Nordatlantik
an allen Gitterpunkten erscheint nicht moglich. Wir werden aber im nachsten Kapitel sehen,
wie man da Abhilfe schaffen kann.
19
2 Multivariat normalverteilte Grundgesamtheiten
Sei fur die weitere Diskussion angenommen m > q, d.h. die Schatzung der Kovarianzmatrix
fuhrt auf eine nichtsingulare Matrix S∗ und wir konnen (S∗)−1 berechnen. Die pdf der Matrix
ZVA (S∗)−1 ist ebenfalls bekannt und heißt die inverse- Wishart Verteilung: Sei A eine
Wishart verteilte Matrix ZVA mit W (A, Σ, m − 1). Dann ist die Matrix B = A−1 verteilt
mit der pdf
W−1(B, Ψ, m − 1) =det(Ψ)
m−12 (det B)
m+q
2 exp(−12Sp(ΨB−1))
2(m−1)q
2 πq(q−1)
4
∏q
i=1 Γ(m−i2
)(2.52)
wobei Ψ = Σ−1 gesetzt wurde. Daraus laßt sich der Erwartungswert fur den Schatzer (S∗)−1
berechnen zu
E((S∗)−1) =m − 1
m − q − 1Σ−1 (2.53)
Obwohl also der Schatzer der Kovarianzmatrix unverzerrt ist, ist die Inverse des Schatzers
eine (moglicherweise stark) verzerrte Schatzung der Inversen der Kovarianzmatrix. Die Ver-
zerrung ist umso großer, je schlechter die Bedingung m > q erfullt ist. Eine unverzerrte
Schatzung der inversen Kovarianzmatrix bildet man aus der inversen geschatzten Kovari-
anzmatrix dann wie folgt:
Σ−1 =m − q − 1
m − 1(S∗)−1 (2.54)
Mit dieser Schatzung ist auch eine unverzerrte Schatzung der Mahalanobisdistanz moglich:
D2 = (~x − ~µ)T Σ−1(~x − ~µ) (2.55)
wohingegen die intuitive Schatzung
D2 = (~x − ~µ)T (Σ∗)−1(~x − ~µ) (2.56)
positiv verzerrt ist, d.h. zu große Abstande liefert, da m−1m−q−1
> 1 ist. Die Verzerrung ist im
ubrigen ein Folge der Schatzung der q Parameter des Erwartungswertes. Setzt man/frau
diesen Erwartungswert als bekannt voraus, so verschwindet die Verzerrung der inversen
geschatzen Kovarianzmatrix.
2.5 Informationskomprimierung
Betrachten wir jetzt den Fall m < q. Auf den ersten Blick scheint man hier nichts machen zu
konnen, wie auch folgendes Zitat aus dem amerikanischen Journal of Atmospheric Sciences
(Chervin, 1981, Vol. 38, p.888) zu belegen scheint:
20
2 Multivariat normalverteilte Grundgesamtheiten
• .., for GCM fields with several thousands grid points, practically unattainable amounts
of computer time and millennia of observations would be required to validate a complete
field ..
Dem ist naturlich nicht so und die Abhilfe lautet:
• Wenn die Stichprobenlange zu kurz ist im Vergleich zur Vektordimension m < q,
mussen wir die Vektorlange adaquat verkurzen: d.h. die Information, die durch q
Datenpunkte dargestellt wird, auf deutlich weniger Datenwerte zu kompri-
mieren.
Diese nennt man/frau auch die ”a-priori” Datenreduktion bzw. allgemeiner die ”a-priori”
Reduktion der raumlichen Freiheitsgrade. Dies ist ein ganz allgemeines Problem der mul-
tivariaten Statistik, der bereits oben erwante Fluch der Dimensionen. Deshalb gibt es z.B.
auch ganze Sonderforschungsbereiche der DFG, die sich mit der Dimensionsreduktion in
komplexen Systemen befassen.
Man laßt naturlich nicht einfach irgendwelche Gitterpunkte wegfallen, dies ware ein etwas
unsystematischer Weg zur Reduktion der Freiheitsgrade, sondern spezifiziert eine Abbildung
H aus dem IRq in den Vektorraum IRq mit q < m � q.
H(~x) = ~y (2.57)
und
dim(~y) = q (2.58)
Diese Abbildungsfunktion H zusammen mit der Festlegung der Dimension q ist die mathe-
matische Zusammenfassung des eingesetzten ”a-priori” Wissens. Beispiel:
• Sei q = 1. Dann ist eine gultige a-priori Datenreduktion die Berechnung des raumlichen
Mittelwertes (Flachenmittelwert) der Vektor ZVA ~X = (X1, . . . , Xq) mit
H =
q∑
j=1
wjXj (2.59)
wobei die wj noch eine raumliche Wichtung (etwa proportional zum Cosinus der geo-
graphischen Breite oder zu der Anzahl der Meßpunkte) darstellt.
21
2 Multivariat normalverteilte Grundgesamtheiten
Besonders beliebt sind naturlich lineare Abbildungen H, d.h. Matrizen der Dimension q×q.
In diesem Fall laßt sich das Datenreduktionproblem auch anders darstellen. Als Vorwissen
definiert man/frau dann sogenannte Guess Vektoren (”Ratevektoren”) ~gk, k = 1, q, die die
ZVA Vektoren ~X des IRq moglichst gut darstellen sollen. Diese Guess-Vektoren sind ebenfalls
Vektoren der Dimension q, sie lassen sich im Fall, daß horizontale Felder ~X betrachtet werden,
wie diese horizontalen Felder als Isolinienkarten darstellen. Die zu betrachtenden Vektor ZVA
~X , d.h. die Karten des Feldes sollen dann so aussehen wie die Guessmuster ~gk, deshalb der
Name. Gesucht sind dann Amplituden ak der Guessmuster ~gk, so daß gilt
~X =
q∑
k=1
~gkak + ~r (2.60)
Der Vektor ~r zeigt an, daß wegen der Einschrankung q � q der Vektor ~X naturlich nicht
vollstandig dargestellt werden kann. Dieses Restglied der Entwicklung stellt den Anteil des
ursprunglichen Datenvektors dar, den man/frau wegen der Einschrankungen der multivaria-
ten Statistik weglassen muß. In ~r werden die Anteile zusammengefaßt, die nicht analysiert
werden konnen. Die Entwicklungskoeffizienten ak sind unbekannt, ebenfalls auch der Vektor
~r. Da man/frau aber moglichst wenig weglassen will und moglichst viel statistisch analy-
sieren will, liegt es nahe, die Forderung aufzustellen, daß die Lange von ~r moglichst klein
ist:
|~r|2 = ~rT~r!= min (2.61)
Der Residuumsvektor ~r ist gegeben durch
~r = ~X −q∑
k=1
~gkak (2.62)
Die Lange |~r|2 ist dann
|~r|2 = ( ~X −q∑
k=1
~gkak)T ( ~X −
q∑
k=1
~gkak) = F(ak) (2.63)
Da ~X bekannt ist, ist die Lange des Residuumsvektor damit eine Funktion F(a‖) der unbe-
kannten Koeffizienten ak. Wegen der Minimumsbedingung folgt dann
∂F∂ak
= 0 ; k = 1, q (2.64)
Dieses sind q Bestimmungsgleichungen fur die Unbekannten ak. Die Losung dieses sogenann-
ten ”least squares” Verfahren (wg. der Forderung Lange2 = Minimum) laßt sich elegant und
22
2 Multivariat normalverteilte Grundgesamtheiten
kompakt in Matrizenschreibweise finden. Die Entwicklungsgleichung (2.60) kann man auch
schreiben als:
~X = G~a + ~r (2.65)
wobei die Spalten der q × q dimensionalen Matrix G die Guessmuster ~gk sind. Einsetzen in
die Minimalfunktion F und ableiten nach den unbekannten Koeffizienten ak fuhrt auf ein
lineares Gleichungssystem
(GT G)~a = GT ~X (2.66)
mit der Losung
~a = (GT G)−1GT ~X (2.67)
wobei die Inverse der Matrix GT G gebildet werden kann, wenn die Guessvektoren ~gk linear
unabhangig sind. Die Matrix (GT G)−1GT ist dann die gesuchte Abbildung H. Die Vektoren
~Xfil = G~a = G(GT G)−1GT ~X (2.68)
stellen raumlich gefilterte Felder der ursprunglichen Felder ~X dar, die eine wohldefinierte An-
zahl von raumlichen Freiheitsgraden besitzten (namlich q viele, da genau q viele Amplituden
ak notwendig sind, das gefilterte Feld ~Xfil darzustellen. Fur die Vektoren ~a gilt nun m > q,
so daß man/frau jetzt fur diese Vektoren die Parameter der multivariaten Normalverteilung
selbst schatzen kann. Dies ist dann die Normalverteilung in dem Unterraum des ursprung-
lichen Vektorraums IRq, der von den q Guess - Vektoren aufgespannt wird. Dadurch blickt
man/frau bei der Analyse der multivariaten Daten nur in Richtung der a-priori festgelegten
Guess Vektoren. Da erkennt man/frau auch sofort das Dilemma: sind die Guess- Vektoren
schlecht vorgeschrieben, so blickt man / frau eventuell an den statistischen Geschehnissen
vorbei und analysiert nur Rauschen, obwohl in anderen Raumrichtungen sehr wohl etwas
statistisch interessantes ablauft. Die wahre Kunst der multivariaten Statistik ist nicht das
Schatzen und das Anwenden von Testverfahren etc., sondern die a-priori Datenreduktion.
2.6 Beispiele fur Guessmuster
Oben haben wir bereits ein Beispiel fur ein geometrisch motiviertes guess Muster zur Da-
tenreduktion gesehen: die Bildung raumlicher Mittelwerte. Weitere Beispiel, die sich aus
Geometrieuberlegungen ergeben sind Fourierzerlegung der Felder in zwei- oder dreidimen-
sionale harmonische Funktionen bei einem rechteckigen Gebiet. Je nach Randbedingungen
23
2 Multivariat normalverteilte Grundgesamtheiten
muss man dann sich nur auf sin oder cos oder aber auch auf beide Formen konzentrieren.
Die Werte dieser Funktionen an den verschiedenen Gitterpunkten sind die Spalteneintrage
in der Matrix G. Ublicherweise numeriert man den Vektor der Amplituden ~a so durch, dass
eine Hierarchie in der Skala entsteht. Die Amplituden fur die grossten Skalen sind die ersten
Eintrage in ~a. Allerdings ist auch jede andere Form der Sortierung zulassig, die sich aus
dem betrachteten Problem ergibt. Auf der Kugel wurde man Kugelfunktion Y ml wahlen, die
charakterisiert sind durch die Eigengleichung
~∇2Y ml = − l(l + 1)
a2Y m
l (2.69)
Eine Weiterfuhrung dieser Idee der Eigenfunktionen legt es nahe, z.B. die Eigenfunktionen
von einfachen, linearen Modellen als Basismuster zu verwenden. In diesem fall kammt man
dem Ziel der Verschmelzung von Statistik und Dynamik ziemlich nahe. Geeignete Kandidatin
ist z.B. die zweidimensionale (im physikalischen Raum) Advektions-Diffusiongleichung bei
einem gegebenen Stromungsfeld ~v fur eine Statthaltervariable Ψ
∂
∂tΨ + ~∇(~vΨ) − ~∇(K~∇Ψ) = 0 (2.70)
Diskretisiert man diese Gleichung in geeigneter Weise z.B. auf einem Gitter durch finite
Differenzen, so erhalt man eine lineare Gleichung der Form
d
dt~Ψ + M(~v, K)~Ψ = 0 (2.71)
Die Systemmatrix M hangt von dem vorgegebenen Parametern ~v unf K ab, so wie von den
Details der Diskretisierung. Der Vektor ~Ψ entsteht durch eine geeignete Numerierung der
Gitterpunkte. Als guess Muster bieten sich nun z.B. die Eigenvektoren der Matrix M an, die
i.A. komplex sind, da die Matrix M nicht notwendigerweise symmetrisch ist. Ist die Matrix
jedoch reell, so kann man auch eine Singularwertzerlegung durchfuhren, die in in jedem Fall
reelle Eigenvektoren erzeugt
M = UΓV T
U tU = I
V tV = I (2.72)
U (V ) nennt man die Matrix der linken (rechten) Eigenvektoren und die Diagonalmatrix
Γ die Matrix der singularen Werte. Ordnet man die Eigenvektoren aufsteigend nach dem
24
2 Multivariat normalverteilte Grundgesamtheiten
Betrag der Eigenwerte (also der vom Betrag her kleinste Eigenwert ist der erste), so findet
man typischerweise, dass die rechten Eigenvektoren dann großskalige Strukturen beschrei-
ben. Als Beispiel findet man in Abb.(2.6) vier ausgewahlte Eigenvektoren des Advektions-
Diffusionsmodells, wenn
• das Stromungsfeld ~v das Windfeld in 850 hPa im Wintermittel 1980-1999 ist und
• die Diffusionskonstante K = 5 × 106 m2sec−1 gross ist
Lasst man den Advektionsanteil gegen Null gehen und betrachtet den Diffusionanteil fur
die gesamte Kugel, erhalten wir eine physikalische Interpretation der Kugelfunktionen als
asymptotische , stationare Losungen der homogenen und isotropen Diffusionsgleichung.
2.7 Bedingte Wahrscheinlichkeitsdichten
Wir betrachten wiederum eine Normalverteilte Zufallsvariable ~X mit der pdf aus Gl.(2.13).
Dies wollen wir von nun an mit N(~µ, Σ) abkurzen. Diese ZVA sei in zwei Unterkomponenten
aufgespalten ~X == ( ~X1, ~X2). Entsprechen kann man den Mittelwert ~µ aufspalten in ~µ =
(~µ1, ~µ2). Die Kovarianzmatrix besteht dann aus drei verschiedenen Untermatrizen
Σ =
Σ11 Σ12
ΣT12 Σ22
(2.73)
Die Randverteilung fur die ZVA ~X1 ist dann N(~µ1, Σ11) und die fur ~X2 N (~µ2, Σ22). Damit
konnen wir nun die bedingten pdf’s
p( ~X1| ~X2 = ~x2) =p( ~X1, ~X2)
p( ~X2)(2.74)
bestimmen.
Da wir die Division einer Normalverteilung mit einer anderen Normalverteilung bestim-
men, ergibt sich fur die bedingte pdf wiederum eine Normalverteilung N (~µ1|2, Σ11|2) mit
~µ1|2 = ~µ1 + Σ12Σ−122 (~x2 − ~µ2)
Σ11|2 = Σ11 − Σ12Σ−122 ΣT
12 (2.75)
Der bedingte Erwartungswert fur ~X1 bei einem gegebenen ~X2 = ~x2 ist also eine lineare
Funktion des vorgegebenen Wertes ~x2. Dieser Wert ist ungleich dem unbedingtem Erwar-
tungswert ~µ1, wenn zwischen den beiden Variablen ~X1 und ~X2 eine Beziehung besteht, die
25
2 Multivariat normalverteilte Grundgesamtheiten
sich durch eine Kovarianzmatrix Σ12 6= 0 ausdruckt. Als ein Maß s fur den Zusammenhang
kann man das Verhaltnis der Varianz der bedingten ZVA ~X1|2 zur Varianz der unbedingten
ZVA ~X1 bestimmen.
s = 1 − Sp(Σ11|2)
Sp(Σ11)
=Sp(Σ12Σ
−122 ΣT
12)
Sp(Σ11)≤ 1 (2.76)
s ist Null, wenn kein Zusammenhang zwischen den beiden ZVA’s besteht und 1, wenn
Σ12 = Σ1211Σ
1222 ist.
26
2 Multivariat normalverteilte Grundgesamtheiten
Abbildung 4 Beispiel fur mogliche guess Muster zur Datenreduktion in einem periodischen Kanal
auf der Kugel zwischen 32.5◦S und 77.5◦N, rechte SVD Moden des Advektions-Diffusionsmodells
mit dem Stromungsfeld ~v in 850 hPa Wintermittel 1980-99 (links oben)
27
3 Bivariate Grundgesamtheiten
3 Bivariate Grundgesamtheiten
3.1 Eigenschaften der Grundgesamtheit
Als einfachsten Fall einer multivariaten Grundgesamtheit wollen wir eine bivariate, normal-
verteilte Grundgesamtheit betrachten.
~X = (X, Y ) = {(~x = (x, y)T , f(x, y)),S = IR2}
f(x, y) =1
2π det Σexp
(−1
2((x, y) − (µx, µy))
T Σ−1 ((x, y) − (µx, µy))
)
Σ =
σ2
x σxy
σxy σ2y
(3.1)
mit den Erwartungswerten ~µ = (µx, µy), den Varianzen σ2x, σ
2y und der Kovarianz σxy.
Betrachten wir die reduzierten ZVA
X∗ =X − µx
σx
Y ∗ =Y − µy
σy
(3.2)
so heißt die Kovarianz zwischen diesen reduzierte ZVA die Korrelation ρxy
ρxy = E(X∗Y ∗) (3.3)
Bilden wir die neue ZVA Z mit Hilfe der Linearkombination (t ∈ IR)
Z = tX∗ + Y ∗ (3.4)
und berechnen die Varianz dieser neuen ZVA
E(Z2) = t2V ar(X∗) + 2tρxy + V ar(Y ∗) = t2 + 2tρxy + 1
= (t + ρxy)2 + (1 − ρ2
xy) ≥ 0, (3.5)
so ergibt sich die Tatsache
ρ2xy ≤ 1
−1 ≤ ρxy ≤ 1 (3.6)
28
3 Bivariate Grundgesamtheiten
Eine andere Moglichkeit, zur Verknupfung der beiden Daten X und Y eine Aussage
zu machen, ist die Annahme einer funktionalen Beziehung zwischen den beiden Variablen
Y = g(X). Gemeint ist hier nicht eine Transformation (wie in der Vorlesung Statistik I be-
trachtet), da dann beide Variablen voneinander strikt abhangig waren, sondern ein Modell-
zusammenhang in dem Sinne, daß g(X) moglichst ahnlich zu Y ist, unter Berucksichtigung
eines bestimmmten Restterms R, also
Y = g(X) + R (3.7)
Eine allgemeine Losung fur ein beliebiges g ist nicht moglich, man muß eine gewisse Voraus-
wahl treffen. Beliebt sind die linearen Funktionen:
Y = αX + β + R (3.8)
mit den Unbekannten α, β. Die Anpassung einer linearen Funktion nennt man auch lineare
Regression. Mit der Annahme, daß das Residuum R eine Normalverteilte ZVA mit Erwar-
tungswert 0 und Varianz σ2r ist, wird die Bedingung ”moglichst ahnlich” durch die Forderung
E(R2) = σ2r
!= min (3.9)
realisiert. Einsetzen der Definition (3.8) ergibt
σ2r = E((Y − αX − β)2)
= E(Y 2) − 2αE(XY ) − 2βµy + α2E(X2) + β2 + 2αβµx (3.10)
Durch Ableitung nach den Unbekannten Variablen α, β erhalten wir dann das Gleichungs-
system
α(E(X2) − µ2x) = E(XY ) − µxµy
β = µy − αµx
α =σxy
σ2x
β = µy −σxy
σ2x
µx (3.11)
Ersetzen wir nun noch die Kovarianz durch die Korrelation, erhalten wir
α =σy
σx
ρxy
β = µy −σy
σx
ρxyµx (3.12)
29
3 Bivariate Grundgesamtheiten
D.h. der Korrelationskoeffizient ρxy ist ein Maß fur den linearen Zusammenhang zwischen
den beiden ZVA X und Y . Einsetzen der Losung in die Minimumsfunktion ergibt
σ2r = σ2
y −σ2
xy
σ2x
= σ2y(1 − ρ2
xy) (3.13)
Dann erkennen wir, daß
• X und Y vollstandig abhangig sind σ2r = 0, wenn ρxy = ±1 ist
• ρ2xy angibt, wieviel von der Varianz in Y durch die Varianz in X ausgedruckt werden
kann (erklarte oder dargestellte Varianz)
Nehmen wir nun noch die Ergebnisse der Rechnung zur bedingten pdf zweier Normalverteilter
ZVA (Gl.(2.75)), so erhalten wir fur den bivariaten Fall
Σ11 = σ2x
Σ22 = σ2y
Σ12 = σxy (3.14)
Damit bestimmt sich der bedingte Erwartungswert µ2|1 zu
µy|x = µy +σxy
σ2x
(x − µx) (3.15)
Dies ist die selbe Losung wie oben, hier jedoch als bedingter Erwartungswert berechnet. Of-
fensichtlich sind die least-squares Losung einer linearen Anpassung und der bedingte Erwar-
tungswert bei Annahme einer Normalverteilung identisch. Das Maß s zum Zusammenhang
zwischen x und y ist dann
s = 1 − Sp(Σ11|2)
Sp(Σ11)
=Sp(Σ12Σ
−122 ΣT
12)
Sp(Σ11)=
σ2xy
σ2xσ
2y
= ρ2xy (3.16)
Das Bestimmtheitsmaß ist als im bivariaten Fall identisch zum quadrierten Korrelationsko-
effizienten.
30
3 Bivariate Grundgesamtheiten
3.2 Schatzer der Eigenschaften der bivariaten GG
Gegeben sei eine Stichprobe der Lange m mit Realisierungen der ZVA X und Y :
x1, . . . , xm, y1, . . . , ym. Wie wir im letzten Kapitel gesehen haben, sind die ML Schatzer fur
Erwartungswert und die Elemente der Kovarianzmatrix gegeben durch
x =1
m
m∑
k=1
xk
y =1
m
m∑
k=1
yk
s2x =
1
m
m∑
k=1
(xk − x)2
s2y =
1
m
m∑
k=1
(yk − y)2
sxy =1
m
m∑
k=1
(xk − x)(yk − y) (3.17)
Als einen Schatzer rxy fur die Korrelation ρxy konnen wir dann schreiben
rxy =sxy
sxsy
(3.18)
Mit Hilfe der Dreiecksungleichung konnen wir sofort zeigen, daß gilt
s2xy ≤ s2
xs2y
r2xy ≤ 1
−1 ≤ rxy ≤ 1 (3.19)
Damit hat rxy ebenfalls die gleiche Eigenschaft wie die Korrelation der Grundgesamtheit.
Auch die Schatzer der Parameter der linearen Regression α und β konnen aus den ML
Schatzern fur Varianz und Kovarianz ermittelt werden.
α =sxy
s2x
β = y − αx (3.20)
Das identische Resultat erhalt man, wenn man α und β aus folgender Minimierungsaufgabe
bestimmt
yk = αxk + β + rk
s2r =
m∑
k=1
r2k
!= min (3.21)
31
3 Bivariate Grundgesamtheiten
Als Wert der Minimumsfunktion s2r(α, β) erhalten wir fur die Schatzer dann die analoge
Beziehung wie im Fall der Grundgesamtheit
s2r = s2
y(1 − r2xy), (3.22)
so daß rxy auch ein Maß fur den linearen Zusammenhang zwischen den beiden Stichproben
xk und yk ist und r2xy den Anteil an der Gesamtvarianz der y Stichprobe angibt, der durch
den linearen Zusammenhang zwischen x und y erklart oder dargestellt wird.
3.3 Konfidenzintervalle
Um ein Konfidenzintervall fur ρxy aus der Stichprobenkorrelation rxy zu bestimmen, mussen
wir bedenken, daß sowohl ρxy als auch rxy auf das Intervall ] − 1, 1[ beschrankt sind. Damit
kann die pdf der Stichprobenkorrelation nicht einer der Standardwahrscheinlichkeitsdichte-
funktionen gehorchen, da diese immer den Raum IR oder IR+ als Ereignisraum haben. Um
diese Problem zu umgehen, hat Fisher (1921) folgende Transformation vorgeschlagen (Fisher
- z Transformation)
z =1
2ln
1 + rxy
1 − rxy
(3.23)
und gezeigt, daß diese Transformation die ZVA Rxy (deren Realisierung die Stichprobenkor-
relation rxy ist) asymptotisch (d.h. m → ∞) in eine Normalverteile ZVA Z uberfuhrt mit
Erwartungswert und Varianz
E(Z) =1
2ln
1 + ρxy
1 − ρxy
V ar(Z) =1
m − 3(3.24)
Damit konnen wir wieder eine reduzierte NV ZVA bestimmen
Z∗ =√
m − 3(Z − E(Z)) (3.25)
und nach Intervallgrenzen ±c suchen, so daß bei a-priori gegebener Wahrscheinlichkeit γ gilt
Prob(−c ≤ Z∗ ≤ c) = 1 − γ
erf(c) =1 + γ
2(3.26)
Damit erhalten wir als Mutungsbereich fur E(Z)
z − c√m − 3
≤ E(Z) ≤ z +c√
m − 3(3.27)
32
3 Bivariate Grundgesamtheiten
Die inverse Operation zur Fisher z Transformation ist gegeben durch
rxy = tanh z =exp z − exp(−z)
exp z + exp(−z)(3.28)
Damit konnen wir ein Mutungsintervall fur ρ berechnen
tanh(z − c√m − 3
) ≤ ρxy ≤ tanh(z +c√
m − 3) (3.29)
Auch fur den Regressionparameter α lasst sich ein Konfidenzintervall angeben. Mit der
Annahme, daß die Residuen R NV mit Erwartungswert Null und Varianz σ2r sind, ist die
Große
t = sx
√(m − 1)(m − 2)
α − α√(m − 1)s2
r
(3.30)
die Realisierung einer Student - t verteilten ZVA mit (m − 2) Freiheitsgraden, da s2r eine
χ2 verteilte ZVA mit m − 2 Freiheitsgraden ist. Das Konfidenzintervall c fur diese Variable
bestimmt sich dann aus
FSt−t,m−2(c) =1 + γ
2(3.31)
und das Mutungsintervall fur α erhalt man dann aus
α − l ≤ α ≤ α + l
l = c
√(m − 1)s2
r√s2
x(m − 2)(m − 1)(3.32)
3.4 Hypothesentest
Ein beliebter Test im Rahmen der hier beschrieben bivariaten Statistik ist die Uberprufung
der Null bzw. Alternativhypothese
H0 : ρxy = 0
HA : ρxy 6= 0
oder HA : ρxy > 0 (3.33)
Bei Gultigkeit der Nullhypothese ist folgende Teststatistik
U = Rxy
√m − 2
1 − R2xy
(3.34)
33
3 Bivariate Grundgesamtheiten
eine Student - t verteilte ZVA mit (m − 2) Freiheitsgraden. Damit konnen wir fur den
zweiseitigen Test ρxy = 0 vs. ρxy 6= 0 Intervallgrenzen ±uα zum Irrtumsniveau α bestimmen,
so daß
Prob(−uα ≤ U ≤ uα|H0) = 1 − α (3.35)
Mit Hilfe der Verteilungsfunktion der Student - t Verteilung folgt dann
FSt−t,m−2(uα) = 1 − α
2(3.36)
Mit Hilfe der aus der Stichprobe bestimmten Realisierung der Testvariablen u
u = rxy
√m − 2
1 − r2xy
(3.37)
wird der Test nun durchgefuhrt
|u| ≤ uα akzeptiere H0
|u| > uα akzeptiere HA (3.38)
34
4 Multivariate Hypothesentests
4 Multivariate Hypothesentests
4.1 Einfuhrung
Die univariaten Hypothesentests auf unterschiedliche Erwartungswerte waren ausfuhrlich in
der Vorlesung Statistik beschrieben worden. War eine Stichproben (x1, . . . , xm) aus einer
Normalverteilten GG mit Erwartungswert µx und Varianz σ gegeben, so sollte uberpruft
werden, ob das Vorauswissen in Form einer Nullhypothese H0 akzeptiert werden kann oder
zugunsten der Alternativhypothes HA abgelehnt wird.
H0 : µx = µ0
HA : µx 6= µ0 (4.1)
Der Test wurde mittels einer sogenannten Teststatistik t durchgefuhrt, die mit Hilfe der
Schatzer von Erwartungswert und Varianz konstruiert wurde und deren pdf bzw. Wahr-
scheinlichkeitsfunktion bekannt war (die bekannte Student - t Verteilung)
t =√
mx − µ0
s
x =1
m
m∑
k=1
xk
s =1
m − 1
m∑
k=1
(xk − x)2 (4.2)
Die Testgrosse t ist im Fall einer gultigen Nullhypothese die Realisierung einer ZVA, die
zentral Student - t verteilt ist mit m − 1 Freiheitsgraden.
Diese Testvariable ist dimensionslos und invariant gegen lineare Transformationen. Schrei-
ben wir
y = ax + b (4.3)
mit a, b ∈ IR, erhalten wir
y = ax + b
µ0,y = aµ0 + b
s2y = a2s2 (4.4)
oder eingesetzt in die Definition fur t
t(y) =√
my − µ0y
sy
= t (4.5)
35
4 Multivariate Hypothesentests
den gesuchten Beweis der Invarianz gegen lineare Transformationen.
Fur den multivariaten Fall sei nun ebenfalls eine multivariate Stichprobe gegeben
~x1, . . . , ~xm, wobei ~xk ∈ IRq q dimensionale Vektorrealsierungen aus einer Normalverteilten
GG mit Erwartungswert ~µ und Kovarianzmatrix Σ sind. Auch das Vorauswissen verallge-
meinern wir zu einer multivariaten Nullhypothese bzw. Alternativhypothese
H0 :~µ = ~µ0
HA :~µ 6= ~µ0 (4.6)
Als Schatzer stehen zur Verfugung die ML Schatzer fur Erwartungswert und Kovarianzmatrix
~x =1
m
m∑
k=1
~xk
S =1
m − 1
m∑
k=1
(~xk − ~x)(~xk − ~x)T (4.7)
Damit ergibt sich die Frage, wie aus diesen Grossen eine brauchbare Testvariable konstru-
iert werden kann.
4.2 Die Hotelling T 2 Variable
Wir betrachten einen beliebigen, aber festen (von Null verschiedenen) Vektor ~a ∈ IRq. Das
Skalarprodukt mit der multivariaten ZVA ~x sei zx:
zx = ~aT~x (4.8)
Der Erwartungswert der univariaten ZVA zx ist dann wegen der Linearitat des Erwartungs-
werts
E(zx) = ~aT ~µx. (4.9)
Die Varianz Var(zx) bestimmt sich mit Hilfe der Rechenregeln aus dem 1. Kapitel zu
Var(zx) = E((zx − E(zx))2)
= E((~aT~x − ~aT ~µx)2)
= E(~aT (~x − ~µx)(~x − ~µx)T~a)
= ~aT Σ~a (4.10)
36
4 Multivariate Hypothesentests
Damit konnen wir aus den ML Schatzern eine Student - t ahnliche Variable bilden
t(~a) =√
m~aT ~x − ~aT ~µ0√
~aT S~a(4.11)
und eine Nullhypothese H0,a als Funktion des Vektors ~a akzeptieren, wenn wir finden
t2(~a) ≤ t2crit (4.12)
mit einem noch zu bestimmenden Quantil t2crit einer zunachst noch unbekannten pdf. Eine
multivariate Uberprufung der ursprunglichen Nullhypothese (4.6) muss nun die obigen Un-
gleichung fur alle ~a 6= 0 uberprufen, wobei wir aber noch berucksichtigen mussen, dass t
dimensionsfrei und invariant gegen Skalentransformationen ist. Wir konnen damit die Skala
fur ~a beliebig festlegen, z.B. durch
~aT S~a = 1 (4.13)
Die Akzeptanzregion fur die multivariate Nullhypothese ist dann die Vereinigungsmenge
aller ~a, die die Normierungsbedingung (4.13) und die Ungleichung (4.12) erfullen. Dann
konnen wir aber auch das Maximum aller t2(~a) Werte nehmen und fur die Akzeptanz der
multivariaten Nullhypothese fordern:
max~a
t2(~a) ≤ t2crit (4.14)
Die passende Teststatistik fur den multivariaten Test ist also
T 2 H0= max~a
m
(~aT ~x − ~aT ~µ0√
~aT S~a
)2
(4.15)
mit der Nebenbedingung (4.13). Dies ist eine Extremwertaufgabe fur den noch unbestimmten
Vektor ~a, die mit Hilfe der Lagrange’schen Multiplikatorenmethode gelost werden kann,
indem das Maximum folgender Funktion gesucht wird:
J = T 2(~a) + λ(1 − ~aT S~a)
T 2(~a) = m~aT (~x − ~µ0)(~x − ~µ0)T~a (4.16)
und λ bezeichnet den Lagrange’schen Multiplikator. Die notwendige Bedingung fur die Exis-
tenz eines Extremums lautet
~∇aJ =∂
∂~aJ = 0
∂J∂λ
= 0 (4.17)
37
4 Multivariate Hypothesentests
Mit Hilfe der im ersten Kapitel vorgestellten Rechenregeln zur Differentation von quadrati-
schen Formen erhalten wir fur den Gradienten
m((~x − ~µ0)(~x − ~µ0)T~a)T + m~aT (~x − ~µ0)(~x − ~µ0)
T − λ(S~a)T − λ~aT S = 0 (4.18)
wobei wir benutzt haben, dass S eine symmetrische Matrix ist S = ST . Die Ableitung
nach der Unbekannten λ ergibt naturlich die NB (4.13), so dass das Gleichungsystem zur
Bestimmung von ~a und λ lautet
(m(~x − ~µ0)(~x − ~µ0)
T − λS)~a = 0
~aT S~a = 1 (4.19)
Multiplikation der oberen Gleichung von links mit ~aT ergibt unter Benutzung der NB
λ = ~aT m(~x − ~µ0)(~x − ~µ0)T~a = T 2 > 0 (4.20)
Existiert die Inverse von S, d.h. S hat den vollen Rang q, konnen wir die obere Gleichung
von (4.19) mit S−1 von links multiplizieren und unter Berucksichtigung der NB schreiben
(S−1m(~x − ~µ0)(~x − ~µ0)
T − λI)~a = 0 (4.21)
Dies ist ein homogenes, lineares Gleichungssystem fur den unbekannten Vektor ~a. Eine nicht-
triviale Losung ~a 6= 0 existiert nur, wenn die Determinante der Systemmatrix verschwindet:
det(S−1m(~x − ~µ0)(~x − ~µ0)
T − λI)
= 0 (4.22)
Dies ist eine zweite Gleichung fur den Lagrange’schen Multiplikator λ und zeigt, dass λ eine
(positive) Nullstelle des charakteristischen Polynoms der Matrix(S−1m(~x − ~µ0)(~x − ~µ0)
T)
ist. Diese Matrix hat aber nur den Rang 1, da zwar S (und damit auch S−1) den vollen Rang
q hat, aber die Matrix (~x − ~µ0)(~x − ~µ0)T nur den Rang 1. Damit muss λ identisch zur Spur
der Systemmatrix sein, bzw. mit den Rechenregeln aus Kapitel 1
λ = Sp(S−1m(~x − ~µ0)(~x − ~µ0)
T)
= m(~x − ~µ0)S−1(~x − ~µ0) = T 2 (4.23)
Diese Testvariable ist auch unter dem Namen Hotelling T 2 Variable bekannt.
Da S und ~x die ML Schatzer fur Erwartungswert und Kovarianzmatrix sind, ist (m− 1)S
eine Wishart verteilte Zufallsmatrixvariable mit m− 1 Freiheitsgraden und Erwartungswert
38
4 Multivariate Hypothesentests
Σ. Gilt die Nullhypothese (4.6), so kann man nun die Verteilungsfunktion fur die Testvariable
angeben. Unter der Voraussetzung, dass S den vollen Rang q hat, ist die univariate ZVA
Z =m − q
q(m − 1)T 2 (4.24)
eine Fisher - F verteilte ZVA mit q, m − q Freiheitsgraden. Damit konnen wir nun den Test
durchfuhren, wenn wir a-priori eine Wahrscheinlichkeit fur die falschliche Ablehnung der
Nullhypothese (Irrtumswahrscheinlichkeit) α festlegen.
Prob(T 2 ≤ T 2crit) = 1 − α
FFisher−F,q,m−q(Zcrit) = 1 − α
T 2crit = Zcrit
q(m − 1)
m − q(4.25)
Bestimmen wir nun aus den Stichprobenwerten eine Realisierung der Hotelling Testvariable
(4.23), so konnen wir den Test (4.6) durchfuhren mit der Entscheidung
akzeptiere H0 wenn T 2 ≤ T 2crit
akzeptiere HA wenn T 2 > T 2crit (4.26)
Sei ~y eine q dimensionale, multivariat normal verteilte ZVA mit Erwartungswert µ und Ko-
varianzmatrix Σ, sowie mS eine Wishart verteilte Zufallsmatrix mit m Freiheitsgraden und
Erwartungswert Σ. Dann konnen wir eine verallgemeinerte Hotelling T 2 Variable betrachten:
T 2 = ~yTS−1~y (4.27)
Die univariate ZVA Z mit
Z =m − q + 1
mqT 2 (4.28)
ist eine Fisher F verteilte ZVA mit q, m−q+1 Freiheitsgraden und dem Nichtzentralitatspa-
rameter δ = ~µT Σ−1~µ. Mit ~µ = 0 und der entsprechenden Wahl von ~y = ~x sowie S und der
Anzahl der Freiheitsgrade erhalten wir sofort den obigen Test. Mit einer entsprechendem
Wahl eines Nichtzentralitatsparameters kann dann zusatzlich auch die Macht des Hotelling
T 2 Tests bestimmt werden.
Ebenso wie der Student - t Test invariant gegen lineare Transformationen der ZVA ist, ist
der Hotelling T 2 Test invariant gegen orthogonale Transformationen der Vektor ZVA. Sei G
die Matrix der Transformation. Dann ist eine orthogonale Transformation gegeben durch
GT G = GGT = I (4.29)
39
4 Multivariate Hypothesentests
Die Vektoren transformieren sich wie
~y = G~x
~x = GT~y (4.30)
Sei Sx der Schatzer der Kovarianzmatrix in der x Darstellung und Sy der entsprechende
Schatzer in der y Darstellung. Dann bestehen folgende Transformationsbeziehungen zwischen
den Matrizen und zwischen den Inversen
Sx = GT SyG
Sy = GSxGT
S−1x = GT S−1
y G
S−1y = GS−1
x GT (4.31)
Die Hotelling T 2 Testvariable in der x Darstellung lasst sich dann wie folgt umformen
T 2x = m(~x − ~µ0)
T S−1x (~x − ~µ0)
= m(GT ~y − ~µ0)T GT S−1
y G(GT ~y − ~µ0)
= m(~y − G~µ0)T GGT S−1
y GGT (~y − G~µ0)
= T 2y (4.32)
womit die Invarianz bewiesen ware.
4.3 Vergleich zweier multivariater Stichproben
Sei nun das Vorauswissen fur die Formulierung der Null/Alternativhypothese durch die Ent-
nahme einer anderen Stichprobe entstanden. D.h. wir haben zwei Stichproben
(~x1, . . . , ~xmx) ∈ N(~µx, Σ)
(~y1, . . . , ~ymy) ∈ N(~µy, Σ) (4.33)
40
4 Multivariate Hypothesentests
Erwartungswerte und Kovarianzmatrizen werden uber die ML Schatzung bestimmt.
~µxest= ~x
~µyest= ~y
Sx =1
mx − 1
mx∑
k=1
(~xk − ~x)(~xk − ~x)T =1
mx − 1Ax
Sy =1
my − 1
my∑
k=1
(~yk − ~y)(~yk − ~y)T =1
my − 1Ay (4.34)
Da beide Kovarianzmatrizen Schatzer fur die gemeinsame Kovarianzmatrix der GG sind,
konnen wir auch eine gemeinsam geschatzte Kovarianzmatrix S angeben (pooled estimate)
S =1
mx + my − 2(Ax + Ay) (4.35)
Dann ist S eine Wishart verteilte Zufallsmatrix mit mx + my − 2 Freiheitsgraden und dem
Erwartungswert Σ. Der Test der Null / Alternativhypothese
H0 : µx = µy
HA : µx 6= µy (4.36)
erfolgt nun uber die Hotelling T 2 Variable
T 2 =mxmy
mx + my
(~x − ~y)TS−1(~x − ~y) (4.37)
Gilt die Nullhypothese, ist die univariate ZVA Z
Z = T 2mx + my − q − 1
q(mx + my − 2)(4.38)
Fisher -F verteilt mit q, mx +my − 1− q Freiheitsgraden, so dass der Test vollig analog, wie
oben im Einstichprobenfall beschrieben, durchgefuhrt werden kann.
Wesentlich fur diesen Test ist die Annahme einer gemeinsamen Kovarianzmatrix Σ der
beiden zu betrachtenden GG. Untersuchen kann man nun den Einfluss auf den Ausgang des
Hotelling T 2 Tests, wenn diese Voraussetzung nicht erfullt ist. Sind die Stichprobenlangen
mx, my gross genug und gleich, gibt es keinen Effekt durch diese Annahme auf die Wahr-
scheinlichkeiten fur Fehlentscheidungen von Typ I und die Macht des Testes (s.u.). Wesentlich
komplizierter wird allerdings der Fall, wenn die Stichprobenlangen (wesentlich) differieren.
Es scheint angebracht, diesen Fall moglichst zu vermeiden (vergl. auch [9]).
41
4 Multivariate Hypothesentests
Es gilt in der GG Es gilt in der GG
Testausgang H0 wahr/H1 falsch H0 falsch/H1 wahr
H0 akzeptieren richtig, (1 − α) Fehler II β
H1 akzeptieren Fehler I, α richtig,1 − β
Tabelle 1 Entscheidungstabelle bei statistischen Hypothesentests, Definition Fehler I und
II Art
4.4 Die Macht des T 2 Tests
Im Fall der Gultigkeit der Alternativhypothese HA von (4.36) ist die ZVA
Z =mx + my − q − 1
q(mx + my − 2)
mxmy
mx + my
(~x − ~y)TS−1(~x − ~y) (4.39)
Nichtzentral Fisher-F verteilt mit (q, mx + my − q − 1) Freiheitsgraden und dem Nichtzen-
tralitatsparameter
δ2 =mxmy
mx + my
(~µx − ~µy)T Σ−1(~µx − ~µy) (4.40)
Wichtig ist hier, dass in den Nichtzentralitatsparameter nicht die Schatzer der Erwartungs-
werte bzw. der Kovarianzmatrizen, sondern die Werte der GG selbst eingehen, d.h. δ2 ist
keine ZVA. Dann lasst sich die Macht β des Hotelling T 2 Tests als Funktion des Nichtzen-
tralitatsparameters und der Entscheidungsvariablen Zcrit(α) berechnen zu
1 − β(δ2) = FFisher,NZ(δ2, Zcrit(α), q, mx + my − q − 1) (4.41)
Damit sind alle Eintrage in Tabelle (1) fur den Hotelling T 2 Test bestimmt, die Macht aller-
dings nur als Funktions eines beliebig vorgegebenen Nichtzentralitatsparameters, da zunachst
kein Schatzer fur δ2 existiert.
4.5 Informationskomprimierung und der T 2 Test
Die Existenz der Testvariablen T 2 erfordert eine nichtsingulare Stichprobenkovarianzmatrix
S, so dass S−1 berechnet werden kann. Dieser Fall ist gegeben, wenn der Stichprobenumfang
m − 1 bzw. mx + my − 2 grosser als die Dimension der betrachteten Vektoren q ist. Im
Standardfall fur klimatologische Probleme ist jedoch der umgekehrte Fall eher die Regel
m � q, so dass vor der Durchfuhrung eines Testes wieder eine Informationskomprimierung
42
4 Multivariate Hypothesentests
stattfinden muss. Dazu fuhren wir wieder guess Muster ~gl, l = 1, . . . , q, q < mx + my −1 ein, die a-priori festgelegt werden mussen und die selbst q dimensionale Vektoren sind
und z.B. typische Raummuster des zu untersuchenden Phanomens sind. Ordnen wir diese
q Spaltenvektoren wieder in einer q × q Matrix G an, so konnen wir die hochdimensionalen
Vektoren der Stichproben ~xk, k = 1, . . . , mx und ~yk, k = 1, . . . , my auch schreiben als
~xk = G~a(x)k + ~rx
~a(x)k = (GT G)−1GT~xk
~yk = G~a(y)k + ~ry
~a(y)k = (GT G)−1GT~yk (4.42)
Wir konnen nun (im Fall zweier Stichproben, analog ist der Fall einer Stichprobe zu behan-
deln) mit Hilfe der q dimensionalen Vektoren die folgende Null/Alternativhypothese uber-
prufen:
H0 : E(~a(x)) = E(~a(y))
HA : E(~a(x)) 6= E(~a(y)) (4.43)
Aufgrund der beschrankten Stichprobenlangen ist also eine fundierte statistische Aussage
nur bezuglich des a-priori ausgewahlten, q dimensionalen Unterraumes S~a ∈ IRq moglich.
Die Auswahl der guess-Muster ist also von entscheidender Bedeutung fur den Ausgang des
Tests. Wird die Auswahl ungeschickt getroffen, so ist die Nullhypothese zu akzeptieren,
obwohl in einem anderen Unterraum moglicherweise die Nullhypothese mit einer geringen
Irrtumswahrscheinlichkeit abgelehnt werden muss.
Aber auch im Fall einer bereits invertierbaren Kovarianzmatrix ist eine Informations-
komprimierung eine durchaus notwendige Massnahme, statistisch sinnvolle Ergebnisse zu
erhalten. Dazu wollen wir annehmen, dass die Kovarianzmatrix Σ bekannt ist und nicht
aus den Daten geschatzt werden muss. Dies erlaubt es uns weiterhin, die Daten ~xk (fur den
Fall einer Stichprobe und der Null/Alternativhypothese (4.6) in der Eigendarstellung der
43
4 Multivariate Hypothesentests
Kovarianzmatrix zu schreiben mit den Eigenvektoren ~el und Eigenwerten λl
~xk = E~bk
~µ = E ~β
ΣE = EΛ
Λ = diag(λ1, . . . , λq)
E = (~e1, . . . , ~eq) (4.44)
Den Erwartungswert ~µ schatzen wir durch den ML Schatzer ~x. Dieser ist fur den Fall, dass
die Stichprobe aus einer multivariaten Normalverteilung entnommen worden ist, ebenfalls
multivariat normalverteilt mit N(~µ, 1m
Σ). Damit ist der Wert der Likelihoodfunktion f(~x)
f(~x) =1√
2πqdet( 1m
Σ)exp(−m
2(~x − ~µ)TΣ−1(~x − ~µ)) (4.45)
In der Eigendarstellung der Kovarianzmatrix erhalten wir dann
f(~x) =
√mq
√2πqdet(Σ)
q∏
i=1
exp(−m
2
(bi − βi)2
λi
) (4.46)
und
det(Σ) =
q∏
i=1
λi (4.47)
Es gelte nun die Nullhypothese H0 : E(~x) = ~µ. Dann konnen wir nach der Wahrscheinlichkeit
fragen, mit der ein Differenzvektor ~x − ~µ langer als eine a-priori vorgebenen Lange R wird:
Prob((~x − ~µ)2 > R2|H0) = 1 −∫
|~x−~µ|<R
f(~x, Σ, ~µ)dq~x (4.48)
Mit Hilfe der Eigendarstellung konnen wir dann schreiben
Prob((~x−~µ)2 > R2|H0) = 1−∫
. . .
∫
V
√mq
√2πqdet(Σ)
q∏
i=1
exp(−m
2
(bi − βi)2
λi
)db1 . . . dbq (4.49)
Dabei wird das Integrationsvolumen V durch das Ellipsoid mit der Bedingung
~b ∈ V ↔q∑
i=1
(bi − βi)2
λi
≤ R2 (4.50)
bestimmt. Diese Wahrscheinlichkeit kann nach oben abgeschatzt werden, indem statt des
Ellipsoidvolumens der einhullenden Hyperkubus V ∗
~b ∈ V ∗ ↔ βi − Ri ≤ bi ≤ βi + Ri
Ri = a
√λi
m(4.51)
44
4 Multivariate Hypothesentests
als Integrationsvolumen genommen wird (mit a beliebig aber fest). Dann erhalten wir
Prob((~x − ~µ)2 > R2|H0) ≤ 1 −q∏
i=1
∫ βi+Ri
βi−Ri
√m√
2πλi
exp(−m
2
(bi − βi)2
λi
)dbi (4.52)
Fur alle Integrale in der letzten Gleichung konnen wir mit Hilfe der error - Funktion schreiben∫
(. . .)dbi = erf(a) − erf(−a) = (2erf(a) − 1) (4.53)
Damit erhalten wir als Abschatzung fur die gesuchte Wahrscheinlichkeit
Prob((~x − ~µ)2 > R2|H0) ≤
1 −q∏
i=1
m1q (2erf(a) − 1)
= 1 − (2erf(a) − 1)q (4.54)
Wahlen wir fur a = 1.96, so entspricht dies im eindimensionalen Fall q = 1 dem 95% Kon-
fidenzintervall bzw. die Wahrscheinlichkeit, dass bei Gultigkeit der Nullhypothese ein Ab-
stand des Mittelwerts vom Erwartungswerts anzutreffen ist, der mehr als doppelt so gross
ist wie die Standardabweichung λ1/m, liegt bei 5% und weniger. Betrachten wir nun aber 10
(unabhangige) Variablen q = 10, so ist die Wahrscheinlichkeit, dass der Schatzer des Erwar-
tungswertvektors ausserhalb des einhullenden Kubus des Ellipsoids liegt, bereits 0.4. D.h.
von 10 Realisierungen des Mittelwertes liegen etwa 4 ”weit” entfernt vom Erwartungswert.
Fur 100 Dimensionen ist diese Wahrscheinlichkeit schon praktisch 1 (genauer (1 - 0.006)). Bei
der Betrachtung vieler Freiheitsgrade besteht deshalb trotz Gultigkeit der Nullhypothese ei-
ne hohe Wahrscheinlichkeit, die Nullhypothese falschlich abzulehnen. Das Signifikanzniveau
nimmt exponentiell mit der Dimension ab. In einem hochdimensionalen Vektorraum ist prak-
tisch jeder Vektor groß. Dies nennt man auch etwas theatralisch den Fluch der Dimensionen.
Somit ist es wegen des Fluchs der Dimensionen auch bei einer nichtsingularen (Schatzung
der) Kovarianzmatrix sinnvoll, eine a-priori Informationskomprimierung durchzufuhren, um
q moglichst gering zu halten.
4.6 Beispiel zum Hotelling T 2 Test
Als ein Anwendungsbeispiel wollen wir und die Frage stellen, welche Struktur die Differenz
der bodennahen Lufttemperatur zwischen 30◦S und 75◦N im Fall eines pazifischen Warm-
Events (El Nino) relativ zu einem Cold-Event (La Nina) hat. Wir definieren somit zwei
Stichproben
45
4 Multivariate Hypothesentests
• die Warm-Falle
• die Kalt-Falle
und werden die Hypothese testen
H0 : E(T2m|Kalt) = E(T2m|Warm)
HA : E(T2m|Kalt) 6= E(T2m|Warm) (4.55)
In Tabelle (4.6) findet man die Jahre (das sogenannte Jahr 0 fur jedes der Ereignisse,
das Maximum eines Warm/Kalt Ereignisses im Pazifik ist jeweils zum Ende des Jahres 0 /
Beginn des Jahres+1), die die Stichproben definieren.
Fur diese Jahre werden die Stichprobenmittelwerte und Stichprobenkovarianzmatrizen
benotigt, um den Test durchzufuhren. Aus dem Stichprobenumfang berechnet sich die An-
zahl der Freiheitsgrade der gemeinsamen Kovarianzmatrix zu 25 + 23− 2 = 46. Damit kann
man also maximal 46 Gitterpunkte in der Ortsraumdarstellung multivariat untersuchen. We-
sentlich angebrachter ist deswegen eine Datenreduktion mit Hilfe von guess-Mustern. Diese
sind im nachfolgenden Beispiel die rechten SVD Moden des Advektions-Diffusionsmodells, so
wie sie in wenigen Beispielen in Abb.(2.6) dargestellt sind. Damit wird gleichzeitig zur statis-
tischen Hypothese auch eine physikalische Erklarung untersucht, namlich ob die gewahlten
Muster, die einem vereinfachten physikalischen Modell entnommen wurden, einen signifikan-
ten Unterschied zwischen den beiden Stichproben aufspannen konnen.
Die Ergebnisse des Tests fur gleitende 3-Monatsmittel jeweils zentriert auf den mittleren
Monat sind in den Abb.(4.6) bis (4.6) dargestellt. In Abb.(4.6) findet man in der oberen
Abbildung das Verhaltnis der Hotelling T 2 Variable zum 5% Quantil der Fisher-F Vertei-
lung dargestellt, wenn die Anzahl der moglichen Modenamplituden von 1 auf 44 erhoht
wird (Ordinate). Eine solche Form von sequentiellen Tests nennt man eine Testhierarchie.
Allerdings sind diese sequentiellen Tests nicht unabhangig voneinander, sie zeigen jedoch,
in welchen Modenamplituden ein Signal verborgen sein kann. Auf der Abzisse ist der Zen-
tralmonat des jeweiligen 3-Monatsfensters dargestellt. Die 1 bezeichnet also die Mittelwerte
Dezember Jahr-1 uber Januar Jahr-0 bis Februar Jahr-0 und die 13 den Zeitraum Dezember
Jahr-0 bis Februar Jahr+1. Man erkennt klar, dass das Signal im April/Mai einsetzt, seine
starkste Auspragung im borealen Spatsommer /Fruhherbst annimmt und zum Jahr+1 hin
46
4 Multivariate Hypothesentests
Warm-Events Cold Events
1884 1886
1887 1889
1891 1892
1896 1898
1899 1903
1902 1906
1905 1908
1911 1916
1914 1920
1918 1924
1923 1931
1925 1938
1930 1942
1932 1949
1939 1954
1941 1964
1951 1966
1953 1970
1957 1973
1965 1975
1969 1978
1972 1988
1976
1982
1986
Tabelle 2 Jahre der pazifischen Warm- und Kalt Ereignisse
47
4 Multivariate Hypothesentests
allmahlich ausklingt. In Abb.(4.6) ist das Gesamtsignal (Differenz aller 250 Modenamplitu-
den, die zur Entwicklung des Feldes ausgewahlt waren) in der oberen Abbildung dargestellt,
wahrend unten der Anteil des Feldes prasentiert wird, der die hochste T 2 Variable in dem
3-Monatszeitraum OND aufweist. Ein Blick in Abb.(4.6) zeigt, dass dies die Testhierarchien
mit ca. 35-37 Modenamplituden sind, die etwa 64 % der raumlichen Varianz des gesamten
Feldes aus der oberen Abbildung erfassen. Diese Differenz kann man als signifikante Differenz
bezeichnen (obwohl sie es streng genommen nicht ist).
48
4 Multivariate Hypothesentests
Abbildung 5 Ergebnisse des Hypothesentests (Gl.(4.55)) mit Hilfe des Hotelling T 2 Tests, in der
oberen Abbildung ist das Verhaltnis der Testvariablen zum 5% Quantil der Fisher F Verteilung
dargestellt fur alle Werte uber 1, hier wird die Nullhypothese mindestens zum Signifikanzniveau
5% abgelehnt; in der unteren Abbildung ist das Rekurrenzniveau in % (siehe Kapitel (5)) bestimmt
nach der OS-Methode bei einer Ablehnung der Nullhypothese zum 5% Signifikanzniveau dargestellt
49
4 Multivariate Hypothesentests
Abbildung 6 Feld der mittlere Differenz aller pazifischen Warm- und Kaltfalle in den letzten 120
Jahren im 3-Monatszeitraum OND (oben), Feld der signifikanten Differenz aller pazifischen Warm-
und Kaltfalle in den letzten 120 Jahren im 3-Monatszeitraum OND (unten) bestimmt aus den
ersten 35 Moden, die zum 5% signifikant von Null verschieden sind
50
4 Multivariate Hypothesentests
Abbildung 7 Wie Abb.(4.6) jedoch fur den 3-Monatszeitraum Dezember Jahr-0 bis Februar
Jahr+1
51
5 Diskriminanzanalyse
5 Diskriminanzanalyse
5.1 Das Problem
Gegeben seien zwei Stichproben (~x1, . . . , ~xmx) und ~y1, . . . , ~ymy
) sowie ein weiteres Stichpro-
benelement ~z. I.A. kann diese Einzelbeobachtung aufgrund ihrer individuellen Struktur nicht
so ohne weiteres als Element einer der beiden Stichproben klassifiziert werden. Vielmehr kann
nur der Vergleich mit den Stichproben diese Frage beantworten. Fur den Fall, dass die bei-
den Stichproben zwei unterschiedlichen Grundgesamtheiten entnommen worden sind, kann
durch diese Entscheidung festgestellt werden, aus welcher GG die Einzelbeobachtung wahr-
scheinlich entnommen worden ist. Bei zwei GG ist es nun das Ziel, den Stichprobenraum,
uber dem beide GG definiert sind, in zwei disjunkte Untervolumina Rx und Ry aufzuteilen,
so dass die Entscheidung ~z ∈ GGx durch die Entscheidung ~z ∈ Rx oder entsprechend fur y
ersetzt werden kann.
Dann kann man folgenden Tabelle aufstellen:
Realitat Entscheidung
~z ∈ Rx ~z ∈ Ry
~z ∈ GGx o.k. C(y|x)
~z ∈ GGy C(x|y) o.k.
Dabei sind C(x|y) bzw. C(y|x) die Kosten fur die Fehlklassifikation. Ein gutes Klassifika-
tionschema liegt also vor, wenn diese Kosten moglichst gering ausfallen. Sei fx(~s) bzw. fy(~s)
die pdf fur die jeweilige GG. Die Wahrscheinlichkeit, eine korrekte Entscheidung zu treffen
bzgl. der Einordnung in die x GG lautet dann:
P (x|x) =
∫
Rx
fx(~s)dqs (5.1)
Die Wahrscheinlichkeit fur eine Fehlklassifikation
~z ∈ Ry obwohl ~z ∈ GGx (5.2)
bestimmt sich dann aus dem Integral
P (y|x) =
∫
Ry
fx(~s)dqs (5.3)
52
5 Diskriminanzanalyse
Analog gilt fur die anderen Entscheidungen
P (y|y) =
∫
Ry
fy(~s)dqs
P (x|y) =
∫
Rx
fy(~s)dqs (5.4)
Sei ferner noch Qx die Wahrscheinlichkeit, dass eine Beobachtung aus der x GG entnommen
wird, und Qy die entsprechende Zahl fur die y GG. Dann gilt Qx +Qy = 1, weil nur zwei GG
zugelassen sind. Wir konnen nun die mittleren Kosten fur die Fehlklassifikationen angeben
C = C(x|y)P (x|y)Qy + C(y|x)P (y|x)Qx (5.5)
Ein Entscheidungsverfahren, dass nun diese mittleren Kosten C minimal halt, heisst ein
Bayes - Entscheidungsverfahren. Die Losung lautet, dass Rx und Ry die Teilmengen des
Stichprobenraums sind, fur die gilt
Rx =⋃
~s
: [C(y|x)Qx]fx(~s) ≥ [C(x|y)Qy]fy(~s)
Ry =⋃
~s
: [C(y|x)Qy]fx(~s) < [C(x|y)Qx]fy(~s) (5.6)
Dies lasst sich folgendermassen beweisen. Fur zwei beliebige Klassifikationsregionen R∗x, R
∗y
sind die mittleren Kosten fur die Fehlklassifikation
C = c(y|x)Qx
∫
R∗
y
fx(~s)dq~s + c(x|y)Qy
∫
R∗
x
fy(~s)dq~s
=
∫
R∗
y
c(y|x)Qxfx(~s)dq~s −
∫
R∗
y
c(x|y)Qyfy(~s)dq~s +
∫
R∗
y
c(x|y)Qyfy(~s)dq~s +
∫
R∗
x
c(x|y)Qyfy(~s)dq~s
=
∫
R∗
y
(c(y|x)Qxfx(~s) − c(x|y)Qyfy(~s))dq~s +
∫
R∗
y∪R∗
x
c(x|y)Qyfy(~s)dq~s
︸ ︷︷ ︸=c(x|y)Qy
(5.7)
da sowohl c, Qxy als auch f positiv definit sind, wird C nur dann minimal, wenn in der Region
R∗y nur die Punkte ~s liegen, fur die gilt
c(y|x)Qxfx(~s) − C(x|y)Qyfy(~s) < 0 (5.8)
und die Punkte ausschliesst, fur die gilt
c(y|x)Qxfx(~s) − C(x|y)Qyfy(~s) ≥ 0 (5.9)
53
5 Diskriminanzanalyse
Schliessen wir die Menge aller einzelnen Punkte aus, fur die gilt
fx(~s)
fy(~s)=
Qyc(x|y)
Qxc(y|x)= 0, (5.10)
ist die Entscheidungsregel
~z ∈ GGx wenn ~z ∈ Rx
~z ∈ GGy wenn ~z ∈ Ry
Rx =⋃
~s
mitfx(~s)
fy(~s)≥ Qyc(x|y)
Qxc(y|x)
Ry =⋃
~s
mitfx(~s)
fy(~s)<
Qyc(x|y)
Qxc(y|x)(5.11)
eine Bayes Entscheidungsmethode bis auf die Punktemenge, fur die (5.10) gilt.
5.2 Klassifikation bei multivariaten Normalverteilungen
Wir betrachten nun zwei Normalverteilte GG mit gemeinsamer Kovarianzmatrix und unter-
schiedlichem Erwartungswert ~µx, ~µy (s.o. T 2 Test). Dann haben die im letzten Unterkapitel
angegebenen pdf’s die Form
fx(~s) =1√
2πqdetΣexp(−1
2(~s − ~µx)
T Σ−1(~s − ~µx))
fy(~s) =1√
2πqdetΣexp(−1
2(~s − ~µy)
T Σ−1(~s − ~µy)) (5.12)
Sei nun das Verhaltnis der Kosten fur die Fehlklassifikation bekannt
c(x|y)Qy
c(y|x)Qx
= k (5.13)
wobei der einfachste Fall c(x|y) = c(y|x) und Qx = Qy mit k = 1 ist. Dann ist die Bayes
Entscheidung gegeben durch die Regionen
Rx =⋃
~s
mit
fx(~s)
fy(~s)= exp
(−1
2(~s − ~µx)
T Σ−1(~s − ~µx) +1
2(~s − ~µy)
T Σ−1(~s − ~µy)
)≥ k
Ry =⋃
~s
mit
fx(~s)
fy(~s)= exp
(−1
2(~s − ~µx)
T Σ−1(~s − ~µx) +1
2(~s − ~µy)
T Σ−1(~s − ~µy)
)< k (5.14)
54
5 Diskriminanzanalyse
oder mit Hilfe der Logarithmus Funktion (streng monoton steigend)
Rx =⋃
~s
mit − 1
2
((~s − ~µx)
T Σ−1(~s − ~µx) − (~s − ~µy)T Σ−1(~s − ~µy)
)≥ ln k
Ry =⋃
~s
mit − 1
2
((~s − ~µx)
T Σ−1(~s − ~µx) − (~s − ~µy)T Σ−1(~s − ~µy)
)< ln k (5.15)
Durch ein einfache Umformung der linken Seite der Ungleichung erhalten wir als aquivalente
Entscheidungsmethode
Rx =⋃
~s
mit U(~s) ≥ ln k
Ry =⋃
~s
mit U(~s) < ln k
U(~s) = ~sT Σ−1(~µx − ~µy) − 1
2(~µx + ~µy)Σ
−1(~µx − ~µy) (5.16)
Damit erhalten wir folgende geometrische Interpretation der Entscheidungsregeln. Die Glei-
chung
U(~s) − ln k = 0 (5.17)
beschreibt eine q − 1 dimensionale Hyperebene (lineare Flache) im IRq. Liegt die zu klassifi-
zierende Beobachtung ~z auf der Seite der Flache, fur die gilt U(~z) ≥ ln k, wird ~z der GG x
zugeordnet, im Fall U(~z) < ln k der GG y (vergl. Abb. (8) fur den ein- und zweidimensionalen
Fall).
Wenn ~z eine Realisierung aus einer multivariaten NV N (~µx, Σ) oder N (~µy, Σ) ist, ist
die ZVA U eine lineare Kombination von Normalverteilten ZVA und damit wiederum eine
Normalverteilte ZVA, von der der Erwartungswert E(U) und die Varianz V ar(U) bestimmt
werden konnen. Ist ~z ∈ N (~µx, Σ), ist
E(U) = E(~zT )Σ−1(~µx − ~µy) −1
2(~µx + ~µy)Σ
−1(~µx − ~µy)
= ~µTx Σ−1(~µx − ~µy) −
1
2(~µx + ~µy)Σ
−1(~µx − ~µy)
=1
2(~µx − ~µy)
T Σ−1(~µx − ~µy) =1
2∆2 (5.18)
wobei ∆ der Mahalanobisabstand ist. Ist dagegen ~z ∈ N (~µy, Σ), so ist
E(U) = E(~zT )Σ−1(~µx − ~µy) −1
2(~µx + ~µy)Σ
−1(~µx − ~µy)
= ~µTy Σ−1(~µx − ~µy) −
1
2(~µx + ~µy)Σ
−1(~µx − ~µy)
= −1
2(~µx − ~µy)
T Σ−1(~µx − ~µy) = −1
2∆2 (5.19)
55
5 Diskriminanzanalyse
Die Varianz von U bestimmt sich aus V ar(U) = E((U − E(U))2), wenn ~z ∈ N (~µx, Σ) zu
U − E(U) = ~zT Σ−1(~µx − ~µy) −1
2(~µx + ~µy)Σ
−1(~µx − ~µy) −1
2∆2
= (~z − ~µx)TΣ−1(~µx − ~µy)
V ar(U) = E((~µx − ~µy)T Σ−1(~z − ~µx)(~z − ~µx)
TΣ−1(~µx − ~µy)
= (~µx − ~µy)T Σ−1 E((~z − ~µx)(~z − ~µx)
T )︸ ︷︷ ︸=Σ
Σ−1(~µx − ~µy)
= (~µx − ~µy)TΣ−1(~µx − ~µy)
= ∆2 (5.20)
Das gleiche Ergebnis erhalten wir fur den Fall ~z ∈ N (~µx, Σ), so dass gilt
~z ∈ N (~µx, Σ) → U ∈ N (+1
2∆2, ∆2)
~z ∈ N (~µy, Σ) → U ∈ N (−1
2∆2, ∆2) (5.21)
Damit lassen sich die Wahrscheinlichkeiten fur Fehlklassifikationen sofort berechnen
P (y|x) =
∫ lnk
−∞
fx(u)du
P (x|y) =
∫ ∞
ln k
fy(u)du (5.22)
Fur den Fall k = 1 ist dann
P (y|x) = P (x|y) = erf(∆
2) (5.23)
Abb (9) zeigt nochmals die entsprechende geometrische Interpretation.
5.3 Klassifikation bei unbekannter pdf
Seien nun (~x1, . . . , ~xmx) und (~y1, . . . , ~ymy
) Stichproben aus zwei multivariaten Normalvertei-
lungen mit Erwartungswerten ~µx, ~µy und der gemeinsamen Kovarianzmatrix Σ. Diese drei
Parameter sind naturlich unbekannt und mussen aus den Stichproben geschatzt werden z.B.
mit Hilfe der ML Schatzer:
~x =1
m
mx∑
k=1
~xk
~y =1
m
my∑
k=1
~yk
(mx + my − 2)S =mx∑
k=1
(~xk − ~x)(~xk − ~xT ) +
my∑
k=1
(~yk − ~y)(~yk − ~yT ) (5.24)
56
5 Diskriminanzanalyse
Einen Schatzer fur die Bayes U Entscheidungsstatistik erhalten wir, wenn wir die Werte der
Grundgesamtheit durch ihre Schatzer ersetzen:
U(~z) = W = ~zT S−1(~x − ~y) − 1
2(~x + ~y)TS−1(~x − ~y) (5.25)
Die W Statistik ist auch unter dem Namen ”Wald - Anderson” Statistik bekannt. Die Ent-
scheidung, zu welcher GG der Stichprobenvektor ~z gehort, erfolgt dann gemass
W (~z) ≥ 0 ⇒ ~z ∈ N (~µx, Σ)
W (~z) < 0 ⇒ ~z ∈ N (~µy, Σ) (5.26)
Um die Wahrscheinlichkeiten fur eine Fehlklassifikation zu bestimmen, haben wir im Fall
bekannter Parameter Σ, ~µx,y die pdf der Statistik U zu benutzen. Im Fall der geschatzten Pa-
rameter S~x, ~y ist aber die pdf der Wald - Anderson Statistik W ausserordentlich kompliziert
und z.B. eine Funktion der Stichprobenlange sowie der unbekannten Mahalanobisdistanz ∆.
Weiterhin ist i.A.
fW (w, z ∈ N (µx, Σ), mx, my, ∆) 6= fW (w, z ∈ N (µy, Σ), mx, my, ∆) (5.27)
ausser im Fall mx = my. Es lasst sich aber zeigen ([10], S.211/212), dass fur den Fall
mx, my → ∞ die Wahrscheinlichkeitsdichte fur W gegen die pdf von U (die Normalverteilung
N (12∆2, ∆2) ) konvergiert. Mit Hilfe einer Reihenentwicklung in ε ∼ 1
mx, 1
my(Okamoto, 1963,
[11]) kann man auch noch die Abweichungen in erster Ordnung fW (w) von fU bestimmen
und damit die Fehlklassifikationwahrscheinlichkeiten berechnen:
Prob
(W (~z) − 1
2∆2
∆≤ u|~z ∈ N (~µx, Σ)
)= erf(u) −
φ(u)(1
2mx∆2(u3 + (q − 3)u − q∆)
+1
2my∆2(u3 + 2∆u + (q − 3 + ∆2)u + (p − 2)∆)
+1
m(4u3 + 4∆u2 + (6q − 6 + ∆2)u + 2(p − 1)∆)
φ(u) =1√2π
exp(−1
2u2)
m = mx + my − 2 (5.28)
Die korrespondierende Beziehung fur die andere Fehlklassifikation Prob(−W (~z)+ 12∆2
∆≥ u|~z ∈
N (~µy, Σ)) erhalten wir aus Gl. (5.28) durch Vertauschen von mx und my. Die oben ange-
gebenen Wahrscheinlichkeiten sind nun noch immer unbekannt, da sie von der unbekannten
57
5 Diskriminanzanalyse
Mahalanobisdistanz ∆ abhangen. Einen Schatzer fur die Wahrscheinlichkeiten der Fehlklas-
sifikation erhalten wir durch die Ersetzung von ∆ durch den Schatzer D mit
D2 = (~x − ~y)TS−1(~x − ~y) (5.29)
Da aber S−1 ein verzerrter Schatzer fur Σ−1 ist, ist auch D2 einen verzerrte Schatzung fur
∆2. Der Erwartungswert fur D2 bestimmt sich aus der Wishartverteilung zu
E(D2) =mx + my − 2
mx + my − q − 3(∆2 + q
mx + my
mxmy
) (5.30)
Damit ist es moglich, die folgende Korrektur zu benutzen (”shrunken” D, DS)
DS2 = D2mx + my − q − 3
mx + my − 2(5.31)
Die Verschiebung −q mx+my
mxmywird i.A. nicht benutzt, da dann nicht mehr garantiert werden
kann, dass DS2 positiv definit bleibt, wie es fur ein Abstandsmaß sein muss.
5.4 Beispiel zur Anwendung der Diskriminanzanalyse
In der meteorologisch - klimatologischen Gemeinde ist statt des Begriffs ”Wahrscheinlichkeit
der Fehlklassifikation” der Begriff ”Rekurrenz” gepragt worden (v. Storch und Zwiers (1989)
[13] bzw. Zwiers und v. Storch (1990) [14]). Die Beziehung zwischen diesen beiden Beziehun-
gen ist jedoch sehr einfach, da man unter Rekurrenz die Wahrscheinlichkeit der korrekten
Klassifikation versteht. Dann gilt Rekurrenz = 1 - Fehlklassifikationswahrscheinlichkeit. Mit
Hilfe der Rekurrenz kann man die Aussagen der statistischen Mittelwerttests (z.B. den Ho-
telling T 2 oder den Student-t Test) erheblich verbessern. Das Problem der Mittelwerttests
ist, dass sie umso scharfer werden,je grossere Stichproben vorliegen. Die Hotelling Testva-
riable nimmt dabei typischerweise mit N zu, die Student -t Variable mit√
N . Damit kann
es passieren, dass bei sehr umfangreichen Stichproben, signifikante Unterschiede detektiert
werden, die von ihrer absoluten Grosse her sehr klein sind. Kritisch wird dieses Problem z.B.
beim einer Modellvalidation: Hier testet man ob der Erwartungswert einer Variablen, die
uni- oder multivariate sein mag, des Modell identisch ist zum Erwartungswert in den Beob-
achtungen. Bei einem hinreichend großem Stichprobenumfang findet man dann jeden noch
so winzigen Modellfehler und muss die Nullhypothese der identischen Erwartungswerte und
damit die Ubereinstimmung von Modell und Realitat ablehnen. Viel geeigneter und besser
58
5 Diskriminanzanalyse
ist jedoch die Frage, ob ein Stichprobenelement,das entweder dem Modell oder der Realitat
entnommen sein kann, mit einer zu bestimmenden Fehlklassifikationswahrscheinlichkeit dem
korrekten Ensemble zugeordnet werden kann. Hierbei werden die Stichprobenelemente direkt
untereinander verglichen und mit ihrer individuellen Standardabweichung (Kovarianzmatrix)
gewichtet. Diese Aussagen sind damit sehr viel realitatsnaher als die eigentlichen Mittelwert-
tests. In der sogenannten Bayesischen Statistik wird diese spezielle Denkweise noch detailliert
ausgefuhrt ([4]).
Die Abb.(4.6) zeigt in ihrem unteren Teil die berechnete Wahrscheinlichkeit fur eine kor-
rekte Klassifikation (”Rekurrenz”) fur den Hypothesentest Gl.(4.55). Angewandt wurde die
Okamoto-Methode (Gl.(5.28)) unter Verwendung des geschrumpften Schatzers fur die Maha-
lanobisdistanz DS2 (Gl.(5.31)). Man erkennt, dass das Differenzsignal im borealen Spatsom-
mer/Fruhherbst und dann wieder im Dezember Jahr-0 bis Februar Jahr+1 einen sehr gute
Zuordnung der individuellen 3-Monatsmittel zu einem der beiden Ensemble Warm oder Kalt
zulasst, wenn man es mit 35-37 Amplituden der rSVD Moden beschreibt. Die Fehlklassifi-
kationswahrscheinlichkeiten (=1-Rekurrenz) liegen z.T. deutlich unter 10%.
59
5 Diskriminanzanalyse
-6 -4 -2 0 2 4 60
0.05
0.1
0.15
0.2
0.25
Eindimensionale Bayes Entscheidung, k=2
s
U(s)
-5 -4 -3 -2 -1 0 1 2 3 4 5-5
-4
-3
-2
-1
0
1
2
3
4
5Bayes Entscheidung 2-dim, k=2
s1
s2
U(s1,s2)-log(k)=0
Abbildung 8 Zur geometrischen Interpretation der Bayes Entscheidungsregel bei Normalverteilten
GG mit identischer Kovarianzmatrix im eindimensionalen (a) und zweidimensionalen (b) Fall
60
5 Diskriminanzanalyse
-6 -4 -2 0 2 4 60
0.05
0.1
0.15
0.2
0.25
Wahrscheinlichkeiten der Fehlklassifikation, k=2
u
f(u)
U(s)
P(y;x) P(x;y)
f(u,z in GGy) f(u,z in GGx)
Abbildung 9 Zur Definition der Fehlklassifikationen mit Hilfe der Bayes Entscheidungsstatistik U
61
6 Empirische Orthogonalfunktionen
6 Empirische Orthogonalfunktionen
6.1 Definition
Nachdem wir im vorletzten Kapitel gesehen haben, dass eine wesentliche Eigenart der multi-
variaten Statistik die Notwendigkeit einer a-priori Datenreduktion ist, erhebt sich die Frage,
ob es eine besonders effektive Art der Datenkomprimierung gibt. Ausserdem waren eventuell
guess - Muster, die orthogonal zueinander stehen und die Lange 1 haben, aus methodischen
Grunden (die Matrix GtG ist dann die Einheitsmatrix) interessant. Um dieses rauszufinden,
muss man/frau erst einmal definieren, was unter ”besonders effektiv” zu verstehen ist.
Sei ~X eine q-dimensionale multivariate ZVA und Σ = E(XX t) die Matrix der zweiten Mo-
mente (nicht notwendigerweise zentriert). Seien λ1 ≥ λ2 ≥ . . . ≥ λq die Eigenwerte und
~u1, . . . , ~uq die dazugehorigen Eigenvektoren der Matrix Σ mit |~uj| = 1. Dann gilt fur alle k
mit 1 ≤ k ≤ q und fur jeden andere Orthonormalbasis ~v1, . . . , ~vq
E(|( ~X −k∑
j=1
( ~X t~uj)~uj)|) ≤ E(|( ~X −k∑
j=1
( ~X t~vj)~vj)|) (6.1)
wobei
( ~X t~uj) = aj( ~X) (6.2)
das Skalarprodukt zwischen den Vektoren ~X und ~uj bzw. analog ~vj bedeutet. Ausserdem
gilt:
E(|( ~X −k∑
j=1
aj( ~X)~uj)|2) = E(| ~X|2) −k∑
j=1
λj =
q∑
j=k+1
λj (6.3)
und falls E( ~X) = 0
E(aj( ~X) = 0
E(aiaj) = λjδi,j (6.4)
Der Beweis wird im Abschnitt (6.5) dieses Kapitels geliefert. Dies bedeutet, dass die
Zerlegung eines Vektors ~X durch die Eigenvektoren der Matrix der zweiten Momente
die (im quadratischen Mittel) effektivste Darstellung liefert. Die Entwicklung nach einem
anderen ONB System liefert mittlere quadratische Abweichungen, die immer grosser sind
als die mittleren quadratischen Abweichungen bei der Zerlegung nach den Eigenvektoren.
Dementsprechend kann man/frau einen hochdimensionalen Vektor ~X u.U. durch sehr viel
62
6 Empirische Orthogonalfunktionen
weniger Entwicklungskoeffizienten aj( ~X) darstellen, ohne viel Restvarianz vernachlassigen
zu mussen. Der nichterklarte mittlere quadratische Rest ist dann gleich der Summe der
Eigenwerte der nicht benutzten Eigenvektoren. Ist die Matrix Σ die Matrix der zentrierten
zweiten Momente, also die Kovarianzmatrix und ist der Erwartungswert der Vektor ZVA
~X = 0, so sind die Erwartungswerte der Entwicklungskoeffizienten ebenfalls null und die
Entwicklungskoeffizienten sind statistisch unabhangig (siehe Gl. (6.4)). Weil also das Eigen-
vektorsystem der Matrix der zweiten Momente derart einzigartige Eigenschaften aufweist,
nennt man/frau die Eigenvektoren auch die Prinzipalen Vektoren, die Entwicklungsko-
effizienten aj( ~X) heißen Prinzipale Komponenten. Hat man/frau eine Schatzung der
Kovarianzmatrix (bzw. der Matrix der zweiten Momente), so kann diese selbstverstandlich
diagonalisiert werden. Die resultierenden Eigenvektoren sind dann Schatzer der Prinzipalen
Vektoren und heiaen empirische Orthogonalfunktionen (EOF), die Schatzer der
Prinzipalen Komponenten nennt man die Amplituden der EOF’s. Diese sind im Bereich
der Meteorologie zuerst von E.N. Lorenz (1956) ([12]) beschrieben worden und erfreuen
sich seither steigender Beliebtheit. Ursprunglichen zweidimensionale Felddaten (etwa die
bodennahen Lufttemperaturen oder der Bodendruck im Nordatlantik) waren ja zu einem q
(q ' 160 beim Bodendruck im Nordatlantik bei einer 5◦ × 10◦ Auflosung oder q ' 1660 bei
den 2m Temperaturen zwischen 30◦S und 75◦N mit einer 5◦ × 5◦ Auflosung) -dimensionalen
Datenvektor umgeordnet worden. Entsprechen konnen jetzt auch die Eigenvektoren wieder
aus der Vektorform zu einem zweidimensionalen Feld zuruckgeordnet werden. Diese zeigt,
dass die prinzipalen Vektoren typische raumliche Muster darstellen, die sozusagen immer
wieder auftreten.
6.2 Schatzung der Prinzipalen Vektoren
Sei D die Stichprobenmatrix der zentrierten Abweichungen (s.o.) einer Stichprobe der
(Zeilen-) Lange m und der (Spalten-) Vektorendimension q. Dann ist der ML Schatzer der
Kovarianzmatrix gegeben durch
Σ = S =1
mDDT (6.5)
Der Rang der Matrix S bestimmt, wieviel Eigenwerte ungleich Null sind. Da rg(S) = rg(D) ≤min(m, q) ist, gibt es also auch hochsten min(m, q) viele Eigenwerte, die von Null verschieden
63
6 Empirische Orthogonalfunktionen
sind. Ist also m > q, so sind wahrscheinlich alle Eigenwerte der geschatzten Kovarianzmatrix
grosser als Null und konnen durch ein Standarddiagonalisierungsprogramm bestimmt wer-
den.
Ist dagegen m < q, so muß man dafur sorgen, dass auch nur m−1 viele Eigenwerte/vektoren
bestimmt werden, die ungleich Null sind. Dazu benutzt man/frau folgendes Verfahren. Wenn
~uj ein Eigenvektor der Matrix S ist, ist der Vektor DT ~uj ein Eigenvektor der Matrix 1m
DT D.
Beweis: Die Eigenwertgleichung lautet
1
mDDT ~uj = λj~uj (6.6)
Multiplikation von links mit DT ergibt die Eigengleichung fur die Behauptung.
1
m(DT D)DT ~uj = λjD
T ~uj (6.7)
Da die Eigenvektoren ~uj linear unabhangig sind, sind auch die neuen Eigenvektoren linear
unabhangig. Beweis: Lineare Unabhangigkeit heißt fur alle a, b ∈ IR
aDT ~uj + bDT ~uk = 0 (6.8)
nur wenn a = b = 0. Multiplikation der rechten Seite der letzten Gleichung von links mit D
ergibt:
aDDT ~uj + bDDT ~uk = aλj~uj + bλk~uk = 0 (6.9)
wobei der letzte Schritt wegen der linearen Unabhangigkeit der ursprunglichen Eigenvektoren
folgt. Der eben bewiesenen Satz zeigt, dass im Fall m < q man nicht die Matrix 1m
DDT
diagonalisieren muss , sondern die Matrix 1m
DTD. Die Eigenvektoren dieser Matrix seien mit
~vj bezeichnet, dann ergeben sich die gesuchten Eigenvektoren ~uj als Losung der Gleichung
~vj = DT ~uj (6.10)
Dies ist ein unterbestimmtes Gleichungssystem fur die q Unbekannten ~uj. Damit gibt es keine
eindeutige Losung, sondern die Losung ist ein Unterraum des IRq. Um die Losung eindeutig
zu machen, muss eine Zusatzbedingung gestellt werden, die ublicherweise wie folgt lautet
1
2~uT
j ~uj!= min (6.11)
Da aber auch Gl. (6.10) gelten muss, mussen wir einen Lagrange’schen Multiplikatorvektor
~γ einfuhren und das Minimum folgender Zielfunktion bestimmen
J =1
2~uT
j ~uj + γT (~vj − DT ~uj)!= min (6.12)
64
6 Empirische Orthogonalfunktionen
Bilden wir den Gradienten nach den beiden unbekannten Vektoren, so erhalten wir
∂J∂~uj
= ~uj − D~γ = 0
∂J∂~γ
= ~vj − DT ~uj = 0 (6.13)
Einsetzen der ersten Gleichung in die zweite ergibt eine Bestimmungsgleichung fur den La-
grange’schen Multiplikatorvektor
DT D~γ = ~vj (6.14)
Da aber ~vj ein Eigenvektor zum Eigenwert λj ist, lautet die Losung
~γ =1
mλj
~vj (6.15)
bzw. die Losung fur den gesuchten Vektor ~uj
~uj =1
mλj
D~vj (6.16)
Der Vorfaktor ist aber unerheblich, da die Vektoren auf den Betrag 1 normiert werden.
Deshalb werden aus den berechneten Eigenvektoren ~vj die Schatzer der prinzipalen Vektoren
(die EOF’s) ~uj berechnet zu
~uj =D~vj
|D~vj|(6.17)
Ist m womoglich sehr viel geringer als q, ist die Diagonalisierung der m × m - Matrix DT D
auch sehr viel weniger aufwendig als die entsprechende Diagonalisierung der q × q Matrix
DDT .
Die Schatzungen der Prinzipalen Komponenten oder die Berechnung der Amplituden der
EOF’s erfolgt wie vermutet. Dazu seien die ersten q Amplituden der EOF’s einer Stichprobe
der Lange m ebenfalls als q × m Matrix A angeordnet. Dann erhalt man/frau diese Matrix
aus der Stichprobenmatrix der zentrierten Abweichungen D
A = UT D (6.18)
wobei U die Matrix der Eigenvektoren der geschatzen Kovarianzmatrix ist. Man sollte sich
aber davor huten anzunehmen, dass
• die geschatzten Prinzipalen Vektoren die Eigenvektoren der Kovarianzmatrix der GG
Σ sind
65
6 Empirische Orthogonalfunktionen
• die Schatzer der prinzipalen Komponenten – die EOF Amplituden – unkorreliert sind.
• die geschatzten Eigenwerte der Kovarianzmatrix unverzerrt sind
D.h. im allgemeinen gilt:
E(aiaj) 6= 0
ΣU 6= UΛ
E(λi) 6= λi
V ar(ai) 6= λi (6.19)
Bei den letzten Beziehungen ist es sogar moglich, Ungleichungen anzugeben fur verschieden
λi:
V ar(ai) < λi < E(λi) λi gross
V ar(ai) > λi > E(λi) λi klein (6.20)
D.h. die grossen Eigenwerte sind positiv verzerrt (uberschatzt) wahrend die kleinen Eigen-
werte unterschatzt werden (wie alles im richtigen Leben).
6.3 EOF und die Singularwert-Zerlegung SVD
Ein Theorem der linearen Algebra besagt, dass eine beliebige q×m Matrix D immer zerlegt
werden kann in ein Produkt aus drei Matrizen U, Λ und V . Sei ohne Einschrankungen q <
m ≤ rang(D). Dann ist U eine Matrix der Dimension q×m, Γ eine Diagonalmatrix mit den
Eintragen γi 6= 0, i = 1 . . . rang(D) und γi = 0, i = rang(D) + 1 . . .m, sowie V eine Matrix
der Dimension m × m. Dann gelten folgende Beziehungen:
D = UΓV T
V V T = V T V = Im
UT U = Im (6.21)
Dabei ist Im die Einheitsmatrix der Dimension m × m. Die Diagonaleintrage der Matrix Γ
nennt man auch die singularen Werte der Matrix D und die obige Aufspaltung der Matrix
D die singular value decomposition (SVD) von D. U nennt man dann die Matrix der linken
Eigenvektoren und V die der rechten. Die Beziehungen zwischen den Spaltenvektoren der
66
6 Empirische Orthogonalfunktionen
Matrizen U und V geben an, dass die Spaltenvektoren sowohl der linken als auch der rechten
Eigenvektoren orthonormal sind. Zusatzlich sind auch die Zeilenvektoren der Matrix V or-
thonormale Vektoren: V ist eine orthogonale Matrix. Die SVD ist die Verallgemeinerung der
Diagonalisierung von quadratischen Matrizen auf beliebige, rechteckige Matrizen. Sei nun
D die Matrix der zentrierten Abweichungen der Stichprobe. Dann ist der ML Schatzer der
Kovarianzmatrix gegeben als
S =1
mDDT (6.22)
Die Schatzer der prinzipalen Vektoren U und und ihrer Eigenwerte sind dann Losung der
Eigengleichung1
mDDT U = U Λ (6.23)
Setzen wir nun die SVD der Datenmatrix D ein, erhalten wir
1
mUΓ V T V︸ ︷︷ ︸
=Ir
ΓUT U = U Λ
1
mUΓ2UT U = U Λ (6.24)
Setzen wir nun als Losung U = U an, erhalten wir sofort
1
mUΓ2 UT U︸ ︷︷ ︸
Ir
= UΛ (6.25)
bzw. fur die Eigenwerte der EOF’s
Λ =1
mΓ2 (6.26)
D.h. die EOF’s als Schatzer der prinzipalen Vektoren der Kovarianzmatrix konnen wir auch
durch die SVD der Datenmatrix bestimmen:
• die EOF’s sind die linken Eigenvektoren der Datenmatrix U = U
• die Eigenwerte erhalten wir aus den singularen Werten durch die Beziehung Λ = 1m
Γ2.
6.4 Schatzung der Konfidenzintervalle
Nachdem wir nun die ”Punkt” Schatzung von prinzipalen Vektoren diskutiert haben, wollen
wir noch zusatzlich bestimmen, wie genau die Schatzungen sind. Im univariaten Fall wurde
dies durch die Konfidenzintervalle der Schatzer angegeben. Etwas ahnliches kann auch fur
die EOF’s hergeleitet werden. Wir erhalten dann (Schatzungen der ) Konfidenzintervalle fur
67
6 Empirische Orthogonalfunktionen
die Eigenwerte und die Eigenvektoren. Wenn S die geschatzte Kovarianzmatrix ist, konnen
wir auch schreiben
S = Σ + εS ′
ε = O(1
m) (6.27)
so dass ε eine kleine Zahl ist. Gesucht werden typische (zufallige) Abweichungen der geschatz-
ten Eigenvektoren ~ ju und der geschatzten Eigenwerte λj von ihren GG Partnern ~uj und λj
aufgrund des Schatzfehlers S ′. Ordnen wir die Vektoren als Spaltenvektoren in eine Matrix
ein und schreiben die Eigenwerte als die Diagonaleintrage einer Diagonalmatrix, so gilt fur
die GG
ΣU = Uλ (6.28)
wahrend die Eigenwertgleichung fur die Schatzer lautet
SU = (Σ + εS ′)U = U λ (6.29)
Da die prinzipalen Vektoren (PV) der Kovarianzmatrix Σ ein orthogonales und vollstandi-
ges Basissystem bilden, konnen wir die Schatzer der PV durch diese Basis darstellen (entwi-
ckeln und die Entwicklungskoeffizienten in eine Matrix A einordnen)
U = UA (6.30)
Da UUT = I ist, gilt auch
A = UT U (6.31)
Setzen wir diesen Ansatz in die Eigenwertgleichung fur die Schatzer (6.29) ein, ergibt sich
mit Hilfe der Eigengleichung fur U und der Orthogonalitat UT U = I
λA + εUT S ′UA = Aλ (6.32)
Wir entwickeln nun sowohl die geschatzten Eigenwerte λ als auch die Entwicklungskoeffi-
zienten der Eigenvektoren in eine Reihe in Ordnung ε
A = I + εA′ + O(ε2)
λ = λ + ελ′ (6.33)
68
6 Empirische Orthogonalfunktionen
Einsetzen dieser Storungsentwicklung ergibt in Ordnung O(ε)
λ′I + (A′λ − λA′)︸ ︷︷ ︸=B
= UT S ′U (6.34)
Dabei ist die Matrix B eine Matrix mit den Diagonalelementen Null, so dass die Abwei-
chungen der Eigenwerte in Ordnung ε gegeben sind durch die Diagonalelemente der Matrix
UT S ′U
λ′j = ~uT
j S ′~uj (6.35)
Da S ′ ein Maß fur die Varianz einer Wishard verteilten Zufallsmatrix S ist, kann man
zeigen (siehe Anderson (1985), [10], p.540), dass (asymptotisch m → ∞) die Eigenwerte
normalverteilt mit Erwartungswert Null und Varianz 2m
λ2j sind
λ′j ∈ N(0,
2
mλ2
j) (6.36)
Damit konnen wir ein approximiertes (95%) Mutungsintervall fur den Eigenwert λj angeben
λj − 2
√2
mλ2
j ≤ λ ≤ λj + 2
√2
mλ2
j (6.37)
Analog kann man/frau zeigen, dass die Kovarianzmatrix der Abweichungen ~u′j = ~uj − ~ ju
gegeben ist durch
E(~u′j(~u
′j)
T ) =
q∑
k=1,k 6=j
λjλk
(λj − λk)2~uj~u
Tk (6.38)
bzw. als Kovarianzmatrix zwischen zwei Abweichungen ~u′j und ~u′
k , (k 6= j)
E(~u′j(~u
′k)
T ) =λjλk
(λj − λk)2~uj~u
Tk (6.39)
Beide Gleichungen zeigen, dass ein geschatzer PV eine hohe Varianz und auch eine grosse
Kovarianz mit anderen geschatzen PV hat, wenn zwei Eigenwerte nahe beieinander liegen
λj ∼ λk. Man / frau spricht in diesem Fall von entarteten (λj = λi) bzw. nahezu entarte-
ten Eigenwerten und -vektoren. Eine Uberprufung auf Entartung erfolgt naturlich mit Hilfe
der Schatzungen λi. Uberlagern sich die Mutungintervalle (6.37) zweier benachbarter Eigen-
wertschatzungen, so kann man/frau eine Entartung stark vermuten. Diese Regel ist in der
meteorologischen Literatur als North’sche Daumenregel (rule-of-thumb) bekannt geworden,
da sie in der Veroffentlichung von North et al. (1982) im Monthly Weather Review einem
weiteren Publikum zuganglich wurde. Allerdings basiert dieser spezielle Artikel auf einer rus-
sischen Publikation. Ausserdem findet man die entsprechende Herleitung auch in Anderson
(1984) ([10]), das in der ersten Auflage bereits 1957 erschien.
69
6 Empirische Orthogonalfunktionen
6.5 Beweis des Satzes uber prinzipale Vektoren
Sei ~v1 . . . ~vk ein beliebiges Orthonormalsystem mit k ≤ q. Dann definieren wir eine Residu-
umsfunktion R(k, ~x) durch
R(k, ~x) = |~x −k∑
j=1
aj(~x)~vj|2 (6.40)
wobei wie oben
aj(~x) = ~xT~vj (6.41)
die Projektion des Datenvektors ~x auf den Basisvektor ~vj darstellt. Ausrechnen der Residu-
umsfunktion ergibt:
R(k, ~x) = (~xT~x) − 2
k∑
j=1
aj(~x)(~xT~vj) +
k∑
i,j=1
ai(~x)aj(~x)(~vTj ~vi) (6.42)
Die Matrix (~vTj ~vi) ist wegen der Orthonormalitat der Basisvektoren aber die Einheitsmatrix
δi,j. Deshalb reduziert sich die letzte Gleichung zu
R(k, ~x) = (~xT~x) −k∑
j=1
aj(~x)(~xT~vj) (6.43)
Der Summand lasst sich auch schreiben als
aj(~x)(~xT~vj) =
(~xT~vj)(~xT~vj) =
~vTj (~x~xT )~vj (6.44)
Nimmt man/frau nun von der Residuumsfunktion den Erwartungswert, folgt
E(R(k, ~x)) = E(|~x|2) −k∑
j=1
~vTj E(~x~xT )~vj =
E(|~x|2) −k∑
j=1
~vTj Σ~vj (6.45)
wobei – wie gehabt – Σ die Kovarianzmatrix der Vektor ZVA ~x ist. Die Behauptung des
zu beweisenden Satzes ist nun, dass der Erwartungswert der Residuumsfunktion den kleins-
ten Wert annimmt, wenn die Orthonormalbasis die Eigenvektoren der Kovarianzmatrix sind.
Wegen der Abhangigkeit von der naturlichen Zahl k bietet sich ein Beweis durch vollstandige
70
6 Empirische Orthogonalfunktionen
Induktion an. Sei also k = 1 (Induktionsanfang). Man/frau sucht ein ~v1, so dass der Erwar-
tungswert der Residuumsfunktion minimal ist und die Nebenbedingung |~v1|2 = 1 erfullt ist.
Dazu fuhrt man/frau einen Lagrange’schen Multiplikator g ein und sucht das Minimum der
Funktion
F (~v1, g) = E(R(1, ~x)) + g(|~v1|2 − 1) (6.46)
bezuglich der Variablen ~v1 und g. Dazu bilden wir den Gradienten von F nach ~v1 und setzten
diesen gleich Null:
∇~v1(E(|~x|2 − ~vT1 Σ~v1) + g~vT
1 ~v1 − g) = −2Σ~v1 + 2g~v1 = 0 (6.47)
Dies fuhrt also auf eine Eigenwertgleichung fur das Tupel ~v1, g
Σ~v1 = g~v1 (6.48)
von der es genau q viele linear unabhangige Losungen gibt. Die Entscheidung, welche Losung
die richtige ist, ergibt der Wert der Minimalfunktion am Minimum. Dieser ist gegeben durch
F = E(|~x2|) − g (6.49)
Damit ist die Losung des Minimalproblems eindeutig bestimmbar: F nimmt den geringsten
Wert an, wenn g der grosste Eigenwert der Kovarianzmatrix Σ ist. Dann ist auch ~v1 = ~u1 der
Eigenvektor der Kovarianzmatrix zum grossten Eigenwert. Damit ist der Induktionsanfang
bewiesen. Gelte nun die Behauptung fur k, zu zeigen ist nun, dass sie dann auch fur k + 1
richtig ist. Wir betrachten die gleiche Minimalaufgabe, ausser dass nun neben der Neben-
bedingung |~vk+1|2 = 1 auch noch die zusatzlichen k NB ~vTj ~vk+1 = 0, j = 1 . . . k zu erfullen
sind. Dazu fuhrt man/frau qk weiter Lagrange’sche Multiplikatoren qj ein und minimiert die
Funktion
F (~vk+1, g, qk) = E(R(1, ~x)) + g(|~vk+1|2 − 1) +k∑
j=1
qj~vTk+1~vj (6.50)
Der Gradient bezuglich ~vk+1 ergibt die Gleichung
−2Σ~vk+1 + 2g~vk+1 +k∑
j=1
qj~vj = 0 (6.51)
Multiplikation dieser Gleichung von links mit ~vTk+1 ergibt unter Berucksichtigung, dass die
ersten k Vektoren ~vj die orthonormalen Eigenvektoren der symmetrischen Kovarianzmatrix
71
6 Empirische Orthogonalfunktionen
sind
2~vTj (Σ − g)~vk+1 + qj =
2((Σ − g)~vj)T~vk+1 + qj =
2((λj − g))~vj)T~vk+1 + qj = qj
⇒ qj = 0 (6.52)
D.h. die Orthogonalitatsforderungen zwischen ~vk+1 und ~vj, j = 1, . . . , k sind automatisch
erfullt. Damit reduziert sich Gl. (6.51) auf die Eigenwertgleichung (6.48) und als Losungen
stehen wiederum nur die Eigenvektoren der Kovarianzmatrix zur Verfugung. Die ersten k
Eigenvektoren erfullen die Orthogonalitatsrelation und um die Funktion F minimal werden
zu lassen, muss g der Eigenwert λk+1 sein und ~vk+1 entsprechend der Eigenvektor ~uk+1.
Damit ist der Induktionschritt bewiesen und der gesamte Beweis abgeschlossen.
6.6 Rotation
In vielen Arbeiten werden die EOF’s noch einer sogenannten Rotation unterworfen. Das Ziel
ist es, die zwei- oder dreidimensionalen Strukturen, die sich in typischen EOF’s wiederfinden
in einfach zu interpretierbare Muster uberzufuhren. Sei zunachst wieder die Eigengleichung
fur die geschatzte Kovarianzfunktion
S =1
mDDt
SU = UΛ
U tU = I
A = U tD (6.53)
Mathematisch bedeutet dies die Einfuhrung einer zunachst unbekannten Rotationsmatrix
R, die auf die Eigenvektoren U wirkt und ein neues System von Basisvektoren P erzeugt
P = RU (6.54)
Eine orthogonale Rotation RtR = I erhalt somit die Orthogonalitat der ursprunglichen
Eigenvektoren aber es sind auch Rotationen vorstellbar, die dies nicht erzwingen (”oblique”).
Die Rotationsmatrix R wird dann durch eine zu minimierende Zielfunktion an die neuen
72
6 Empirische Orthogonalfunktionen
Vektoen P bestimmt
V(P )!= min (6.55)
Dabei ist V eine Vorschrift an die einfache Struktur der neuen Basisvektoren. Eine beliebte
Funktion ist die Varimax Vorschrift
V(~p1, . . . , ~pK) =K∑
i=1
fV(~pi)
fV =1
q
q∑
j=1
(pi(j)
sj
)4 − 1
q2
(q∑
j=1
(pi(j)
sj
)2
)2
(6.56)
Hierbei ist pi(j) die j-te Komponente des Vektors ~pi und sj ein vom Benutzer vorzugebender
Wichtungswert z.B.
sj = 1
sj =K∑
i=1
(pi(j))2 (6.57)
Die Meinungen zur Rotation sind allerdings sehr weit auseinander. Ein Teil der Commu-
nity schwort heftigst auf die Rotation (” You have not rotated the EOF’s, therefore we can
not accept your paper for publication”), ein anderer Teil (zu dem auch ich gehore) sind die
zum Teil wenig fundierten und nur heuristisch begrundeten Grundlagen der Rotation zu
vage, um das mathematisch wohl fundierte Optimierungskalkul ausser Kraft zu setzen. Ro-
tation scheint angebracht, wenn mehrer Eigenwerte nahe beieinander liegen (”entartet”). In
diesem Fall definieren die Eigenvektoren nur einen Unterraum, in dem jede beliebige Linear-
kombination eine Losung des Minimierungsproblems ist. In diesem Fall kann die zusatzliche
Forderung nach einfachen Strukturen unter einer orthogonalen Rotation Ordnung in den
Unterraum bringen und bestimmte Richtungen auszeichnen.
6.7 Singular Spectral Analysis
Ein Abart der EOF Analyse geht auf Fraedrich (1987) zuruck und ist durch Vautard et al.
(1992) verfeinert worden. Diese nennt sich Singular Spectral Analysis (SSA). Ursprunglich
nur fur univariate Zeitreihen gedacht, ist diese Methode dann auch von Vautard (1995) auf
Felder angewandt worden und firmiert hier unter dem titel Multichannel Singular Spectral
Analysis. Wir wollen uns aber hier nur der einfachen SSA zuwenden. Dazu betrachten wir
73
6 Empirische Orthogonalfunktionen
eine univariate Zeitreihe xt, t = 1, T . Ein q dimensionaler (Spalten)-Vektor ~Xt entsteht durch
die Anordnung von q sequentiellen Werten aus der Zeitreihe
~Xt = (xt, xt+1, . . . , xt+q−1)t (6.58)
Eine derartige Konstruktion wird nach Takens (1980) ein delay-Koordinaten Vektor ge-
nannt.
Als nachstes wird dieser Vektor wie ein normaler multivariater Vektor behandelt. Die
Kovarianzmatrix Σ = E( ~Xt~XT
t ) ist eine sogenannte Toplitzmatrix, da entlang der Bander
parallel zur Hauptdiagonalen alle Werte identisch sind zum Wert der Autokovarianzfunktion
der univariaten Zeitreihe zur Zeitverschiebung τ = 0, 1, . . . , q − 1. Die Eigenwerte dieser
Kovarianzmatrix nennt man die singularen Spektralkoeffizienten und die Eigenvektoren kann
man als typische Zeitstrukturen in der Zeitreihe interpretieren. Fur verschiedene typische
Zeitprozesse kann man aus der Struktur der Toplitzmatrix die Struktur der EOF’s direkt
erschliessen. So ist ein Weisses Rauschen der Varianz σ2ε )
xt = εt (6.59)
mit
E(εt) = 0
E(εtεs) = σ2ε δts (6.60)
durch die Kovarianzmatrix
Σ = σ2ε I (6.61)
mit q identischen Eigenwerten und den Eigenfunktionen
~ui = (0, . . . , 1︸︷︷︸i
, 0, . . .) (6.62)
ausgezeichnet.
Ein random walk Prozess (Wiener Prozess)
xt+1 = xt + εt (6.63)
(εt wie oben)
hat die Eigenfunktionen
~ui =2√
2q
(2i + 1)πsin(
(2i + 1)πt
2q) t ∈ [0, q − 1] (6.64)
(siehe Kloeden und Platen, 1999, [17] ).
74
6 Empirische Orthogonalfunktionen
6.8 Beispiel fur eine EOF Analyse
Die wichtigste Entscheidung bei der Berechnung von EOF’s ist die Wahl der Schatzmethode
nach Gl.(6.6) fur den Fall m > q oder nach Gl.(6.7) und (6.10) fur den Fall m ≤ q. Da q in
den allermeisten Fallen die Anzahl der Gitterpunkte ist, ist die Berechnung meist nur nach
dem letzteren Verfahren moglich. Wir werden an dem folgenden Beispiel sehen, welche Feh-
ler man sich einhandelt, wenn die falsche Methode gewahlt wird. Die zweite Frage, die man
sich vor Beginn einer EOF - Analyse stellen muss, ist die nach der Geometrie des Problems.
Man kann anhand der raumlich kontinuierlichen Formulierung der Karhunen-Loeve Funktio-
nen namlich zeigen, dass die Eigen-Gleichungen diskretisierte Integralgleichungen sind. Die
diskreten Summen sind approximierte Integrale ubder den jeweilig betrachteten Raum. Bei
zweidimensionalen, globalen horizontalen Feldern ist dies die Kugel, so dass die Integrale in
Kugelgeometrie bestimmt werden mussen und dies sich auch in der diskreten Darstellung
wiederfinden muss. Allgemein kann man dies durch eine Metrik-Matrix W berucksichtigen.
Dann lautet die Eigengleichung (6.6):
1
mDDTW~uj = λj~uj (6.65)
oder kompakt in Matrixform1
mDDT WU = U Λ (6.66)
Multiplikation von links mit DT W ergibt
1
mDTWD (DT WU)︸ ︷︷ ︸
=U
= DT WU Λ (6.67)
so dass wir die Eigengleichung
1
mDtWDU = UΛ (6.68)
losen. Zur Ruckrechnung mussen wir dann noch die Minimierungsaufgabe
J =1
2~uT
j W~uj + ~γT (~vj − DT W~uj)!= min (6.69)
75
7 EOF’s und ihre physikalische Interpretation
losen, was aber wie oben zur gleichen Form der Losung fuhrt:
D~Γ = ~vj
~vj = DT W~uj
~Γ =1
λj
~vj
~uj =1
λD~vj (6.70)
so dass wir wie bereits fur den Fall ohne Metrik die Eigenvektoren der Matrix DTWD von
links mit der Datenmatrix multiplizieren mussen.
Die Anwendung auf die bodennahen Lufttemperaturen zwischen 1881 und 1998 (m = 118
Jahre) mit einer Auflosung von 5◦ × 5◦ zwischen 32.5◦S und 77.5◦N (q = 72× 23 = 1656) ist
in den folgenden Abbildungen zusammengefasst. In Abb.(10) sehen wir die normalisierten
Spektren der Eigenwerte (Eigenwerte normalisiert mit der Gesamtvarianz, Eigenmodenum-
mer normalisiert mit der maximalen Zahl der von Null verschiedenen Eigenwerte m) fur
gleitende 3-Monatsmittel beginnend in DJF bis NDJ. Man erkennt, dass die Spektren fur die
Sommermonate flacher abfallen (”weisser”) als die fur die Winterjahreszeiten. In Abb.(11)
finden wir fur den 3-Monatszeitraum August-September-Oktober ASO die Verteilung der
Standardabweichung der bodennahen Lufttemperaturen sowie die ersten drei EOF’s.
7 EOF’s und ihre physikalische Interpretation
Eine Warnung vor EOF’s sollte man auf jeden Fall aussprechen. Sie sind ob ihrer Berech-
nung zunachst nur reine mathematische Konstruktionen, die einem bestimmten optimalen
Charakter aufweisen. Man darf aber auf keinen Fall ohne zusatzliche Informationen diesen
Moden einen physikalischen Inhalt geben. Dies gilt insbesondere fur die zweite und jede
weitere Mode, die durch die Zwangsbedingung der Orthogonalitat massiv eingeschrankt in
ihrer Interpretation sind. Man kann anhand eines einfachen Beispiels vorfuhren, dass generell
der Zwang zur Orthogonalitat auch in der ersten Mode zu Problemen furhen kann. Dazu
stellen wir uns einen Prozess ~x vor aus drei nicht-orthogonalen Mustern ~mi, i = 1, 3, deren
Amplituden ai, i = 1, 3 mit einer bestimmten Varianz schwanken und untereinander nicht
korreliert sind. Zur Vereinfachung wollen wir uns auf einen dreidimensionalen Vektorraum
76
7 EOF’s und ihre physikalische Interpretation
Abbildung 10 Normalisierte Eigenwertspektren der EOF Analyse der bodennahen Lufttempera-
turen zwischen 1881 und 1998
konzentrieren:
~x =
3∑
i=1
ai ~mi
Σ = E(~x~xT ) =3∑
i=1
σ2(ai)~mi ~mTi (7.1)
Geben wir die Moden ~mi und die Varianzen σ2(ai)vor, konnen wir aus der Summe der
Elementarmatrizen ~mi ~mTi die Kovarianzmatrix Σ bestimmen und diese diagonalisieren. Dere
Vergleich der Eigenvektoren mit den vorgegebenen Moden zeigt dann, wie die Moden durch
die Orthogonalisierung verandert werden.
Wahlen wir
~m1 = (1, 0, 0)T
~m2 = (0, 0, 1)T
~m3 =1√3(1, 1, 1)T (7.2)
77
7 EOF’s und ihre physikalische Interpretation
und
σ2(a1) = 0.436
σ2(a2) = 0.353
σ2(a3) = 0.209 (7.3)
so ergeben sich folgende Eigenvektoren
~e1 = (0.84, 0.19, 0.51)T
~e2 = (−0.53, 0.07, 0.85)T
~e3 = (0.13,−0.97, 0.16)T (7.4)
und den Eigenwerten
σ2(a1) = 0.564
σ2(a2) = 0.385
σ2(a3) = 0.05 (7.5)
Man erkennt deutlich, das die blosse Anwesenheit eines verhaltnismassig schwachen aber
nichtorthogonalen Modes ~m3 zu Eigenvektoren fuhrt, die mit dem ursprunglichen Mustern
nicht im geringsten zu vergleichen sind.
78
7 EOF’s und ihre physikalische Interpretation
Abbildung 11 Standardabweichung der bodennahen Temperaturen im 3 Monatszeitraum ASO
zwischen 1881 und 1998 (links oben) sowie die raumliche Struktur der ersten drei EOF’s
79
8 Kanonische Korrelationsanalyse
8 Kanonische Korrelationsanalyse
Im zweiten Kapitel dieser Vorlesung haben wir als eine spezielle Eigenschaft bivariater
Grundgesamtheiten die Korrelation zwischen zwei ZVA X und Y eingefuhrt als ein Mass
des linearen Zusammenhangs (Modells) zwischen den beiden Grossen. Die Frage, die sich
nun nach der Betrachtung allgemeiner multivariater GG aufdrangt, ist: Kann man die Kor-
relation zwischen zwei ZVA verallgemeinern auf vektorwertige ZVA ~X und ~Y , so dass diese
Verallgemeinerung wiederum ein Mass fur den linearen Zusammenhang (lineare Abbildung)
zwischen diesen beiden Vektoren darstellt? Dies ist die Aufgabe der sogenannte kanonischen
Korrelationsanalyse, die von Hotelling in den 30’er Jahren entwickelt wurden und die sich
ebenfalls steigender Beliebheit in der Klimatologie erfreut.
8.1 Definition
Wir betrachten eine normalverteilte ZVA ~X mit dem Erwartungswert ~µx und Kovarianz-
matrix Σxx, sowie eine zweite Normalverteilte ZVA ~Y mit E(~Y ) = ~µy und Kovarianzmatrix
Σyy. Die Kovarianzmatrix zwischen den beiden ZVA (als Verallgemeinerung der Korrelation
ρxy) kann man dann berechnen aus
Σxy = E(( ~X − ~µx)(~Y − ~µy)T ) (8.1)
Ohne Einschrankung der Allgemeinheit sei nun angenommen, dass die Erwartungswerte
verschwinden ~µx = ~µy = 0. Gesucht sind nun (guess-) Muster ~α fur ~X und ~β fur ~Y , so dass
die Projektionen (Entwicklungskoeffizienten)
U = ~αT ~X
V = ~βT ~Y (8.2)
maximal korreliert sind. Die Korrelation ρuv ist definiert als
ρuv =E(UV )√
E(U2)E(V 2)(8.3)
Dieser Ausdruck wird besonders einfach, wenn zusatzlich zur maximalen Korrelation gefor-
dert wird
E(U2) = E(V 2) = 1 (8.4)
80
8 Kanonische Korrelationsanalyse
Die Definitionen fur die Projektionen (8.2) sind als Euklidisches Skalarprodukt definiert.
Im Allgemeinen wird ein Skalarprodukt durch eine Metrik W festgelegt, so dass wir Gl. (8.2)
noch schreiben konnen
U = ~αT Wx~X
V = ~βTWy~Y (8.5)
mit den zunachst noch frei wahlbaren (symmetrischen) Metrikmatrizen Wx, Wy. Die Nor-
mierungsbedingungen lauten dann
E(U2) = E((~αTWx
~X)(~αT Wx~X)T
)
= E(~αT Wx
~X ~XT Wx~α)
= ~αTWxE( ~X ~XT )Wx~α
= ~αT WxΣxxWx~α
E(V 2) = ~βT WyΣyyWy~β (8.6)
Die Extremalbedingung lautet dann analog
E(UV ) = ~αT WxΣxyWy~β
!= max (8.7)
Wir haben hier eine Extremwertaufgabe mit zwei Nebenbedingungen zu losen. Dies erfolgt
mit der Methode der Lagrange’schen Multiplikatoren, in dem das Maximum der Funktion
J gesucht wird mit
J = ~αT WxΣxyWy~β − γ1
2(~αT WxΣxxWx~α − 1) − γ2
2(~βTWyΣyyWy
~β − 1) (8.8)
Unbekannt sind die Vektoren ~α und ~β sowie die beiden skalaren Lagrange-Multiplikatoren
γ1,2. Bilden wir die Ableitungen der Funktion J nach den Unbekannten (bzgl. der beiden
Vektoren nach den Regeln vom Beginn der Vorlesung), so ergibt sich
~∇αJ = WxΣxyWy~β − γ1WxΣxxWx~α = 0
~∇βJ = WyΣTxyWx~α − γ2WyΣyyWy
~β = 0
∂J∂γ1
= ~αTWxΣxxWx~α − 1 = 0
∂J∂γ2
= ~βTWyΣyyWy~β − 1 = 0 (8.9)
81
8 Kanonische Korrelationsanalyse
Zur Berechnung der γ1,2 multiplizieren wir die erste Gleichung von links mit ~αT und die
zweite mit ~βT und benutzen die beiden Nebenbedingungen
~αTWxΣxyWyβ − γ1 ~αT WxΣxxWxα︸ ︷︷ ︸=1
= 0
γ1 = ~αT WxΣxyWyβ
βTWyΣTxyWxα − γ2 βTWyΣyyWyβ︸ ︷︷ ︸
=1
= 0
γ2 = βT WyΣTxyWxα = γ1 = γ (8.10)
Die beiden Lagrange’schen Multiplikatoren sind also gleich. Damit lautet das Gleichungs-
system zur Bestimmung der unbekannten Muster ~α und ~β nun:
WxΣxyWy~β − γWxΣxxWx~α = 0
WyΣTxyWx~α − γWyΣyyWy
~β = 0 (8.11)
Da sowohl Σxx als auch Σyy den vollen Rang haben und die Matrizen Wx bzw. Wy Metriken
sind und damit auch den vollen Rang haben, konnen wir beide Gleichungen jeweils mit dem
Inversen
(WxΣxxWx)−1
(WyΣyyWy)−1 (8.12)
von links multiplizieren.
(WxΣxxWx)−1WxΣxyWy
~β − γWxΣxxWx~α = 0
W−1x Σ−1
xx ΣxyWy~β = γ~α
(WyΣyyWy)−1WyΣ
TxyWx~α − γWyΣyyWy
~β = 0
W−1y Σ−1
yy ΣTxyWx~α = γ~β (8.13)
Multiplizieren wir beide Gleichungen nun noch mit γ, so konnen wir dann wechselweise
einsetzen und erhalten als Losung
W−1x Σ−1
xx ΣxyΣ−1yy ΣT
xyWx~α = γ2~α
W−1y Σ−1
yy ΣTxyΣ
−1xx ΣxyWy
~β = γ2~β (8.14)
82
8 Kanonische Korrelationsanalyse
Die gesuchten Muster sind also Eigenvektoren der entsprechenden Matrixkombinationen zum
gemeinsamen Eigenwert γ2. Die Eigenwerte sind also alle positiv und sind gegeben durch
Gl. (8.10)
γ2 = (~βTWyΣTxyWx~α)2 (8.15)
Da aber U = ~αT Wx~X und V = ~βWy
~Y ist, ist γ2 identisch (E(UV ))2, dem Quadrat der
Funktion, die maximiert werden sollte. Entsprechend wird die gesuchte Losung ~α, ~β durch
die Eigenvektoren mit dem grossten Eigenwert festgelegt, da wir ja die Korrelation zwischen
U und V maximieren wollten.
Seien nun diese Losungen bestimmt. Dann ergibt sich die Frage, ob man diese Maximie-
rungsprozedur nochmals durchfuhren kann, um auf ein weiteres Paar von Muster zu kommen,
die die Daten ~X und ~Y abzuglich der im ersten Schritt bestimmten Varianzanteile optimal
korreliert. Dazu nehmen wir die im ersten Schritt bestimmten Eigenvektoren ~α1 und ~β1 und
die neu zu bestimmenden Muster ~α2 bzw. ~β2 mit
U1 = ~αT1 Wx
~X
V1 = ~βT1 Wy
~Y
U2 = ~αT2 Wx
~X
V2 = ~βT2 Wy
~Y (8.16)
Wie im ersten Schritt sollen auch die neuen Muster so bestimmt werden, dass
E(U22 ) == ~αT
2 WxΣxxWx~α2 = 1
E(V 22 ) == ~βT
2 WyΣyyWy~β2 = 1 (8.17)
gilt. Zusatzlich mussen wir aber noch spezifizieren, wie die Beziehung zwischen U1 und U2
bzw. V1 und V2 aussehen soll. Dafur bietet sich an, die Projektionen so zu wahlen, dass die
Entwicklungskoeffizienten unkorrelliert (statistisch orthogonal) sind.
E(U1U2) = ~αT1 WxΣxxWx~α2 = 0
E(V1V2) = ~βT1 WyΣyyWy
~β2 = 0 (8.18)
zu setzen. Mit Hilfe der Gleichungen (8.11) lasst sich dann sehr einfach zeigen, dass dann
auch gilt
E(U2V1) = ~αT2 WxΣxyWy
~β1 = γ~αT2 WxΣxxWx~α1 = 0
E(V2U1) = ~βT2 WyΣ
TxyWx~α1 = γ~βT
2 WyΣyyWy~β1 = 0 (8.19)
83
8 Kanonische Korrelationsanalyse
Damit lautet die Extremalbedingung zur Berechnung der Muster ~α2 und ~β2
J2 = ~αT2 WxΣxyWy
~β2 −γ1
2(~αT
2 WxΣxxWx~α2 − 1) − γ2
2(~βT
2 WyΣyyWy~β2 − 1)
+γ3~αT1 WxΣxxWx~α2 + γ4
~βT1 WyΣyyWy
~β2!= max (8.20)
Durch die Berechnung der Gradienten der Funktion J2 nach den Unbekannten erhalten
wir den Satz von Gleichungen
WxΣxyWy~β2 − γ1WxΣxxWx~α2 + γ3WxΣxxWxα1 = 0
WyΣTxyWx~α2 − γ2WyΣyyWy
~β2 + γ4WyΣyyWyβ1 = 0
~αT2 WxΣxxWx~α2 − 1 = 0
~βT2 WyΣyyWy
~β2 − 1 = 0
~αT1 WxΣxxWx~α2 = 0
~βT1 WyΣyyWy
~β2 = 0
(8.21)
Multiplizieren wir die erste Gleichung skalar von links mit ~αT1 und die zweite mit ~βT
1 ,
benutzen die Gl. (8.19) sowie die vorletzte bzw. letzte Gleichung, so erhalten wir
γ3 = γ4 = 0 (8.22)
Multiplizieren wir dagegen die beiden Gleichungen von links skalar mit ~αT2 bzw. ~βT
2 , be-
nutzen die gerade bewiesene Identitat und die Normierungsnebenbedingungen, folgt sofort
γ1 = γ2 = µ2 (8.23)
Damit erhalten wir auch fur die zweite Stufe der kanonischen Korrelationsrechnung die
Eigenwertgleichungen (8.11) fur die Unbekannten ~α2, ~β2 und µ2. Iterieren wir nun den For-
malismus (z.B. mit Hilfe der vollstandigen Induktion, siehe [10]), so werden alle gesuchten
Muster ~αi und ~βi als Eigenvektoren des Gleichungssystems (8.11) zu den Eigenwerten µi
bestimmt. Dabei sind die Eigenwerte der Grosse nach sortiert. Wollen wir wissen, wieviel
Eigenwerte und Eigenvektoren bestimmbar sind, so mussen wir den Rang der Matrizen be-
stimmen. Entscheidend ist dabei der Rang der Kovarianzmatrix Σxy, wie man/frau leicht
aus den aufgelosten Gleichungen (8.14) sieht.
rg(Σxy) = min(q, r) (8.24)
Damit gibt es also hochstens min(q, r) viele Eigenwerte, die von Null verschieden sind.
84
8 Kanonische Korrelationsanalyse
8.1.1 Bedeutung der Metrikmatrizen
Bisher haben wir noch nicht die Metrikmatrizen Wx, Wy spezifiziert. Die naheliegende Wahl
Wx = Iq und Wy = Ir bestimmen die kanonische Korrelationsanalyse bezuglich der Eukli-
dischen Norm. Die Losungen ~αi, ~βi des Systems
Σ−1xx ΣxyΣ
−1yy ΣT
xy~αi = µ2i ~αi
Σ−1yy ΣT
xyΣ−1xx Σxy
~βi = µ2i~βi (8.25)
heissen kanonische Muster und die Amplituden Ui = ~αT~x bzw. V = ~βi~y die kanonischen
Variablen.
Die Abstandberechnung kann aber auch bezuglich des Mahalanobisabstands gemacht wer-
den. Dann wahlen wir
Wx = Σ−1xx
Wy = Σ−1yy (8.26)
Einsetzen in das Gleichungssystem (8.11) ergibt
ΣxyΣ−1yy ΣT
xyΣ−1xx
~iα = µ2
i~
iα
ΣTxyΣ
−1xx ΣxyΣ
−1yy
~iβ = µ2
i~
iβ (8.27)
Die Muster ~iα,
~iβ heissen die kanonischen Korrelationsmuster, da zwischen den kanoni-
schen Mustern und den kanonischen Korrelationsmustern die folgende Beziehung besteht:
~iα = Σxx~αi = E(U~x)
~iβ = Σyy
~βi = E(V ~y) (8.28)
d.h. die kanonischen Korrelationsmuster sind die Muster, die die lokale (komponentenweise)
Kovarianz des betrachteten Datenvektors mit der jeweiligen kanonischen Variablen zusam-
menfassen. Dies sind die ”typischen” Muster der Daten, die optimal im Sinne der kanonischen
Korrelationsanalyse zusammenhangen. Die kanonischen Muster sind dagegen die optimalen
Gewichtungsmuster fur die Datenvektoren, die benotigt werden, die kanonischen Korrelati-
onsmuster zu erzeugen.
Eine dritte Variante der kanonischen Korrelationsanalyse erhalten wir mit der Wahl
Wx = Σ− 1
2xx
Wy = Σ− 1
2yy (8.29)
85
8 Kanonische Korrelationsanalyse
Das Gleichungssystem (8.14) konnen wir dann umformen zu
Σ− 1
2xx ΣxyΣ
− 12
yy︸ ︷︷ ︸ρxy
Σ− 1
2yy ΣT
xyΣ− 1
2xx︸ ︷︷ ︸
ρTxy
~iα = µ2
i~
iα
Σ− 1
2yy ΣT
xyΣ− 1
2xx Σ
− 12
xx ΣxyΣ− 1
2yy
~iβ = µ2
i~
iβ (8.30)
Dabei ist die Matrix ρxy eine verallgemeinerte Korrelationsmatrix, die nur im Falle dia-
gonaler Kovarianzmatrizen Σxx, Σyy mit der bisher definierten Korrelationsmatrix identisch
sind. Benutzt man die symmetrische Darstellung der kanonischen Korrelationsanalyse, so
sind die Eigenwerte µ2 und Eigenvektoren ~iα,
~iβ durch eine SVD Zerlegung der Matrix ρxy
zu gewinnen.
ρxy = UΛV T
ρTxy = V ΛUT
mit UT U = 1
V T V = 1
ρxyρTxy = UΛV T V ΛUT = UΛ2UT
ρTxyρxy = V ΛUT UΛV T = V Λ2V T (8.31)
Somit sind die Linkseigenvektoren (die Spalten der Matrix U) die gesuchten Eigenvektoren ~iα
und die Rechtseigenvektoren (die Spalten der Matrix V ) die Eigenvektoren~β. Die Eigenwerte
µ2 berechnen sich als die Quadrate der singularen Werte Λ. Durch eine einfache Umformung
der Gleichung (8.14) erhalten wir dann die Beziehung zwischen den kanonischen Mustern
~αi, ~βi bzw. den kanonischen Korrelationsmustern ~iα,
~iβ und den singularen Vektoren ~
iα,~
iβ
zu:
~iα = Σ
12xx~αi
~iβ = Σ
12yy
~βi
~iα = Σ
− 12
xx~
iα
~iβ = Σ
− 12
yy~
iβ (8.32)
Damit kann also auch aus der SVD Zerlegung der Matrix ρxy der vollstandige Satz der
verschiedenen Typen von Mustern der kanonischen Korrelationsanalyse erstellt werden.
86
8 Kanonische Korrelationsanalyse
8.2 Schatzung der kanonischen Korrelation
Seien Dx und Dy die zentrierten Datenmatrizen der beiden Stichproben, aus denen die
kanonische Korrelation geschatzt werden soll. D.h. Dx ist eine q × m Matrix und Dy eine
r × m Matrix, wenn m die Stichprobenlange ist.
Die ML Schatzer bei Annahme eine Normalverteilten GG sind dann
Sxx =1
mDxD
Tx
Syy =1
mDyD
Ty
Sxy =1
mDxD
Ty (8.33)
Die Schatzer fur die verschiedenen kanonischen Eigenvektoren erhalten wir nun, indem wir
die Matrixschatzer in die entsprechenden Gleichungen (8.25), (8.27) oder (8.30) einsetzen.
Erforderlich ist dabei, dass die Kovarianzmatrizen Sxx und Syy den vollen Rang haben, da
die Inversen bzw. die Quadratwurzel der entsprechenden Matrizen bestimmt werden mussen.
Dies erfolgt ublicherweise durch eine a-priori Datenreduktion mit Hilfe vorgegebener guess-
Vektoren. Die Erfahrung zeigt, dass eine erhebliche Datenreduktion i.A. notwendig ist, um
robuste Ergebnisse zu erhalten q, r ∼ 13. . . 1
4m.
8.3 Signifikanzschatzungen der kanonischen Korrelationen
In der Klimatologie sind parametrische Signifikanzschatzungen fur die CCA nicht sehr ver-
breitet. Eher findet man die Versuche, die Nullhypothese durch Monte Carlo Experimente
zu generieren und den Test somit nicht-parametrisch und numerisch durchzufuhren. Hierbei
ergibt sich allerdings immer noch die Frage, welche Form die univariate Testvariable anneh-
men soll. Die Antwort liefern die parametrischen Tests, so dass man dann die Monte Carlo
Experimente gezielt durchrechnen kann und die Ergebnisse garantiert fur grosse Stichproben
und Normalverteilte Variablen gegen die parametrischen Testverfahren konvergiert.
Die Nullhypothese bei der CCA ist die Annahme, dass die beiden multivariaten ZVA ~X
und ~Y unabhangig sind:
Σxy = 0 (8.34)
87
8 Kanonische Korrelationsanalyse
Die gemeinsame Kovarianzmatrix hat dann folgendes Aussehen:
Σ =
Σxx 0
0 Σxx
(8.35)
Die Schatzung der gemeinsamen Kovarianzmatrix S sieht aber so aus:
S =
Sxx Sxy
STxy Syy
(8.36)
Die allgemeine Testvariable erhalt man als sogenanntes likelihood Verhaltnis U
U =det S
det Sxx det Syy
(8.37)
8.4 Beispiel fur eine CCA
Das folgende Beispiel ist der Promotionsarbeit von Petra Friederichs entnommen. Dabei
wurde der (lineare) Zusammenhang von Ozeanoberflachentemperaturen (SST) und atmo-
spharischen Stromungen untersucht. Datengrundlage sind einmal 500 hPa Geopotential Mo-
natsmittel aus NCEP Reanalysen und ECHAM3 Modellsimulationen aus der einen Seite
und beobachtete SST aus dem GISST Datensatz des britischen Hadley Centers auf der
anderen Seite (die ECHAM3 Simulationen sind mit den SST angetrieben worden). Die Da-
tenreduktion erfolgte fur beide Datensatze durch EOF’s. Die Berechnung der kanonischen
Korrelationsmuster ist identisch zur dritten Berechnungsmethode, die oben angegeben wur-
de. In der EOF Darstellung sind die Kovarianzmatrizen Sxx und Syy diagonal. Somit liefert
die SVD der Korrelationsmatrix ρxy die gewunschten Vektoren. Die Anzahl der EOF’s kann
in der CCA Analyse systematisch variiert werden und der hypothesentest somit wieder in
Form einer (doppelten) Hierarchie durchgefuhrt werden (Abb.(14)).
88
8 Kanonische Korrelationsanalyse
Abbildung 12 CCA Analyse zwischen 500 hPa Geopotential und SST im Atlantik, Daten NCEP
Reanalysen
89
8 Kanonische Korrelationsanalyse
Abbildung 13 Wie Abb.(12 jedoch fur die Modellsimulationen ECHAM3
90
8K
anonisch
eK
orrela
tionsa
naly
se
3 6 9 3 6 9 3 6 9 3 6 9 3 6 9 3 6 9 3 6 9 3 6 9 3 6 9 0
0.2
0.4
0.6
0.8
1
no EOF Euro−Atlantic Z500
3 4 5 6 7 8 9 10 11
0
0.2
0.4
0.6
0.8
1
no EOF North−Atlantic SST
3 6 9 3 6 9 3 6 9 3 6 9 3 6 9 3 6 9 3 6 9 3 6 9 3 6 9 −0.2
0
0.2
0.4
0.6
0.8
no EOF Euro−Atlantic Z500
3 4 5 6 7 8 9 10 11
−0.2
0
0.2
0.4
0.6
0.8
no EOF North−Atlantic SST
Abbild
ung
14
Hierarch
ische
Tests
zur
Uberp
rufu
ng
der
Nullh
ypoth
eseΣ
xy
=0
91
Literatur
Literatur
[1] Kolmogoroff A., Grundbegriffe der Wahrscheinlichkeitsrechnung, Berlin, Springer,
1933
[2] Schonwiese, C.D., Praktische Statistik, Gebr. Borntraeger, Berlin 1985
[3] Brandt, S., Datenanalyse, BI Wissenschaftsverlag 1981
[4] Berger J.O. An Introduction to Bayesian Statistics, Springer, 1985
[5] Kreyzig, E., Statistische Methoden und ihre Anwendungen, Vandenhoeck und
Ruprecht, 1975
[6] Taubenheim, J., Statistische Auswertung geopysikalischer und meteorologischer
Daten, Leipzig, Akademische Verlagsgesellschaft, 1979 (wird nicht mehr aufgelegt)
[7] Press, W.H., Flannery, B.P., Teukalsky, S.A., Vetterling W.T., Numerical Recipes,
Cambridge University Press, 1986
[8] Schuster, Deterministic Chaos, An Introduction, Physik - Verlag, Weinheim
[9] Morrison, D.F., Multivariate Statistical Methods, McGraw Hill Series in Proba-
bility and Statistics
[10] Anderson, T.W., An Introduction to Multivariate Statistical Analysis, 2nd Editi-
on, J. Wiley & Sons,
[11] Okamoto M (1963), An asymptotic expansion for the distribution of the linear
discriminant function, Ann.Math.Stat.,34,1286-1301
[12] Lorenz E.N. (1956), Empirical orthogonal functions and statistical weather predic-
tion, Technical report, Statistical Forecast Project Report No.1, Dept. of Meteor.,
MIT, 49pp (in der Bibliothek des Meteorologischen Instituts Universit2at Bonn
verfugbar)
[13] v Storch H. und FW Zwiers (1988) Recurrence analysis of climate sensitivity
experiments, J.Climate, 1, 17-171
92
Literatur
[14] Zwiers FW und H v Storch (1989) Multivariate recurrence analysis, J.Climate,2,
1538-1553
[15] Hans von Storch and Francis Zwiers (1999) Statistical Analysis in Climate Rese-
arch, Cambridge University Press, 483pp
[16] Hartung, J. und Elpelt, B. Multivariate Statistik, Lehr- und Handbuch der ange-
wandten Statistik
[17] Kloeden PE und Platen E (1999) Numerical solution of stochastic differential
equations, Springer, 636pp
93