Skript zur Vorlesung - meteo.uni-bonn.de · Meteorologisches Institut der Universit at Bonn Skript...

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.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


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


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.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


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


-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


(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


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


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


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


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


• .., 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


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


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


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


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


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


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.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


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


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.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


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


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


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.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


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


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


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


Σ. 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


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


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


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


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


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


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


• 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


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


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


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


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


Abbildung 7 Wie Abb.(4.6) jedoch fur den 3-Monatszeitraum Dezember Jahr-0 bis Februar

Jahr+1

51

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


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


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


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


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


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


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


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


-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


-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.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


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


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


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


• 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


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


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


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.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


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


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


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


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.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.


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


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


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


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


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


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


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


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


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.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


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.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


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


Abbildung 12 CCA Analyse zwischen 500 hPa Geopotential und SST im Atlantik, Daten NCEP

Reanalysen

89


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

Skript zur Vorlesung - meteo.uni-bonn.de · Meteorologisches Institut der Universit at Bonn Skript...

Documents

Transcript of Skript zur Vorlesung - meteo.uni-bonn.de · Meteorologisches Institut der Universit at Bonn Skript...