• Reguläres,D-dimensionales Gitter aufgrund der exponentiellen Skalierung nur für bis zu dreidimensionale Probleme möglich.
• Dafür wurden spezielle Methoden wie Triangulierung (in 2D) entwickelt.
• Hier kurz behandelt: Interpolation von verstreuten Punkten im Mehrdimensionalen (keine exponentielle Skalierung).
10.6.1 Interpolation mit radialen Basisfunktionen
• Approximative Annahme: Jeder Stützpunkt xi, y(xi)
„beinflusst“ seine lokale Umgebung unabhängig von der Richtung.
• Der Funktionswert in der Umgebung vom Stüzpunkt hängt dann nur vom radialen Abstand kx−xik2 ab.
• Dementsprechend lautet der Interpolationsansatz
y(x) =
n−1
X
i=0
ωiφ(kx−xik2) (10.89)
bei gegebener 1D(!) Basisfunktionφund zu bestimmenden Wichtungenωi.
• Ermittlung der Wichtung durch Einsetzten der Interpolationsbedingung
y(xj)=!
n−1
X
i=0
ωiφ(kxj −xik2), j ∈[0, n−1]. (10.90) Dies führt zu einem linearen Gleichungssystem. Lösung dieses Systems skaliert also kubisch(!) mit der Zahl der Punkte.
• Diese Art kann also nur für wenige Punkte bis zun ∼104 noch vernünftig verwendet werden.
• Weitere Variante: Normierte Basisfunktion:
y(x) =
• Viele Basisfunktionenφsind möglich. Am Häufigsten verwendet:
– Multiquadratische Funktion:φ(r) = p
r2+r20mit Skalierungsparameterr0(sollte größer sein als der typische Abstand zwischen zwei Stützstellen).
– Genauso möglich: Inverse multiquadratische Funktionφ(r) = √ 1
r2+r02.
145
10.6.2 Shepard Interpolation/Inverse Distanzwichtung
• Nehme normierte, radiale Basisfunktionen (Gl.(10.91)) undφ(r)sodassφmitr >0monoton abnimmt und fürr→0gegen Unendlich geht.
• Dann ergibt die Interpolationsbedingung
y(xj)= lim!
x→xj
Pn−1
i=0 ωiφ(kx−xik2) Pn−1
i=0 φ(kx−xik2) = lim
x→xj
ωjφ(kx−xjk2)
φ(kx−xjk2) =ωj. (10.93)
• Die Wichtungenωj sind also die Funktionswerte:
y(x) =
(Pn−1
i=0y(xi)φ(kx−xik2) Pn−1
i=0 φ(kx−xik2) ,∀j ∈[0, n−1],x6=xj,
y(xj) ,x=xj. (10.94)
• Diese Methode erfordert kein Lösen eines Gleichungssystems und der Aufwand skaliert linear mit der Zahl der Stützstellen, allerdings weniger genau als normale Interpolation per radialer Basisfunktion.
• Oftmals verwendete Basisfunktion:φ(r) =r−p, mitpalsreelleZahl zwischen1und3.
• Verbesserung möglich, indem für die Wichtung,ωi, anstelle der Funktionswertey(xi)eine Taylorentwicklung (bis zur 2. Ordnung) verwendet wird.3Die Wichtung ist dann (wie die radiale Basisfunktion) abhängig vom Abstand zur Stützstelle. Diese Art der Shepard-Interpolation ist dann letztlich eine gewichtete Mittelung von Taylorreihen an den Stützstellen.
• Großer Vorteil: Zusätzliche Stützstellen können ohne zusätzlichen Aufwand hinzugenommen werden.
3Genauer: Taylorentwicklung in reziproken Längen (1/x), siehe J. Ischtwan, M. A. Collins,J. Chem. Phys.,100, 8080 (1994).
11 Ausgleichungsrechnung
• Gegeben:N Datenpunkte(xi, yi), i= 1, . . . , N, oftmals mitFehlernbehaftet.
• (möglichst sinnvolles)Modell:y(x) =y(x;a1, a2, . . . , aM)mitM anzupassenden Parametern.
• Gesucht: Parameterai, quantifizierbares Maß für Güte des Ausgleichs.
• Nomenklatur:yi ist ein Datenpunkt,y(xi)der entsprechende Wert der Modellfunktion.
• Normalerweise istN M.
• Generell gilt: Der Ausgleich kann bei Datenpunkten mit Fehlern nie perfekt sein.
• Daher Unterschied zur Interpolation: Die Modellfunktiony(x)geht nicht durch die Daten-punkte. Stattdessen versucht man, die quadrierte Differenz[y(xi)−yi]2unter Beachtung des Fehlers vony˜i zu minimieren: Methode der kleinsten Quadrate,least squares fitting.
• Interpolation ist aufgrund der Fehler der Datenpunkte und aufgrund der mehr oder weniger guten Modellfunktion auchnicht sinnvoll:
• Für spezielle Modellfunktionen gibt es optimierte Verfahren:
– Linear:y(x) =a1x+a2 – Polynomial:y(x) =PN
i=0aixi – Rational:y(x) =
PN j=0ajxj
/
PN i=0aixi
– Gauß-Funktion:y(x) =PN
i=1aiexp
−x−x
0i
σi
2 – ...
• Für hochdimensionale Probleme sehr viel aktuelle Forschung (Neuronale Netzwerke,Machine Learning, ...).
147
11.1 Hintergrund: Statistik
• Naheliegende Frage (A): „Wie wahrscheinlich ist es, dass ein Parametersatz{ai}Ni=1korrekt ist?“ ist strenggenommen sinnlos, denn:
– nicht nur die Parameter, auch das Modell kann falsch sein;
– es gibt nureinenwahren funktionalen Zusammenhang miteinemParametersatz; beide sind prinzipiell unbekannt.
• Besser ist die umgedrehte Frage (B): „Gegeben sei ein Modell und ein Parametersatz. Wie wahrscheinlich ist dann der vorliegende Datensatz?“
• Mögliche Antworten darauf erlauben Folgerungen:
– (sehr) wahrscheinlich: Parametersatz o.k., Modell (evtl.) auch;
– (sehr) unwahrscheinlich: wenn keine vernünftigen Parameter gefunden werden können, ist das Modell zu verbessern.
• DefinitionFit:Maximiertdie obige Wahrscheinlichkeit (B).
• Annahmen:Die Fehler an jedem Messpunkt seien – zufällig (nicht systematisch),
– unabhängig,
– normalverteilt (Gauß-verteilt)1, mit Standardabweichungσi.
• Dann ist die zu maximierende GesamtwahrscheinlichkeitP (Daten|Modell)für diesen Da-tensatz das Produkt aus den Einzelwahrscheinlichkeiten:
P (Daten|Modell) =
N
Y
i=1
( exp
"
−1 2
yi−y(xi;a1, . . . , aM) σi
2#
∆y )
. (11.1)
Dabei müssen wir für jeden Messpunkt eine (feste) Abweichungsbreite∆yzulassen, weil es beliebig unwahrscheinlich ist, dass alle Datenpunkte mathematisch exakt auf der Modellkurve liegen (P = 0für∆y= 0).
• Erinnerung:P (A|B)ist die bedingte Wahrscheinlichkeit, dass EreignisAeintritt, unter der Bedingung, dass das EreignisBschon eingetreten ist.
• Erinnerung: Die Standardabweichungσiist ein Maß für die Breite der Gaußkurve an jedem Messpunkti, nach folgenden Faustregeln1: 68% aller Datenpunkte sind im Intervall±σi, 95%
in±2σi, 99.7% in±3σi, usw.
1Nur korrekt im Limit sehr großer Zahlen; typische Meßpunkte sind eher Poisson-verteilt = höhere Wahrscheinlichkeit für größere „Ausreißer“.
11.1.1 χ
2-Statistik
• Maximierung vonP ist äquivalent zu Minimierung des negativen Logarithmus vonP: 1
• Weglassen der konstanten Terme liefert die neue Größe χ2 =
• Mitσi =const. ist Minimierung vonχ2äquivalent zur Minimierung der Fehlerquadratsumme (least squares fit).
• χ2 ist auch wieder eine Wahrscheinlichkeitsverteilung und beschreibt die Verteilung der Summe vonνquadrierten, unabhängigen Zufallsvariablen. Die Dichte der Verteilung sieht wie folgt aus2
• Hier istνdie Zahl der Freiheitsgrade (N −M).
• Für diese funktionale Form lässt sich analytisch die WahrscheinlichkeitQdafür berechnen, dassχ2größer als der aus dem Fit resultierendeχ2-Wert ist:
Q χ2
Dabei istQΓ(a, x)die unvollständigeΓ-Funktion:
QΓ(a, x) = 1−PΓ(a, x) = Γ(a, x)
wobei die vollständigeΓ-Funktion gegeben ist durch Γ(a) =
• Bei einem akzeptablen Fit giltχ2 ≈ν. Vorsicht:χ2 ν bedeutet meistensover-fitting: Die Fehler der Daten wurden entweder überschätzt, oder das Modell beinhaltet zu viele Variablen, sodass das Rauschen der Daten (noise) mitbeschrieben wird.
• WennQgroß ist, dann ist die Wahrscheinlichkeit, ein größeresχ2zu finden, groß (Definition vonQ). Daχ2 möglichst klein sein soll, soll entsprechendQmöglichst groß sein.
• Q-Werte liefern daher ein quantitatives Maß für die Fit-Güte. Faustregeln:
– Q <0.001: Modell falsch?
– Q& 0.001Fit evtl. akzeptabel, oderσi unterschätzt, oder Abweichungen nicht Gauß-verteilt (unsichere Folgerung, da das hier vorausgesetzt wurde)
– Q&0.1: Fit vermutlich o.k.
– Q-1: Fehler überschätzt, oder Betrug.
Abschätzung der Standardabweichung
• Die Standardabweichungσist oftmals nicht bekannt und wird oft nachträglich abgeschätzt:
– Setzeσi =σ(also gleich für alle Meßpunkte).
– Führe die Regression durch.
– Berechneσ2 =PN
i=1[yi−y(xi)]2/N.
• Dabei sollte man jedoch beachten, daß
– Diesesσnur dann sinnvoll ist, wenn der Fit in Ordnung ist;
– Dannχ2nicht mehr unabhängig von diesemσals Fit-Gütemaß verwendet werden kann.
Durchführung der Minimierung
• Minimierung vonχ2 erfolgt in der üblichen Weise:
Setze die 1. Ableitung von χ2 nach den Parametern ak, k = 1,2, . . . , M, gleich Null. Wir erhalten (Kettenregel):
0 =
N
X
i=1
yi−y(xi;a1, . . . , aM) σi2
∂y(xi;a1, . . . , aM)
∂ak
, fürk = 1,2, . . . , M. (11.7)
• Das sindM Gleichungen für dieM unbekanntenak-Werte.
• Je nach der Art des Modellsy(x) =y(x;a1, a2, . . . , aM)können diese Gleichungen sein – nicht-linear
– linear.