Lineare Ausgleichsprobleme

Datenreduktion und Parameteranpassung

  

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 x₁, x₂, ..., x_N einer unabhängigen Variablen x zugehörige Werte y₁, y₂, ..., y_N 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 x_i 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

  

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

{xi,yi,si}i=1,¼, N

gegeben. Dabei wird angenommen, daß die x_i exakt bekannt, die y_i aber mit einem statistischen Meßfehler s_i behaftet sind. (Die Größenordnung von s_i kann man z.B. abschätzen, indem man für jedes x_i mehrere Messungen von y durchführt und den Standardfehler des Mittelwerts bestimmt.) Weiters wird angenommen, daß die Meßwerte y_i statistisch unabhängig und normalverteilt sind, d.h. die Wahrscheinlichkeit (eigentlich Wahrscheinlichkeitsdichte), als i-ten Meßpunkt den Wert y_i 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(x<sub>i</sub>) wäre der ``wahre'' Funktionswert an der Stelle x_i, bei einer tatsächlichen Messung wird aber ein davon abweichender Wert y_i beobachtet. Zusätzlich wird noch vorausgesetzt, daß die s_i unabhängig von den im konkreten Fall gemessenen Werten y_i sein sollen.

Das theoretische Modell hängt im allgemeinen von M Parametern l₁, l₂, ..., l_M ab

y = y(l1,¼,lM,x)

Nach dem Prinzip der Maximum Likelihood-Schätzung werden die l_a 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 y_i statistisch unabhängig sein sollen, ist die Wahrscheinlichkeitsdichte für das Auftreten der Meßwerte y₁, y₂, ..., y_N 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 l₁, l₂, ..., l_M 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 l₁, l₂, ..., l_M 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 c²-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 l_a nicht a priori, sondern bestimmt sie erst durch Anpassung an die Daten. In diesem Fall sind zwar die y_i (bzw. h_i) nicht mehr alle voneinander unabhängig, man kann aber für den Fall normalverteilter Variablen zeigen, daß s dann noch immer c²-verteilt ist, allerdings mit

n = N - M

(Anzahl der Beobachtungen minus Anzahl der Parameter) Freiheitsgraden. Da für große n der Median der c²-Verteilung von der Größenordnung n ist, geht man in der Praxis oft nach der Faustregel vor, daß man ein ``reduziertes c²'', 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 l₁, l₂, ..., l_M 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 j_a(x) beliebige Funktionen (häufig Polynome) sein können. D.h. y ist eine Linearkombination der j_a, in der die Parameter l_a als Koeffizienten auftreten. Differenzieren von

s = N å i=1  1 si2   é ë yi- M å b = 1  lbjb(xi) ù û 2

nach l_a 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(xi) yi

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 l_a aufgelöst werden:

l = A-1b

Falls die j_a(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 y_a(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 j_a(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 x_i und den (als von den Daten unabhängig angenommenen) Fehlern s_i abhängt, aber nicht von den eigentlichen Meßwerten y_i. Für feste x_i und s_i sind daher die l_a lineare Funktionen der y_i,

la = M å b = 1  (A-1)ab N å i=1  1 si2  jb(xi) yi

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 {y₁,y₂,¼,y_N} einen zugehörigen Satz von Parametern {l₁,l₂,¼,l_M}. Der Mittelwert des Parameters l_a über alle diese Meßserien ist dann

álañ = M å b = 1  (A-1)ab N å i=1  1 si2  jb(xi) áyiñ

Ein Maß für die Abweichung von l_a 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 y_i und y_j 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.