Die Methode der Kleinsten Fehlerquadrate
▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀
Angenommen, es werden bei einem Versuch n Meßwerte yi, i=1,n,
erfaßt. Wir stellen die yi als Spaltenvektor y (ohne Index) dar,
y habe den Erwartungswert E(y) ≡ y^, und auch die Kovarianzmatrix
Cy des Vektors y sei bekannt:
┌──────────────────────────┐
│ ┌ T ┐ │
│ Cy = E│ (y-y^)∙(y-y^) │ │
│ └ ┘ │
└──────────────────────────┘ .
Nehmen wir zunächst an, daß die y normalverteilt seien.
Wir setzen voraus, daß Cy invertierbar ist und benutzen statt der
Kovarianzmatrix Cy besser ihre Inverse, die Gewichtsmatrix P:
-1
Gewichtsmatrix P := Cy .
┌ dW ┐
Dann lautet die Wahrscheinlichkeitsdichtefunktion f(y) │= ───│:
└ dⁿy┘
┌────────────────────────────────────────────┐
│ ┌ ┐½ ┌ T ┐ │
│ │ │P│ │ -½∙│ (y-y^) ∙P∙(y-y^) │ │
│ f(y) = │──────│ ∙ e └ ┘ │
│ │(2∙π)ⁿ│ │
│ └ ┘ │
└────────────────────────────────────────────┘ .
Da y^ im allgemeinen nicht aus den Meßwerten y bestimmbar ist, ist
auch f(y) nicht bekannt.
Es seien nun die gemessenen y in einem linearen Zusammenhang mit
m < n Koeffizienten a derart, daß
┌─────────────┐
│ y = F∙a + v │ , und es gebe ein a^ ≡ E(a), so daß
│ │
│ y^= F∙a^ │
└─────────────┘ ;
die Größe v nennt man "scheinbare Fehler":
┌───────────────────┐
│ v = y - F∙a │
│ scheinbare Fehler │
└───────────────────┘ .
Dann können wir eine Wahrscheinlichkeitsdichtefunktion f(y;a)
einführen, die von den Koeffizienten a abhängt und im Falle a = a^
identisch ist mit f(y):
┌────────────────────────────────────────────────────┐
│ ┌ ┐½ ┌ T ┐ │
│ │ │P│ │ -½∙│ (y - F∙a) ∙P∙(y - F∙a) │ │
│ f(y;a) = │──────│ ∙ e └ ┘ │
│ │(2∙π)ⁿ│ │
│ └ ┘ │
└────────────────────────────────────────────────────┘ ,
was kürzer aus so geschrieben werden kann:
┌────────────────────────────────────┐
│ ┌ ┐½ ┌ T ┐ │
│ │ │P│ │ -½∙│ v ∙P∙v │ │
│ f(y;a) = │──────│ ∙ e └ ┘ │
│ │(2∙π)ⁿ│ │
│ └ ┘ │
└────────────────────────────────────┘ .
Nach dem Prinzip der Maximum Likelihood ist nun a so zu bestimmen,
daß für die gemessenen y die sog. Likelihoodfunktion L, in diesem
Fall also f(y;a), maximal wird. Statt das Maximum von L = f(y;a)
zu bestimmen ist es einfacher, das der sog. logarithmischen
Likelihoodfunktion l = ln(L) zu bestimmen. Die logarithmischen
Likelihoodfunktion l lautet:
┌──────────────────────────────────────────────┐
│ ┌ ┐ │
│ │ │P│ │ ┌ T ┐ │
│ l = ln(f(y;a)) = ½∙ln│──────│ - ½∙│ v ∙P∙v │ │
│ │(2∙π)ⁿ│ └ ┘ │
│ └ ┘ │
└──────────────────────────────────────────────┘ ,
und offenbar gilt:
┌────────────────────────────────────────────────────────┐
│ Die Likelihood L = f(y;a) ist genau dann maximal, wenn │
│ │
│ T │
│ die Fehlerquadratsumme Q := v∙P∙v minimal ist. │
│ │
│ Das ist die Forderung der "kleinsten Fehlerquadrate"! │
└────────────────────────────────────────────────────────┘ .
Bestimmen wir deshalb a so, daß Q minimal ist. Es ist:
T T
Q = v ∙P∙v = (y-F∙a) ∙P∙(y-F∙a)
T T T
= (y - a ∙F )∙P∙(y - F∙a) , also
T T T T T T
Q = y ∙P∙y - y ∙P∙F∙a - a ∙F ∙P∙y + a ∙F ∙P∙F∙a .
Q ist ein Skalar, also auch die einzelnen Glieder des obigen
Ausdrucks. Deshalb gilt:
┌ ┐T
T T │ T T │ T
a ∙F ∙P∙y = │a ∙F ∙P∙y│ = y ∙P∙F∙a , so daß die beiden
└ ┘
mittleren Glieder des obigen Ausdrucks für Q identisch sind und
wir daher auch schreiben können:
T T T T T
Q = y ∙P∙y - 2∙a ∙F ∙P∙y + a ∙F ∙P∙F∙a .
Schreiben wir die partiellen Ableitungen dQ/daj als Spalte auf, so
erhalten wir unter Berücksichtigung, daß
T
F ∙P∙F symmetrisch
ist:
dQ T T
──── = -2∙F ∙P∙y + 2∙F ∙P∙F∙a .
da
Das Minimum erhalten wir durch Nullsetzen, was uns zu den sog.
Normalgleichungen führt:
┌───────────────────┐
│ T T │
│ F ∙P∙F∙a = F ∙P∙y │
└───────────────────┘ .
Formal erhalten wir damit als Schätzung a˜͂ für den Vektor der
unbekannten Koeffizienten:
┌────────────────────────┐
│ ┌ T ┐-1 T │ ┌ ┐
│ a˜͂ = │F ∙P∙F│ ∙F ∙P∙y │ │ = A∙y │
│ └ ┘ │ └ ┘
└────────────────────────┘ .
Offenbar ist die Schätzung a˜͂ erwartungstreu, denn es ist:
┌ ┐
│┌ T ┐-1 T │ ┌ T ┐-1 T
E(a˜͂) = E││F ∙P∙F│ ∙F ∙P∙y│ = │F ∙P∙F│ ∙F ∙P∙E(y)
│└ ┘ │ └ ┘
└ ┘
┌ T ┐-1 T
= │F ∙P∙F│ ∙F ∙P∙y^
└ ┘
┌ T ┐-1 T
= │F ∙P∙F│ ∙F ∙P∙F∙a^ = a^ , also
└ ┘
┌────────────┐
│ E(a˜͂) = a^ │
└────────────┘ .
Bestimmen wir sogleich auch die Kovarianzmatrix Ca˜͂ der Schätzung
a˜͂:
┌ T┐ ┌ T┐
Ca˜͂ = E│(a˜͂-a^)∙(a˜͂-a^) │ = E│(A∙y-a^)∙(A∙y-a^) │
└ ┘ └ ┘
┌ T T T ┐
= E│(A∙y-a^)∙(y∙A -a^ )│
└ ┘
┌ T T T T T T ┐
= E│ A∙y∙y ∙A - A∙y∙a^ - a^∙y ∙A + a^∙a^ │
└ ┘
┌ T┐ T T T T
= A∙E│y∙y │∙A - a^∙a^ - a^∙a^ + a^∙a^
└ ┘
┌ T┐ T T
= A∙E│y∙y │∙A - a^∙a^ .
└ ┘
-1 ┌ T┐ ┌ T┐ T
Aus Cy = P = E│(y-y^)∙(y-y^) │ = E│y∙y │ - y^∙y^ folgt ferner:
└ ┘ └ ┘
┌ T┐ -1 T
E│y∙y │ = P + y^∙y^ .
└ ┘
Setzen wir dieses in den Ausdruck für Ca˜͂ ein, so erhält man:
-1 T T T T
Ca˜͂ = A∙P ∙A + A∙y^∙y^ ∙A - a^∙a^
┌ T ┐-1 T -1 ┌ T ┐-1 T T
= │F ∙P∙F│ ∙F ∙P∙P ∙P∙F∙│F ∙P∙F│ + a^∙a^ - a^∙a^ ,
└ ┘ └ ┘
┌──────────────────┐
│ ┌ T ┐-1 │
│ Ca˜͂ = │F ∙P∙F│ │
│ └ ┘ │
└──────────────────┘ .
Entsprechend erhalten wir Aussagen zur Schätzung der Meßwerte y˜͂:
┌───────────┐
│ y˜͂ = F∙a˜͂ │
└───────────┘
und deren Kovarianzmatrix Cy˜͂:
E(y˜͂) = F∙E(a˜͂) = F∙a^ = y^ ,
┌──────────────────────────────────┐
│ T ┌ T ┐-1 T │
│ Cy˜͂ = F∙Ca˜͂∙F = F∙│F ∙P∙F│ ∙F │
│ └ ┘ │
└──────────────────────────────────┘ .
T
Eine wichtige im folgenden oft gebrauchte Größe ist Q = v ∙P∙v.
Betrachten wir zunächst den Erwartungswert E(Q). Für die
Herleitung werden wir die folgenden Eigenschaften des
Spuroperators sp benutzen:
1. Unter dem Spuroperator sp verhalten sich Matrizen kommutativ:
┌───────────────────┐
│ sp(A∙B) = sp(B∙A) │ , da:
└───────────────────┘
n m m n
sp(A∙B) = Σ Σ aik∙bki = Σ Σ bki∙aik = sp(B∙A) .
i=1 k=1 k=1 i=1
2. Sei B eine beliebige quadratische n x n -Matrix und ε ein
n-dimensionaler Zufallsvektor mit
T
Erwartungswert E(ε) = 0 und Kovarianzmatrix Cε = E(ε∙ε ).
Dann gilt
┌──────────────────────┐
│ T │
│ E(ε ∙B∙ε) = sp(B∙Cε) │ , denn es ist
└──────────────────────┘
T T T
E(ε ∙B∙ε) = E(sp(ε ∙B∙ε)) da ε ∙B∙ε ein Skalar,
T
= E(sp(B∙ε∙ε )) da kommutativ unter sp,
T
= sp(E(B∙ε∙ε )) da E(Summe) = Summe(E),
T
= sp(B∙E(ε∙ε )) da B konstant ist, also:
= sp(B∙Cε) .
Erinnern wir uns, daß galt:
-1
Y = Y^ + ε mit ε = 'wahrer' Fehler mit Cε = P bekannt,
Y = Y˜͂ + v mit v = scheinbare Fehler aus Ausgleichung,
┌ T ┐-1 T
Y˜͂ = F∙a˜͂ = F∙│F ∙P∙F│ ∙F ∙P∙Y als Ausgleichslösung.
└ ┘
Schreiben wir für die n- bzw. m- dim. Einheitsmatrix In bzw. Im,
so ist:
┌ T ┐-1 T
v = Y - Y˜͂ = (In - F∙│F ∙P∙F│ ∙F ∙P)∙Y , und mit Y = Y^ + ε:
└ ┘
┌ T ┐-1 T ┌ T ┐-1 T
v = Y^ - F∙│F ∙P∙F│ ∙F ∙P∙Y^ + (In - F∙│F ∙P∙F│ ∙F ∙P)∙ε ,
└ ┘ └ ┘
┌ T ┐-1 T ┌ T ┐-1 T
v = Y^ - F∙│F ∙P∙F│ ∙F ∙P∙F∙a^ + (In - F∙│F ∙P∙F│ ∙F ∙P)∙ε ,
└ ┘ └ ┘
┌ T ┐-1 T
v = Y^ - F∙a^ + (In - F∙│F ∙P∙F│ ∙F ∙P)∙ε ,
└ ┘
┌ T ┐-1 T
v = Y^ - Y^ + (In - F∙│F ∙P∙F│ ∙F ∙P)∙ε ,
└ ┘
┌ T ┐-1 T
v = (In - F∙│F ∙P∙F│ ∙F ∙P)∙ε .
└ ┘
Nun bilden wir:
T T ┌ T ┐-1 T ┌ T ┐-1 T
v ∙P∙v = ε ∙(In-P∙F∙│F ∙P∙F│ ∙F )∙P∙(In-F∙│F ∙P∙F│ ∙F ∙P)∙ε
└ ┘ └ ┘
T ┌ T ┐-1 T
= ε ∙(P - P∙F∙│F ∙P∙F│ ∙F ∙P)∙ε
└ ┘
T T ┌ T ┐-1 T
= ε ∙P∙ε - ε ∙P∙F∙│F ∙P∙F│ ∙F ∙P∙ε .
└ ┘
T
Jetzt können wir den Erwartungswert E(v ∙P∙v) unter Benutzung der
zuvor erwähnten Eigenschaften des Spuroperators bestimmen:
T T
E(v ∙P∙v) = E(sp(v ∙P∙v))
T T ┌ T ┐-1 T
= E(sp(ε ∙P∙ε - ε ∙P∙F∙│F ∙P∙F│ ∙F ∙P∙ε))
└ ┘
T T ┌ T ┐-1 T
= E(sp(ε ∙P∙ε)) - E(sp(ε ∙P∙F∙│F ∙P∙F│ ∙F ∙P∙ε))
└ ┘
T ┌ T ┐-1 T T
= E(sp(P∙ε∙ε )) - E(sp(P∙F∙│F ∙P∙F│ ∙F ∙P∙ε∙ε ))
└ ┘
T ┌ T ┐-1 T T
= sp(E(P∙ε∙ε )) - sp(E(P∙F∙│F ∙P∙F│ ∙F ∙P∙ε∙ε ))
└ ┘
T ┌ T ┐-1 T T
= sp(P∙E(ε∙ε )) - sp(P∙F∙│F ∙P∙F│ ∙F ∙P∙E(ε∙ε ))
└ ┘
┌ T ┐-1 T
= sp(In) - sp(P∙F∙│F ∙P∙F│ ∙F )
└ ┘
┌ T ┐-1 T
= n - sp(│F ∙P∙F│ ∙F ∙P∙F)
└ ┘
= n - sp(Im) = n - m , also
┌───────────────────────────┐
│ T │
│ E(Q) = E(v ∙P∙v) = n - m │
└───────────────────────────┘ .
* * *
Zur Ableitung einer möglichst optimalen Methode zur Bestimmung von
Parametern aus Messungen gingen wir von der Maximum-Likelihood-
Methode aus, die uns unter der Annahme gaußverteilter Meßwerte zur
Methode der kleinsten Fehlerquadrate führte.
Die Methode der kleinsten Fehlerquadrate ist aber auch anwendbar,
wenn die Meßwerte nicht gaußverteilt sind! Voraussetzung ist nur,
daß die Varianz der Verteilung existiert. In aller Regel ist das
der Fall. Dennoch sei als warnendes Beispiel an die
Cauchyverteilung erinnert, deren Varianz nicht existiert.
Es bleibt noch zu klären, ob es eventuell eine andere Schätzung
als die Methode der kleinsten Fehlerquadrate gibt, die kleinere
Varianzen liefert.
Wir werden jetzt zeigen, daß die oben hergeleiteten Kovarianz-
matrizen Ca˜͂ und Cy˜͂ die Schätzungen kleinster Varianz im linearen
Modell (y^ = F∙a^) darstellen. Wir führen die Herleitung für Ca˜͂
(Fall a) und die für Cy˜͂ (Fall b) gleichzeitig durch.
Behauptung: Die lineare, erwartungstreue Schätzung mit kleinster
Varianz lautet:
a: ┌ T ┐-1 T
für a^ : a˜͂ = │F ∙P∙F│ ∙F ∙P∙y,
└ ┘
b: ┌ T ┐-1 T
für y^ : y˜͂ = F∙│F ∙P∙F│ ∙F ∙P∙y .
└ ┘
Beweis durch Widerspruch: Annahme, es existiere eine weitere
lineare erwartungstreue Schätzung mit kleinerer Varianz:
a: t˜͂ = L∙y mit
a^ = E(t˜͂) = L∙y^ = L∙F∙a^
Da a^ = L∙F∙a^ und a^ den vollen Zeilenrang hat, folgt:
L∙F = I (Einheitsmatrix)
b: z˜͂ = M∙y mit
y^ = E(z˜͂) = M∙y^ = M∙F∙a^ .
Da andererseits auch y^ = F∙a^ gilt, so folgt:
M∙F∙a^ = F∙a^ .
Da a^ den vollen Zeilenrang hat, folgt daher:
M∙F = F .
Mit diesem Ergebnis kann man zeigen, daß gilt:
a: ┌ ┐ ┌ ┐T
-1 T │┌ T ┐-1 T │ -1 │┌ T ┐-1 T │
Ct˜͂ = L∙P ∙L = ││F ∙P∙F│ ∙F ∙P│∙P ∙││F ∙P∙F│ ∙F ∙P│
│└ ┘ │ │└ ┘ │
└ ┘ └ ┘
┌ ┐ ┌ ┐T
│ ┌ T ┐-1 T │ -1 │ ┌ T ┐-1 T │
+ │L - │F ∙P∙F│ ∙F ∙P│∙P ∙│L - │F ∙P∙F│ ∙F ∙P│ ,
│ └ ┘ │ │ └ ┘ │
└ ┘ └ ┘
Ct˜͂ = R + S .
b: ┌ ┐ ┌ ┐T
-1 T │ ┌ T ┐-1 T │ -1 │ ┌ T ┐-1 T │
Cz˜͂ = M∙P ∙M = │F∙│F ∙P∙F│ ∙F ∙P│∙P ∙│F∙│F ∙P∙F│ ∙F ∙P│
│ └ ┘ │ │ └ ┘ │
└ ┘ └ ┘
┌ ┐ ┌ ┐T
│ ┌ T ┐-1 T │ -1 │ ┌ T ┐-1 T │
+ │M - F∙│F ∙P∙F│ ∙F ∙P│∙P ∙│M - F∙│F ∙P∙F│ ∙F ∙P│ ,
│ └ ┘ │ │ └ ┘ │
└ ┘ └ ┘
Cz˜͂ = U + W .
Da die Diagonalelemente der Matrizen R,S bzw. U,W immer größer
oder gleich Null sind, können die Diagonalelemente von Ct˜͂ bzw.
Cz˜͂ nur minimal sein (und damit die Varianzen), wenn
a: ┌ T ┐-1 T
L = │F ∙P∙F│ ∙F ∙P ,
└ ┘
b: ┌ T ┐-1 T
M = F∙│F ∙P∙F│ ∙F ∙P .
└ ┘
Das ist aber ein Widerspruch zur Annahme, womit obige Behauptung
bewiesen ist.
* * *
Wir wollen noch einige nützliche Beziehungen zwischen den y, y˜͂, v
sowie den zugehörigen Kovarianzmatrizen Cy, Cy˜͂, Cv herleiten.
Gehen wir aus von der Gleichung:
┌ T ┐-1 T
y˜͂= F∙│F ∙P∙F│ ∙F ∙P∙y = H∙P∙y und
└ ┘
┌ ┐
│ ┌ T ┐-1 T │ ┌ ┐
v = y - y˜͂ = │I - F∙│F ∙P∙F│ ∙F ∙P│∙y = │I - H∙P│∙y .
│ └ ┘ │ └ ┘
└ ┘
Es gelten die Beziehungen:
┌ ┐ ┌ ┐ ┌ ┐ ┌ ┐ ┌ ┐ ┌ ┐
│H∙P│∙│H∙P│ = │H∙P│ und │I - H∙P│∙│I - H∙P│ = │I - H∙P│ .
└ ┘ └ ┘ └ ┘ └ ┘ └ ┘ └ ┘
Wir können nun schreiben:
T T ┌ ┐ ┌ ┐
v ∙P∙v = y ∙│I - P∙H│∙P∙│I - H∙P│∙y
└ ┘ └ ┘
T T T T
= y ∙P∙y - y ∙P∙H∙P∙y - y ∙P∙H∙P∙y + y ∙P∙H∙P∙H∙P∙y
T T
= y ∙P∙y - y ∙P∙H∙P∙y
T T
= y ∙P∙y - y ∙P∙H∙P∙H∙P∙y
T ┌ ┐T
= y ∙P∙y - │H∙P∙y│ ∙P∙H∙P∙y ,
└ ┘
T T T
v ∙P∙v = y ∙P∙y - y˜͂ ∙P∙y˜͂ oder
┌───────────────────────────────────────────────────────────┐
│ T T T T T T │
│ y ∙P∙y = y˜͂ ∙P∙y˜͂ + v ∙P∙v ( = a˜͂ ∙F ∙P∙F∙a˜͂ + v ∙P∙v ) │
└───────────────────────────────────────────────────────────┘.
Zur Herleitung einer Beziehung zwischen den Kovarianzmatrizen
gehen wir aus von:
T
y = y^ + ε mit Cy ≡ Cε ≡ E(ε∙ε ), E(ε) = 0 .
┌ T┐ ┌ T T T T┐
E│y∙y │ = E│y^∙y^ + y^∙ε + ε∙y^ + ε∙ε │ ,
└ ┘ └ ┘
┌ T┐ T T
E│y∙y │ = y^∙y^ + Cε = y^∙y^ + Cy .
└ ┘
Betrachten wir nun:
-1 T
y˜͂ = H∙P∙y mit P = Cy und H = H .
Wir schreiben:
┌ ┐
┌ T┐ │┌ ┐ T ┌ ┐│
Cv = E│v∙v │ = E││I - H∙P│∙y∙y ∙│I - P∙H││
└ ┘ │└ ┘ └ ┘│
└ ┘
┌ ┐ ┌ T┐ ┌ ┐
= │I - H∙P│∙E│y∙y │∙│I - P∙H│
└ ┘ └ ┘ └ ┘
┌ ┐ ┌ T ┐ ┌ ┐
= │I - H∙P│∙│y^∙y^ + Cy│∙│I - P∙H│
└ ┘ └ ┘ └ ┘
┌ ┐ T ┌ ┐ ┌ ┐ ┌ ┐
= │I - H∙P│∙y^∙y^ ∙│I - P∙H│ + │I - H∙P│∙Cy∙│I - P∙H│ ,
└ ┘ └ ┘ └ ┘ └ ┘
und da (siehe Vorseite oben):
┌ ┐ ┌ T ┐-1 T
│I - H∙P│∙y^ = y^ - F∙│F ∙P∙F│ ∙F ∙P∙y^ = Y^ - Y^ = 0
└ ┘ └ ┘
gilt, erhalten wir:
Cv = Cy - H - H + H∙P∙H = Cy - H (H∙P∙H = H ≡ Cy˜͂),
Cv = Cy - Cy˜͂ oder
┌───────────────┐
│ Cy = Cy˜͂ + Cv │
└───────────────┘ .
* * *
Zur Schätzung der Kovarianzmatrix C˜͂y der Beobachtungen y. Nehmen
wir zunächst an, es sei statt der wahren Kovarianzmatrix Cy nur
eine Matrix C'y (mit Cy = σ²∙C'y) bekannt. Die Kovarianzmatrix Cy
ist also bis auf den skalaren Faktor σ² bekannt. Besteht die
Möglichkeit σ² zu schätzen und damit auch C˜͂y = σ²∙C'y als
Schätzung für Cy zu nehmen?
Wenn Cy = σ²∙C'y, so haben wir eine C'y entsprechende Gewichts-
matrix P':
P' = σ²∙P .
Zunächst zeigt sich, daß die Formeln zur Bestimmung von a˜͂ und y˜͂
invariant gegen einen skalaren Faktor von P sind. Das gilt aber
nicht für die Formeln für Ca˜͂, Cy˜͂ und Q! Jedoch finden wir bei
Betrachtung der Größe E(Q'):
┌ T ┐
E(Q') = E│v ∙P'∙v│
└ ┘
┌ T ┐
= E│v ∙σ²∙P∙v│
└ ┘
┌ T ┐
= σ²∙E│v ∙P∙v│
└ ┘
= σ²∙E(Q) ,
E(Q') = σ²∙(n-m) ,
so daß wir σ²˜͂ als Schätzung für σ² offenbar nehmen können:
┌─────────────────────────┐
│ T │
│ Q' v ∙P'∙v │
│ σ²˜͂ = ───── = ───────── │
│ n - m n - m │
└─────────────────────────┘ .
Damit haben wir aber auch erwartungstreue Schätzungen C˜͂y, C˜͂a˜͂
und C˜͂y˜͂ für die Kovarianzmatrizen Cy, Ca˜͂ und Cy˜͂, nämlich:
┌────────────────────────────────────────────────────────┐
│ C˜͂y = σ²˜͂∙C'y │
│ │
│ ┌ T ┐-1 │
│ C˜͂a˜͂ = σ²˜͂∙C'a˜͂ = σ²˜͂∙│F ∙P'∙F│ │
│ └ ┘ │
│ │
│ T ┌ T ┐-1 T │
│ C˜͂y˜͂ = σ²˜͂∙C'y˜͂ = σ²˜͂∙F∙C'a˜͂∙F = σ²˜͂∙F∙│F ∙P'∙F│ ∙F │
│ └ ┘ │
└────────────────────────────────────────────────────────┘ .
Betrachten wir:
Cy = σ²∙C'y ,
und wählen wir C'y so, daß ein Diagonalelement, das ja die Varianz
eines Beobachtungswertes von y darstellt, gleich 1 ist, so stellt
σ² offensichtlich die Schätzung der Varianz für diesen Datenwert
dar.
In der Praxis ist es meist so, daß Cy diagonal und damit auch P
diagonal ist. Dann kann man P' so skalieren, daß das geometrische
Mittel der Diagonalelemente gleich 1 ist. Dann stellt σ² die
Schätzung der Varianz für einen Datenwert dar, der die mittlere
Varianz besitzt.
* * *
Die bisherige Betrachtung der "Varianzschätzung" macht keine
Voraussetzung über die Wahrscheinlichkeitsdichte der Messungen y.
Daher sind auch keine Aussagen über die Zuverlässigkeit der
Varianzschätzungen möglich.
Gehorchen die Meßwerte jedoch einer Gaußverteilung, so gilt (ohne
Beweis) folgendes:
T
1) Q = v ∙P∙v gehorcht einer X²-Verteilung mit n-m Freiheits-
graden, d.h., die Wahrscheinlichkeitsdichte lautet:
1 [(n-m)/2 - 1] -Q/2
f(Q) = ───────────────────────∙Q ∙e .
[(n-m)/2]
Γ[(n-m)/2]∙2
2) E(Q) = n - m (gilt auch sonst) ,
3) die Varianz von Q ist:
var(Q) = 2∙(n-m) ,
4) die Streuung von Q ist also:
½
σ(Q) = (2∙(n-m)) ,
5) für die Streuung der Schätzung der Streuung σ˜͂ (= √(σ²˜͂))
gilt:
┌ T ┐½ ┌ ┐
│v ∙P'∙v│ │ -½ │
√(σ²˜͂) = │───────│ ∙ │ 1 ± (2∙(n-m)) │
│ n - m │ │ │
└ ┘ └ ┘ .
* * *
Befassen wir uns nun mit der Schätzung der Streuung von Meßgruppen
innerhalb einer Messung.
Die Kovarianzmatrix Cy bestehe nun aus g diagonal angeordneten
Untermatrizen Cyi, i=1,g, die ihrerseits die Kovarianzmatrizen von
g Meßgruppen sind. Es sei also:
┌┌───┬───┬─ ┬───┐┐
││Cy1│ * │ * ∙ *│ * ││
│├───┼───┼─ ┼───┤│
││ * │Cy2│ * ∙ *│ * ││
│├───┼───┼─ ┼───┤│
Cy = │ * * ∙ * │
│ ∙ ∙ ∙ ∙ │
│ * * ∙ * │
│├───┼───┼─ ┼───┤│
││ * │ * │ * ∙ *│Cyg││
└└───┴───┴─ ┴───┘┘ , und es gilt außerdem:
Cy = Cy˜͂ + Cv und für die Untermatrizen entsprechend:
Cyi = Cy˜͂i + Cvi , i = 1,g ,
dabei möge die Gruppe i aus ni Messungen bestehen.
Wir führen willkürlich normierte Kovarianzmatrizen ein:
Cyo und Cyoi, i=1,g,
und damit zugleich willkürlich normierte Gewichtsmatrizen:
-1 -1
Po = Cyo und Poi = Cyoi .
Dabei muß aber beachtet werden, daß die Korrelationen der
Gesamtmatrix Cyo erhalten bleiben.
Der Unterschied zwischen "wahren" und "normierten" Kovarianz-
matrizen werde beschrieben durch:
Cy = σo²∙Cyo und Cyi = σo²i∙Cyoi,
daher gilt für die Gewichtsmatrizen entsprechend:
Po = σo²∙P und Poi = σo²i∙Pi .
Von besonderer Bedeutung ist im folgenden die Betrachtung der
Größen:
T T
Q = v ∙P∙v und Qi = vi ∙Pi∙vi, i = 1,g bzw.:
T T
Q' = v ∙Po∙v und Qi' = vi ∙Poi∙vi, i = 1,g .
Der Erwartungswert für Qi ergibt, da wir es mit Diagonalmatrizen
zu tun haben:
T T
E(Qi') = E(vi ∙Poi∙vi) = E(sp(vi ∙Poi∙vi))
T
= E(sp(Poi∙vi∙vi )) = sp(Poi∙Cvi) = σo²i∙sp(Pi∙Cvi)
= σo²i∙sp(Pi∙Cyi - Pi∙Cy˜͂i) = σo²i∙(ni - sp(Pi∙Cy˜͂i)) ,
so daß wir schließlich erhalten:
T
E(vi ∙Poi∙vi)
σo²i = ────────────────── bzw.
ni - sp(Cy˜͂i∙Pi)
T
E(vi ∙Poi∙vi)
σo²i = ──────────────────── , da Cy˜͂∙P = Cy˜͂o∙Po und
ni - sp(Cy˜͂oi∙Poi)
┌ T ┐-1 T
Cy˜͂o = F∙│F ∙Po∙F│ ∙F .
└ ┘
Ersetzen wir den Erwartungswert durch die Summe, so erhalten wir
eine Schätzung für σo²i:
┌──────────────────────────────┐
│ T │
│ vi ∙Poi∙vi │
│ σo²i˜͂ = ──────────────────── │
│ ni - sp(Cy˜͂oi∙Poi) │
└──────────────────────────────┘ .
Diese Überlegung gilt nur, wenn die für alle Meßwerte zuständige
Gewichtsmatrix P bis auf einen Faktor unbekannt ist. Sind jedoch
die Untermatrizen Pi jeweils um untereinander verschiedene
Faktoren falsch, so stellen die σo²i˜͂ Faktoren dar, mit denen die
Gewichtsmatrizen Poi korrigiert werden können. Im Falle, daß die
Messungen nur eine einzige Gruppe bilden, liefert die obige Formel
exakt die Schätzung σo²i˜͂ ≡ σo²˜͂. Im Falle mehrerer Gruppen
erhält man wegen der relativ zueinander inkorrekten
Gruppengewichte nur eine Näherung. Wir müssen daher einen
iterativen Prozeß zur Schätzung der Gruppengewichtsmatrizen bzw.
der mittleren Streuung der Gruppen anwenden.
Für den Fall, daß Cy eine reine Diagonalmatrix ist, sieht die
Iteration folgendermaßen aus:
┌─────────────────────────────────────────────────┐
Anfang: │ 0. τ sei der Iterationsindex (anfangs τ=0). │
│ Die diagonale Gewichtsmatrix Pτ aller Mes- │
│ sungen wird in g Gruppengewichtsmatrizen Piτ │
│ aufgeteilt und mittlere Gruppenvarianzen │
│ σ²i˜͂τ werden berechnet: │
│ -1 │
│ ┌ ni ┐── │
│ Piτ, i=1,g ; σ²i˜͂τ = │ π Pi[jj]τ│ni . │
│ └j=1 ┘ │
└────────────────────────┬────────────────────────┘
│ <--«
├───────────────────────────┐
│ │
┌────────────────────────┴────────────────────────┐ │
Iteration: │ 1. Mit Pτ wird eine Ausgleichung gerechnet. Man │ │
│ erhält u.a. die scheinbaren Fehler aller │ │
│ Messungen und relativ zu Pτ die Kovarianz- │ │
│ matrix der ausgeglichenen Meßwerte (ent- │ │
│ sprechende Ergebnisse auch für die Gruppen): │ │
│ │ │
│ vτ --> viτ, i=1,g und Cy˜͂τ --> Cyi˜͂τ . │ │
│ │ │
│ 2. Schätzung neuer Gruppenvarianzen, Gruppen- │ │
│ gewichtsmatrizen und der Gewichtsmatrix │ │
│ aller Messungen: │ │
│ │ │
│ T │ │
│ vi ∙Piτ∙vi │ │
│ σ²i˜͂(τ+1) = ──────────────────── , i=1,g , │ │
│ ni - sp(Cyi˜͂τ∙Piτ) │ │
│ │ │
│ Pi(τ+1) = Piτ/σ²i˜͂(τ+1), i=1,g , -> P(τ+1) │ │
└────────────────────────┬────────────────────────┘ │
┌───────────────────────┴───────────────────────┐ │
│ g │ σ²i˜͂(τ+1) - σ²i˜͂τ │ │»->│
│ Σ │ ────────────────── │ > Schranke ? ├───┘
│ i=1│ σ²i˜͂(τ+1) │ │ ja
└───────────────────────┬───────────────────────┘
│ nein
┌──────────┴──────────┐
│ STOP der Iteration │
└─────────────────────┘
Kleinste Fehlerquadrate im nichtlinearen Fall
▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀
Analog zum bisher ausschließlich betrachteten linearen Zusammen-
hang zwischen Meßwerten y und Koeffizienten a mit gegebener
Gewichtsmatrix P haben wir nun einen funktionellen Zusammenhang
y=y(x,a), so daß mindestens ein Koeffizient nichtlinear vorkommt.
Es seien also die gemessenen yi; i=1,n in einem funktionellen
Zusammenhang mit m < n Koeffizienten aj; j=1,m derart, daß mit
bekannter Gewichtsmatrix P der y gilt:
┌───────────────────┐
│ yi = y(xi,a) + vi │ , und es gebe ein a^ ≡ E(a), so daß
│ │
│ yi^= y(xi,a^) │
└───────────────────┘ .
Wie im linearen Fall früher beschrieben fordern wir nun, daß
T T n ┌ ┐2
Q = v ∙P∙v = (y-y(a)) ∙P∙(y-y(a)) = Σ Pi∙│yi-y(xi,a)│ min!
i=1 └ ┘
Folgende Vektor- bzw. Matrix-Notationen wollen wir benutzen:
1) y = [yi] Spaltenvektor, n Meßwerte
2) x = [xi] Spaltenvektor, n Meßorte
3) P = [Pi] n x n -Matrix, n Gewichte (diagonal)
4) a = [aj] Spaltenvektor, m unbekannte Koeffizienten
5) a0 = [a0j] Spaltenvektor, Ort für Entwicklung nach a
6) ⌂a = a - a0 Spaltenvektor, Differenz zu a0
7) y(a) = [y(xi,a)] Spaltenvektor, berechnete Werte bei a
8) y(a0) = [y(xi,a0)] Spaltenvektor, berechnete Werte bei a0
9) v = y - y(a) Spaltenvektor, scheinbare Fehler
A) z = y - y(a0) Spaltenvektor, Hilfsvektor
┌ ┐ ┌ ┐
┌ ┐ │dy(xi,a)│ │dy(xi,a0)│
B) │Fij│ = │────────│ bzw. │─────────│ n x m -Matrix
└ ┘ │ daj │ │ daj │
└ ┘ └ ┘
┌ ┐
┌ ┐ │ dQ│
C) │Lj│ = │───│ Spaltenvektor der Ableitungen nach aj
└ ┘ │daj│
└ ┘
┌ ┐
┌ ┐ │ d²Q │
D) │Djk│ = │───────│ m x m -Matrix der 2. part. Abl. von Q
└ ┘ │daj∙dak│
└ ┘
Wir bilden zunächst die ersten und zweiten Ableitungen von Q nach
den aj an der Stelle a0:
dQ n ┌ ┐
─── = Σ -2∙Pi∙│yi-y(xi,a0)│∙Fij ,
daj i=1 └ ┘
T T
L = -2∙F ∙P∙(y-y(a0)) = -2∙F ∙P∙z , und
d²Q n n ┌ ┐ d²y(xi,a)
─────── = Σ 2∙Pi∙Fij∙Fik - Σ 2∙Pi∙│yi-y(xi,a0)│∙───────── ,
daj∙dak i=1 i=1 └ ┘ daj∙dak
T
D = 2∙F ∙P∙F mit Vernachlässigung des zweiten Terms.
Den zweiten Term in der zweiten Ableitung können wir vernach-
lässigen, weil er wegen kleiner vi = yi - y(xi,a) mit statistisch
wechselnden Vorzeichen bei der Summenbildung die Tendenz zum
Verschwinden hat.
Zur Auffindung des Minimums werden wir Q = Q(a) bis zur zweiten
Ordnung an der Stelle a0 entwickeln:
m dQ m,m d²Q
Q = Q(a0) + Σ ───∙⌂aj + ½∙ Σ ⌂aj∙⌂ak∙─────── + Rest ,
j=1 daj j=1 daj∙dak
k=1
T T T
Q = Q(a0) + ⌂a ∙L + ⌂a ∙F ∙P∙F∙⌂a
mit Vernachlässigung des Rests.
Als notwendige Bedingung für das Minimum von Q muß die erste
Ableitung nach den aj als Spalte aufgeschrieben gleich Null sein:
dQ T T T
── = L + 2∙F ∙P∙F∙⌂a = -2∙F ∙P∙z + 2∙F ∙P∙F∙⌂a = 0 , also:
da
┌─────────────────────┐
│ T T │
│ F ∙P∙F∙⌂a = F ∙P∙z │
└─────────────────────┘ .
Wir haben also wieder die Normalengleichungen, die uns schon vom
linearen Fall der Ausgleichung bekannt sind. Man kann nun dieses
Gleichungssystem (für ein vorgeschätztes) a0 lösen und ersetzt
dann a0 durch a0 + ⌂a um die Rechnung daraufhin solange zu
wiederholen, bis ⌂a unter eine gewisse Schranke fällt.
Da naturgemäß das Minimum von zweiter Ordnung ist und zur Lösung
bis zur zweiten Ordnung entwickelt wurde, konvergiert dieses
Verfahren sehr schnell, falls man sich hinreichend nahe am Minimum
befindet und die Koeffizienten a nicht korreliert sind. Es gibt
aber leider auch den Fall, daß man so weit vom Minimum liegt, daß
die Entwicklung der Funktion Q bis zur zweiten Ordnung
unzureichend ist, oder daß die Koeffizienten stark korreliert
sind. In solchen Fällen hat sich das Verfahren nach Marquard
bewährt.
Offenbar ist -L der Gradient von Q an der Stelle a0. Man nähert
sich dem Minimum, wenn man einen hinreichend kleinen Schritt von
a0 in Richtung des negativen Gradienten, also in Richtung
T
F ∙P∙z
bewegt.
Ändert man die obigen Normalengleichungen folgendermaßen ab:
┌ T ┐ T
│F ∙P∙F + I∙δ│∙⌂a = F ∙P∙z ,
└ ┘
mit δ≥0 und I der Einheitsmatrix, so erhält man für δ=0 wieder das
ursprüngliche System. Macht man jedoch δ hinreichend groß, so
wird die Matrix I∙δ dominant und man erhält als Lösung für ⌂a:
1 T
⌂a = ─∙F ∙P∙z ,
δ
was einen Schritt in Richtung des negativen Gradienten und damit
in Richtung Minimum bedeutet. Dabei geht man nun so vor:
wann immer man feststellt, daß Q abgenommen hat, verringert man δ
um den Faktor 10, und im Falle, daß Q gleich blieb oder zunahm,
vergrößert man δ um den Faktor 10.
Der Iterationsablauf kann also aussehen wie auf folgender Seite
dargestellt: