# Die Gauß'sche Methode der kleinsten Quadrate

Die folgenden Beispiele verwenden die Gauß'sche Methode der kleinsten Quadrate, um eine Linearkombination gegebener Funktionen zu bestimmen, die eine bestimmte Funktion möglichst gut approximiert.

## Das Problem

Angenommen, wir beobachten ein Objekt, das sich auf einer Geraden durch die Ebene bewegt. Drei aufeinanderfolgende Messungen liefern die Bahnpunkte (3,3), (6,3) und (9,6).
Wie die Abbildung zeigt, gibt es keine Gerade durch diese drei Messpunkte. Man könnte nun einfach einen Messwert ignorieren und bekäme je nach Wahl eine der drei roten Geraden. Bei einem fehlerbehafteten Messgerät werden aber alle Messungen ähnliche Abweichungen haben, so dass eine vermittelnde Gerade in der Regel zu einem besseren Ergebnis führt.
In der Abbildung ist die maximale Abweichung der blauen Geraden von den Messpunkten kleiner als bei jeder der drei roten Geraden.

Konkret suchen wir eine Gerade

\def\yellow{\bbox[#ffffa0,2pt]} \def\green{\bbox[#c0ffc0,2pt]} \green{f(x)} = a\yellow x + b

mit den unbekannten Koeffizienten a0 und a1.
Durch Einsetzen der drei Messwerte bekommen wir:

\def\yellow{\bbox[#ffffa0,2pt]} \def\green{\bbox[#c0ffc0,2pt]} \begin{eqnarray} \yellow 3 a + b & = \green 3 \cr \yellow 6 a + b & = \green 3 \cr \yellow 9 a + b & = \green 6 \end{eqnarray}

In Matrixform bekommen wir Ax = b mit

A = \pmatrix{3 & 1 \cr 6 & 1 \cr 9 & 1 } \quad \textbf x = \pmatrix{a \cr b} \quad \textbf b = \pmatrix{3 \cr 3 \cr 6}

Dieses Gleichungssystem ist überbestimmt und nicht lösbar.

## Die Lösung

In der Vorlesung Lineare Algebra für Informatiker wird der folgende Satz gezeigt:

Satz

Das Normalsystem ATAx = AT b eines linearen Gleichungssystems Ax = b ist konsistent.
Seine Lösungen sind die Näherungslösungen von Ax = b mit projW(b) = Ax
wobei W der Spaltenraum von A ist.

Wir wenden nun den Satz auf unser Beispiel an.
Für AT schreibt man in mathGUIde `A.transp()` oder kürzer `~A`

    A = Matrix([[3,1], [6,1], [9,1]])
    b = ~Vector([3, 3, 6])
    (~A*A).solve(~A*b)

Damit erhalten wir die Gerade f(x) = 0.5 x + 1
Wir plotten diese Funktion und zeigen dazu die Messpunkte an:

    plot(["1+0.5*x"], 3.0, 9.0, mark=[(3,3),(6,3),(9,6)])

## Mehr Komfort: Die Funktion `fit`

Um uns den Matrixansatz zu ersparen, bietet mathGUIde die Funktion `fit` an, die aus den Messwerten und dem Funktionstyp direkt die fertige Funktion berechnet.

Die Funktion `fit` erwartet zwei Parameter

1. Eine Liste mit den Datenpunkten, jeweils (x,y)
2. Eine Liste mit Elementarfunktionen (als Zeichenkette geschrieben!), aus denen die gesuchte Funktion als Linearkombination zusammengesetzt wird (Variable x)

Für unser Beispiel übergeben wir:

    fit([(3,3),(6,3),(9,6)], ["1", "x"])

## Weitere Beispiele

### Beispiel 1

Gesucht ist eine Gerade der Form f(x) = ax + b, die die drei Punkte (3,3), (6,4) und (9,6) möglichst gut approximiert (Regressionsgerade).

    points = [(3,3), (6,4), (9,6)]
    # Gerade ist Linearkombination der beiden folgenden Funktionen:
    functions = ["x", "1"]
    print(fit(points, functions))

mathGUIde hat (hier in etwas vereinfachter Form) die Funktion f(x) = x/2 + 4/3 geliefert.
Zur Kontrolle der Approximation schauen wir uns einen Funktionsplot an:

    points = [(3,3), (6,4), (9,6)]
    functions = ["x", "1"]
    plot(fit(points, functions), 3.0, 9.0, mark=points)

### Beispiel 2

Eine Parabel soll an vier Punkte angenähert werden:

    points = [(-1,2), (0,2), (1,1), (2,0)]
    functions = ["1", "x", "x^2"]
    print(fit(points, functions))

Kontrolle des Ergebnisses:

    points = [(-1,2), (0,2), (1,1), (2,0)]
    functions = ["1", "x", "x^2"]
    plot(fit(points, functions), -1.0, 2.0, mark=points)

### Beispiel 3

Transzendente Funktion: `f(x) = a + b*x*log(x) + c*exp(x)`
Gesucht sind die Koeffizienten a, b, c

    functions = ["1", "x*log(x)", "exp(x)"]
    # Die Funktion sollte die folgenden 4 Punkte approximieren:
    points = [(1,1), (2,1), (3,3), (4,8)]
    print(fit(points, functions))

Kontrolle des Ergebnisses:

    functions = ["1", "x*log(x)", "exp(x)"]
    points = [(1,1), (2,1), (3,3), (4,8)]
    plot(fit(points, functions), 1.0, 4.0, mark=points)