Lineare Ausgleichsprobleme

Datenreduktion und Parameteranpassung

 vexfit.png 

Sowohl bei der Auswertung von Laborexperimenten als auch bei der Analyse von numerischen Berechnungen ergibt sich oft die Situation, daß man für eine Reihe von Werten x1, x2, ..., xN einer unabhängigen Variablen x zugehörige Werte y1, y2, ..., yN einer abhängigen Variablen y gemessen (oder berechnet) hat und an diese Daten eine analytische Funktion y(x) anpassen will.

Der Zweck eines solchen ``Fits'' kann es z.B. sein, die Streuung von mit Meßfehlern behafteten Daten auszugleichen, die Meßwerte für beliebige x zwischen den xi zu interpolieren, oder einen umfangreichen Datensatz auf einige wenige Parameter zu reduzieren, die zur Beschreibung der Funktion y(x) notwendig sind. In all diesen Fällen muß die Fit-Funktion zunächst keine physikalische Bedeutung haben, und man wird für y(x) die einfachste mathematische Form wählen, die es erlaubt, die Daten im Rahmen der Meßfehler zu reproduzieren. Beliebte Funktionen dieser Art sind Polynome, Exponentialfunktionen, rationale Funktionen, usw.

Andererseits kann es sein, daß man, auf der Grundlage eines konkreten physikalischen Modells, sehr wohl eine Vorstellung von der funktionalen Form von y(x) hat. Meist hängt dieses Modell von einer Reihe von unbekannten Parametern ab, und die Aufgabe des Fits ist es, die Werte dieser Parameter aus den Messungen zu bestimmen. Zusätzlich möchte man auch noch beurteilen können, ob das theoretische Modell überhaupt mit den Daten verträglich ist. Es muß also ein Gütekriterium für den Fit formuliert werden.

Beide Fragestellungen sind insofern äquivalent, als sie im allgemeinen mit demselben Verfahren gelöst werden, nämlich der sogenannten ``Methode der kleinsten Quadrate''. Während sich dieses Verfahren im Fall der Parameteranpassung unmittelbar aus den Vorstellungen der Maximum Likelihood-Schätzung ableiten läßt, könnte man bei bloßen Datenreduktionsproblemen aber natürlich auch nach alternativen Strategien vorgehen.


Maximum Likelihood-Schätzung

 mle.png 

Es seien N Paare von Meßwerten mit zugehörigen Meßfehlern


{xi,yi,si}i=1,, N
gegeben. Dabei wird angenommen, daß die xi exakt bekannt, die yi aber mit einem statistischen Meßfehler si behaftet sind. (Die Größenordnung von si kann man z.B. abschätzen, indem man für jedes xi mehrere Messungen von y durchführt und den Standardfehler des Mittelwerts bestimmt.) Weiters wird angenommen, daß die Meßwerte yi statistisch unabhängig und normalverteilt sind, d.h. die Wahrscheinlichkeit (eigentlich Wahrscheinlichkeitsdichte), als i-ten Meßpunkt den Wert yi zu erhalten, ist


p(yi) =

1
2psi2


1/2

 
exp

- [yi-y(xi)]2
2si2


Es wird hier also die Gültigkeit eines ``theoretischen Modells'' y=y(x) angenommen, d.h. y(xi) wäre der ``wahre'' Funktionswert an der Stelle xi, bei einer tatsächlichen Messung wird aber ein davon abweichender Wert yi beobachtet. Zusätzlich wird noch vorausgesetzt, daß die si unabhängig von den im konkreten Fall gemessenen Werten yi sein sollen.

Das theoretische Modell hängt im allgemeinen von M Parametern l1, l2, ..., lM ab


y = y(l1,,lM,x)
Nach dem Prinzip der Maximum Likelihood-Schätzung werden die la so bestimmt, daß die Wahrscheinlichkeit, bei einer Messung genau jene Werte zu erhalten, die man tatsächlich gemessen hat, maximal wird.

Da vorausgesetzt wurde, daß die yi statistisch unabhängig sein sollen, ist die Wahrscheinlichkeitsdichte für das Auftreten der Meßwerte y1, y2, ..., yN das Produkt der Einzelwahrscheinlichkeitsdichten


p(y1,y2,,yN)
=


1
2ps12


1/2

 
exp

- [y1-y(x1)]2
2s12


 

1
2ps22


1/2

 
exp

- [y2-y(x2)]2
2s22


 
×

1
2psN2


1/2

 
exp

- [yN-y(xN)]2
2sN2


=
N

i=1 


1
2psi2


1/2

 
exp

- N

i=1 
[yi-y(xi)]2
2si2


Diese Wahrscheinlichkeitsdichte wird maximal, wenn die Summe im Exponenten (die als Summe von Quadraten ja nicht negativ sein kann) minimal wird. Man wird also die Parameter l1, l2, ..., lM so bestimmen, daß die Größe


s = N

i=1 
1
si2
[yi-y(l1,l2,,lM,xi)]2
ein Minimum wird. Dies ist die Methode der kleinsten Quadrate.


Gütekriterium

Die aus dem Maximum Likelihood-Prinzip abgeleitete Methode der kleinsten Quadrate hat den Vorteil, daß man unter den gemachten Voraussetzungen sofort ein Gütekriterium für den Fit formulieren kann. Wären nämlich die exakten Parameter l1, l2, ..., lM des theoretischen Modells bekannt, dann wären die Variablen


hi = [yi-y(l1,l2,,lM,xi)]/si
N(0,1)-verteilt (normalverteilt mit Mittelwert 0 und Varianz 1) und s, als Summe der Quadrate von N unabhängigen N(0,1)-verteilten Zufallsvariablen, wäre c2-verteilt mit N Freiheitsgraden. Die Wahrscheinlichkeit für das Auftreten des beobachteten oder eines noch größeren Wertes der Abweichungsquadratsumme s ist dann


P(t > s) = 1 - P(t s) = 1 - 1
G(N/2)

s/2

0 
dt tN/2-1e-t
wobei G(z) die Gamma-Funktion ist, d.h. G(1/2)={p}, G(1)=1 und G(z+1)=zG(z). Damit das theoretische Modell mit den Messungen verträglich ist, sollte diese Wahrscheinlichkeit nicht zu klein sein, z.B. nicht kleiner als 50 % (sonst hätten die Meßergebnisse nie zustande kommen können).

In Wirklichkeit kennt man jedoch die la nicht a priori, sondern bestimmt sie erst durch Anpassung an die Daten. In diesem Fall sind zwar die yi (bzw. hi) nicht mehr alle voneinander unabhängig, man kann aber für den Fall normalverteilter Variablen zeigen, daß s dann noch immer c2-verteilt ist, allerdings mit


n = N - M
(Anzahl der Beobachtungen minus Anzahl der Parameter) Freiheitsgraden. Da für große n der Median der c2-Verteilung von der Größenordnung n ist, geht man in der Praxis oft nach der Faustregel vor, daß man ein ``reduziertes c2'', s/(N-M), berechnet. Ist s/(N-M) wesentlich größer als 1, so ist das theoretische Modell nicht mit den Daten verträglich (oder man hat die Meßfehler unterschätzt); ist s/(N-M) wesentlich kleiner als 1, dann hat man zu viele Parameter angepaßt (oder die Meßfehler überschätzt).


Normalgleichungen

Zur tatsächlichen Bestimmung der Parameter müssen die simultanen Extremumsbedingungen


s
la
= 0,    a = 1,,M
nach l1, l2, ..., lM aufgelöst werden. Das ist im allgemeinen Fall ein System von M gekoppelten transzendenten Gleichungen, dessen Lösung sich beliebig schwierig gestalten kann.

Wesentlich einfacher ist die Situation, wenn y(x) ein lineares Modell ist,


y = M

a = 1 
la ja(x)
wobei die ja(x) beliebige Funktionen (häufig Polynome) sein können. D.h. y ist eine Linearkombination der ja, in der die Parameter la als Koeffizienten auftreten. Differenzieren von


s = N

i=1 
1
si2
 
yi- M

b = 1 
lbjb(xi)
2
 
nach la liefert in diesem Fall


s
la
= - 2 N

i=1 
1
si2
 
yi- M

b = 1 
lbjb(xi)
 ja(xi) = 0
oder, umgeformt,


M

b = 1 


N

i=1 
1
si2
ja(xi)jb(xi)

 lb = N

i=1 
1
si2
 ja(xi)yi
Mit den Abkürzungen


Aab
=
N

i=1 
1
si2
ja(xi)jb(xi)
ba
=
N

i=1 
1
si2
 ja(xiyi
ist das ein lineares Gleichungssystem der Form


Al = b
Dieses System der Normalgleichungen kann mit dem Gaußschen Eliminationsverfahren oder einem anderen gängigen Verfahren zur Lösung linearer Gleichungssysteme nach den la aufgelöst werden:


l = A-1b

Falls die ja(x) Polynome sind, kann schon für M ~ 5-10 das Lösungsverfahren numerisch problematisch werden. In diesem Fall empfiehlt es sich, zunächst orthogonale Polynome ya(x) mit


N

i=1 
1
si2
ya(xi)yb(xi) = dab
zu konstruieren, das Ausgleichsproblem in diesen zu lösen (A ist dann eine Diagonalmatrix) und die Lösung auf die ja(x) umzurechnen.


Fehlerabschätzung

Für ein lineares Modell und unter der Voraussetzung, daß das Modell mit den Messungen verträglich ist, kann man auch noch Aussagen über die Unsicherheit der mit Hilfe eines Fits bestimmten Parameter machen.

Dazu benützt man, daß die Matrix A in den Normalgleichungen zwar von den Stützstellen xi und den (als von den Daten unabhängig angenommenen) Fehlern si abhängt, aber nicht von den eigentlichen Meßwerten yi. Für feste xi und si sind daher die la lineare Funktionen der yi,


la = M

b = 1 
(A-1)ab N

i=1 
1
si2
 jb(xiyi
und als solche ebenfalls Zufallsvariable.

Stellt man sich nun vor, daß man den Prozeß des Messens und Anpassens viele Male wiederholt, so bekommt man für jeden Satz von Messungen {y1,y2,,yN} einen zugehörigen Satz von Parametern {l1,l2,,lM}. Der Mittelwert des Parameters la über alle diese Meßserien ist dann


la = M

b = 1 
(A-1)ab N

i=1 
1
si2
 jb(xiyi
Ein Maß für die Abweichung von la in einem einzelnen Fit vom Mittelwert ist die Varianz


sla2
=
(la-la)2
=
M

b = 1 
(A-1)ab M

g = 1 
(A-1)ag N

i=1 
1
si2
 jb(xi) N

j=1 
1
sj2
 jg(xj)(yi-yi)(yj-yj)
Wegen der statistischen Unabhängigkeit von yi und yj für i j ist


(yi-yi)(yj-yj) = dij si2
und daher


sla2
=
M

g = 1 
(A-1)ag M

b = 1 
(A-1)ab N

i=1 
1
si2
 jb(xi)jg(xi)
=
M

g = 1 
(A-1)ag M

b = 1 
(A-1)ab Abg
=
M

g = 1 
(A-1)agdag
Die Unsicherheit der Parameter ist also durch die Diagonalelemente der Matrix A-1 gegeben:


sla2 = (A-1)aa

Eine analoge Rechnung ergibt


cov(la,lb) = (A-1)ab
Dabei ist die Kovarianz


cov(la,lb) = (la-la)(lb-lb)
ein Maß für die Korrelation zwischen den Parametern.




File translated from TEX by TTH, version 2.86.
On 9 Jun 2005, 11:18.