Statistische Methoden der Datenanalyse
Vorlesung PHY523, Fakultat Physik, Technische Universitat Dortmund
Wolfgang Rhode
11. Oktober 2010
Inhaltsverzeichnis
Vorbemerkung i
Motivation iii
I. Grundlagen 1
1. Numerische Grundlagen 31.1. Arithmetische Ausdrucke . . . . . . . . . . . . . . . . . . . . . . . . . . . 41.2. Zahlen, Operationen und elementare Funktionen am Computer . . . . . . 7
1.2.1. Ganze Zahlen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71.2.2. Gleitpunktzahlen . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81.2.3. Operationen und Funktionen . . . . . . . . . . . . . . . . . . . . . 11
1.3. Stabilitat . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161.3.1. Entwicklung von Kriterien fur die numerische Stabilitat . . . . . . 16
1.4. Fehlerfortpflanzung und Kondition . . . . . . . . . . . . . . . . . . . . . . 221.4.1. Vergleich von Kondition und Stabilitat . . . . . . . . . . . . . . . . 24
2. Wahrscheinlichkeit und Generatoren 272.1. Vorbemerkungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.1.1. Definitionen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 272.1.2. Kombination von Wahrscheinlichkeiten . . . . . . . . . . . . . . . . 28
2.2. Zufallsvariable und deren Verteilung . . . . . . . . . . . . . . . . . . . . . 292.3. Allgemeine Eigenschaften einer Zufallsvariablen: Erwartungswert, Streu-
ung, Momente, etc. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 312.3.1. Regeln uber Mittelwerte und Varianzen . . . . . . . . . . . . . . . 33
2.4. Gleichverteilung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 362.5. Erzeugung von gleich- und beliebig verteilten Zufallszahlen auf dem Com-
puter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 372.5.1. Lineare kongruente Generatoren (LCG) . . . . . . . . . . . . . . . 372.5.2. Multiplikativ linear kongruente Generatoren MLCG . . . . . . . . 382.5.3. Spektraltest . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 392.5.4. Erzeugung beliebig verteilter Zufallszahlen (Teil 1) . . . . . . . . . 41
3. Spezielle Wahrscheinlichkeitsdichten 453.1. Gleichverteilung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 453.2. Die Binomialverteilung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3
Inhaltsverzeichnis
3.3. Die Normal- oder Gauß-Verteilung . . . . . . . . . . . . . . . . . . . . . . 47
3.4. Die Poisson-Verteilung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.5. Die Gamma-Verteilung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.6. Die χ2-Verteilung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.7. Die Cauchy-Verteilung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.8. Die t-Verteilung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.9. Die F-Verteilung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4. Beliebige verteilte Zufallszahlen (Teil 2) 554.1. Transformation der Gleichverteilung . . . . . . . . . . . . . . . . . . . . . 55
4.2. Das Neumannsche Ruckweisungsverfahren . . . . . . . . . . . . . . . . . . 55
4.3. Erzeugung normalverteilter Zufallszahlen . . . . . . . . . . . . . . . . . . 57
4.4. Erzeugung Poisson-verteilter Zufallszahlen . . . . . . . . . . . . . . . . . . 59
4.5. Erzeugung χ2-verteilter Zufallszahlen . . . . . . . . . . . . . . . . . . . . . 60
5. Mehrdimensionale Verteilungen 615.1. Problemstellung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
5.2. Erwartungswert, Varianz, Kovarianz und Korrelation bei zwei Variablen . 62
5.3. Mehrere Veranderliche . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
5.4. Die mehrdimensionale Gauß-Verteilung . . . . . . . . . . . . . . . . . . . . 68
6. Einfache statistische Methoden 716.1. Trennung von Datensatzen: Diskriminanzanalyse . . . . . . . . . . . . . . 71
6.2. Theoreme und Satze . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
6.2.1. Tschebyscheff-Ungleichung . . . . . . . . . . . . . . . . . . . . . . 74
6.2.2. Gesetz der großen Zahl . . . . . . . . . . . . . . . . . . . . . . . . 75
6.2.3. Der Zentrale Grenzwertsatz . . . . . . . . . . . . . . . . . . . . . . 76
6.3. Methode der kleinsten Quadrate . . . . . . . . . . . . . . . . . . . . . . . 77
6.3.1. Vorbemerkungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
6.3.2. Kleinste Quadrate in linearen Modellen . . . . . . . . . . . . . . . 78
6.3.3. Gaußverteilte Meßfehler . . . . . . . . . . . . . . . . . . . . . . . . 82
6.3.4. Nichtlineare kleinste Quadrate . . . . . . . . . . . . . . . . . . . . 86
6.4. Nachtrag und Exkurs: Fehlerfortpflanzung . . . . . . . . . . . . . . . . . . 87
6.4.1. Transformation einer Variablen . . . . . . . . . . . . . . . . . . . . 87
6.4.2. Transformation mehrerer Variablen . . . . . . . . . . . . . . . . . . 89
6.5. Numerische Optimierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
6.5.1. Vorbemerkungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
6.5.2. Eindimensionale Minimierung . . . . . . . . . . . . . . . . . . . . . 95
6.5.3. Mehrdimensionale Minimierung . . . . . . . . . . . . . . . . . . . . 97
7. Spezielle Verfahren zur Datenanalyse 997.1. Die Maximum-Likelihood-Methode . . . . . . . . . . . . . . . . . . . . . . 99
7.1.1. Problemstellung . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
7.1.2. Aufgabe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
4
Inhaltsverzeichnis
7.1.3. Ansatz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 997.1.4. Maximum-Likelihood-Prinzip . . . . . . . . . . . . . . . . . . . . . 997.1.5. Log-Likelihood-Funktion . . . . . . . . . . . . . . . . . . . . . . . . 1007.1.6. Beispiele . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
7.2. Fehlerbestimmung bei der Maximum-Likelihood-Methode . . . . . . . . . 1017.2.1. Ein Parameter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1017.2.2. Mehrere Parameter . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
7.3. Die Maximum-Likelihood-Methode, Eigenschaften . . . . . . . . . . . . . 1037.3.1. Konsistenz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1037.3.2. Erwartungstreue? . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1037.3.3. Gaußahnlichkeit . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1037.3.4. Varianz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
7.4. Bayesische Statistik . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1037.5. Entfaltung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
7.5.1. Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1067.5.2. Akzeptanzkorrektur . . . . . . . . . . . . . . . . . . . . . . . . . . 1087.5.3. Diskretisierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1087.5.4. Anwendung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1107.5.5. Entfaltung (ohne Regularisierung) . . . . . . . . . . . . . . . . . . 1117.5.6. Problemanalyse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1127.5.7. Regularisierung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
5
Vorbemerkung
Das vorliegende Skript basiert auf der Beschreibung des Moduls PHY523 im Bachelor-bzw. Masterstudiengang Physik an der Technischen Universitat Dortmund. Die Veran-staltung richtet sich an Physikstudierende im funften (B.Sc.) bzw. ersten und dritten(M.Sc) Semester. Ziel dieser dreistundigen einsemestrigen Veranstaltung ist eine Ver-mittlung von Kompetenzen, die fur die Erstellung der Bachelor- bzw. Masterarbeit unddas spatere Berufsleben relevant sind. Vermittelt werden soll ein geeigneter Umgangmit statistischen Methoden zur Analyse von moderaten bis sehr großen Datenmengen.Ein Teil der Ubungsaufgaben, die in den zugehorigen zweistundigen Ubungen bespro-chen werden, konnen und sollen unter Einbeziehung des von in der Experimentalphysikbenutzten Datenverwaltungs- und Analysesystems root auch am Computer gelost wer-den. In diesem Zusammenhang werden Grundkenntnisse der Programmiersprache C++benotigt bzw. erlernt.
Im Rahmen des Bachelor- bzw. Masterstudiengangs werden mit einem generischenArbeitsaufwand von 240 Stunden, von denen 81 Stunden auf Ubungen sowie die Ab-schlussprufung entfallen, 8 Credits erworben.
i
Motivation
In der modernen Experimentalphysik werden Messdaten in der Regel auf elektronischemWeg erhoben und verarbeitet. Vor ihrer Durchfuhrung werden kostenintensive Experi-mente und Messungen in der Regel mit statistischen Verfahren (Monte Carlo) geplant.Dieselben Verfahren werden auch zur Interpretation und Analyse eingesetzt. In dieserVorlesung wird Basiswissen vermittelt, das zur Planung eines Experimentes, der Analy-se der erfassten Daten sowie schließlich zur Extraktion von physikalischen Parameternbenotigt wird. Der Aufbau der Vorlesung folgt einer gedanklichen Planung und Ana-lyse eines Experiments. Daher werden zunachst die numerischen Randbedingungen dernumerischen Mathematik, dann der Wahrscheinlichkeitstheorie bis hin zu Monte Car-lo Simulationen von Experimenten und schließlich der Datenanalyse mit einfachen undkomplexen Verfahren behandelt.
Die Vorlesung umfasst folgende Lehrinhalte:
• Numerische Methoden der Datenverarbeitung,
• Datenbehandlung und Programmierung,
• Algorithmen und Datenstrukturen,
• Methoden der linearen Algebra,
• Wahrscheinlichkeitsrechnung,
• Ein- und mehrdimensionale Verteilungen,
• Zufallszahlen und Modellieren von Experimenten mit Monte Carlo Methoden,
• Parameterschatzung,
• Optimierungsprobleme,
• Die Methode der kleinsten Quadrate,
• Die Maximum Likelihood-Methode,
• Konfidenzintervalle und Hypothesentests,
• Parametrisierung von Daten,
• Bayes’sche Verfahren,
iii
Motivation
• Fehlerfortpflanzung
• Entfaltung
• Neuronale Netze
• Zufallige Walder
• Boosted Decision Trees
• Supported Vector Machines (SVN)
Zusatzlich zu dem vorliegenden Skript seien folgende Lehrbucher empfohlen:
• R.J. Barlow, Statistics, Wiley;
• V. Blobel, E. Lohrmann, Numerische und Statistische Methoden der Datenanalyse,Teubner;
• S. Brand, Datenanalyse, Spektrum Verlag;
• T. Butz, Fouriertransformation fur Fußganger, Teubner;
• G.D. Cowan, Statistical Data Analysis, Oxford University Press;
• W.T. Eadie et al., Statistical Methods in Experimental Physics, North-Holland;
• H.L. Harney, Bayesian Inference, Springer;
• F. James, Telling the truth with Statistics, CERN Academic Training Programm;
• F. James, Statistical Methods in Experimental Physics, World Scientific;
• D.E. Knuth, The Art of Computer Programming, Addison Wesley;
• L. Lyons, Statistics for nucelar and particle physicists, Cambridge University Press;
• W.T. Press et al., Numerical Recipes, Cambridge University Press;
• D.S. Sivia, Data Analysis - A Bayesian Tutorial, Oxford University Press
• to be completed
Von IP-Adressen der Technischen Universitat Dortmund sind manche Lehrbucher desSpringer-Verlages frei zuganglich.
Das deutschsprachige Paket Naturwissenschaften besteht aus der folgenden Samm-lung:
Mathematics and Statistics:
http://www.springerlink.com/mathematics-and-statistics/?Content+Type=Books&Copyright=2006&sortorder=asc&Language=German
Earth and Environmental Science:
http://www.springerlink.com/earth-and-environmental-science/?Content+Type=Books&Copyright=2006&sortorder=asc&Language=German
Chemistry and Materials Science:
http://www.springerlink.com/chemistry-and-materials-science/?Content+Type=Books&Copyright=2006&sortorder=asc&Language=German
Biomedical and Life Science:
http://www.springerlink.com/biomedical-and-life-sciences/?Content+Type=Books&Copyright=2006&sortorder=asc&Language=German
Physics and Astronomy:
http://www.springerlink.com/physics-and-astronomy/?Content+Type=Books&Copyright=2006&sortorder=asc&Language=German
iv
Dazu gehoren:
• H.R. Schwarz, N. Kockler, Numerische Mathematik, Teubner-Verlag 2006
http://www.springerlink.com/content/h58j45
• J.-P. Kreiß, G. Neuhaus, Einfuhrung in die Zeitreihenanalyse, Springer 2006
http://www.springerlink.com/content/p47746
• und andere
Die Ubungsaufgaben sollten -sofern eine numerische Behandlung erforderlich ist- mitHilfe des objektorientierten Datenverarbeitungspaketes root gelost werden. Hierzu sindGrundkenntnisse in C++ erforderlich. In der Teilchen- und Astroteilchenphysik wirdroot als Standard eingesetzt. Fur die gangigen Betriebssysteme kann root bei Beachtungder dort gegebenen copy right Informationen von der Webseite
http://root.cern.ch
heruntergeladen werden. Das root-Handbuch ist unter
http://root.cern.ch/root/doc/RootDoc.html
zu finden. Entsprechende Tutorials stehen unter
http://root.cern.ch/root/Tutorials.html
bereit.
v
Teil I.
Grundlagen der DatenanalyseBruderchen, Du schwatzest zu sub-til: Du grubelst und grubelst undhast am Ende nichts als Unruhe undUngewisheit zum Lohne. Ich glau-be frisch weg, ohne mich links oderrechts umzusehen, daß alles gut undweise angeordnet ist. Ich komme ambesten dabey zu rechte. Ist auchwirklich alles Nothwendigkeit undZufall; muß ich mich von diesen bey-den Machten herumstoßen laßen -wohlan! ich wills gar nicht wissen,daß sie mich blind herumstoßen. DerKopf wird so dadurch wirblicht ge-nug, soll ich mir ihn noch durchGrubeleyen wirblicht machen?Johann Karl Wezel, Belphegor, 1776
1
1. Numerische Grundlagen
Bei der computergestutzten Analyse von Daten geht es letztlich immer um eine moglichstgenaue Extraktion von physikalisch zu interpretierenden funktionalen Zusammenhangenaus einer erfassten Datenmenge. Die Genauigkeit der Analyse ist dabei durch drei Fak-toren eingeschrankt:
• Durch den statistischen Fehler. Wie dieser Fehler bestimmt und unter den gegebe-nen Umstanden minimiert wird, wird im Hauptteil dieser Vorlesung diskutiert.
• Durch den systematischen Messfehler. Dieser Fehler kann mit rein statistischenMitteln nur bedingt abgeschatzt und nicht reduziert werden. In dieser Vorlesungwerden daher nur Darstellungsfragen behandelt.
• Durch den numerischen Fehler. Dieser Fehler entsteht durch die Bearbeitung derDaten auf dem Computer und kann durch geeignete Maßnahmen bei der Pro-grammierung der Datenanalyse minimiert werden. Es ist daher notwendig, vor derDiskussion und Programmierung numerisch aufwandiger Prozeduren, die Grund-regeln der numerischen Datenverarbeitung zu besprechen.
3
1. Numerische Grundlagen
1.1. Arithmetische Ausdrucke
Zusammenhange in der Physik und Statistik werden durch mathematische Formeln be-schrieben.
Beispiele:
|~v| =√x2 + y2 Bahnradius (1.1)
E = E0(1 + ε)n Energiegewinn bei stochastischer Beschleunigung (1.2)
σ =
√∑x2i − (
∑xi)
2 /n
nStandardabweichung (1.3)
x1,2 =−b±
√b2 − 4ac
2aLosung einer quadratischen Gleichung (1.4)
h = arcsin(sinψ sin δ + cosψ cos δ cos t) Hohe eines Sterns (1.5)
θ =1
2π
x∫−x
e−t2
2 dt Flache unter der Normalverteilungskurve (1.6)
Formeln, in denen nur die Grundrechenarten und elementare Funktionen vorkommen,werden als arithmetische Ausdrucke bezeichnet.
Definition der arithmetischen Ausdrucke:
Es gelten folgende Bezeichnungen:
Variable: x1, x2, ..., xn ⊂ R .
Zweistellige Operationen: O = {+,−, ∗, /, ∗∗}.
Elementare Funktionen: F = {sin, cos, exp, ln, sqrt, abs, ...}.
Damit kann ein arithmetischer Ausdruck wie folgt definiert werden:
Die Menge A=A(x1, x2, ..., xn) der arithmetische Ausdrucke in x1, x2, ..., xn ist defi-niert durch
i) R ⊆ A
ii) xl ∈ A, fur l = 1, 2, ..., n
iii) g ∈ A ⇒ (−g) ∈ A
4
1.1. Arithmetische Ausdrucke
iv) g, h ∈ A, · ∈ O ⇒ (g · h) ∈ A
v) g ∈ A, φ ∈ F ⇒ φ(g) ∈ A
vi) A(x1, x2, ..., xn) ist minimal unter den Mengen A, die (i)− (v) erfullen.
Um spater Regeln fur eine geeignete Programmierung formulieren zu konnen, sollenzunachst Beispiele betrachtet werden. 1
Beispiel: Losung einer quadratischer Gleichung
y = [−b+ sqrt(b ∗ b− 4 ∗ a ∗ c)]/(2 ∗ a) ∈ A(a, b, c).
Zur Auswertung werden die Variablen durch Zahlen ersetzt.
Beispiel: Das Horner–Schema
Ein Polynom n-ten Grades ist eine Funktion der Gestalt
f(x) = a0xn + a1x
n−1 + ...+ an−1x+ an
mit vorgegebenen Koeffizienten a0, a1, ..., an. Die Auswertung moge an der Stelle xerfolgen.
• naive Vorgehensweise:
– Bildung aller Potenzen xk
– Multiplikation mit den Koeffizienten ai
– Addition
• Horner–Schema: Wir schreiben statt dessen:
fi := a0xi + a1x
i−1 + ...+ ai−1x+ ai fur i = 1, 2, ..., n.
Dann giltfn = f(x).
Der Vorteil dieser Darstellung besteht darin, dass die fi rekursiv definierbar sind.Es ist also im
Horner− Schema :
f0 = a0 ;fi = fi−1 · x+ ai, i = 1, 2, ..., n;f(x) = fn.
Da ein Polynom zu den arithmetischen Ausdrucken gehort, kann auch die Ableitungrekursiv gebildet werden.
1Es ist empfohlen, diese Beispiele in C++ bzw. root zu programmieren und in geeigneten Grafiken zuvisualisieren.
5
1. Numerische Grundlagen
Horner− Schema der ersten Ableitung
f ′0 = 0 ;f ′i = f ′i−1 · x+ fi−1, i = 1, 2, ..., n;f ′(x) = f ′n.
Beispiel:
Fur das Polynom f(x) = 4x2 + 2x+ 3 ist
f0 = 4 f ′0 = 0f1 = 4 · x+ 2 f ′1 = (0 · x) + 4f2 = (4 · x+ 2) · x+ 3 f ′2 = 4 · x+ (4 · x+ 2)
Das vollstandige Horner-Schema liefert den Funktionswert und samtliche Ableitungen:
a0 a1 a2 ... an1 an+0 +xf0 +xf1 ... +xfn−2 fn−1
f0 f1 f2 ... fn1 fn
+0 +xf ′1 +xf ′2 ... +xf ′n−1
f ′1 f ′2 f ′3 ... f ′n
Beispiel:
Betrachte f(x) = (1− x)6 und berechne die Losung dann:
• (1− x)6 einfach genau.
• 1− 6x+ 15x2 − 20x3 + 15x4 − 6x5 + x6 naiv.
• mit dem Horner-Schema.
• doppelt genau.
In Abbildung 1.1 ist die relative Abweichung zwischen der naiv berechneten Funktionund der in diesem Fall geeigneten zusammengefassten Berechnung dargestellt. Offenbarzerstort das Rechnen mit Maschinenzahlen hier die Monotonieeigenschaft von f . Ins-besondere bei naiver Berechnung treten viele Vorzeichenwechsel in der Umgebung vonx = 1 auf.
6
1.2. Zahlen, Operationen und elementare Funktionen am Computer
0.99 0.992 0.994 0.996 0.998 1 1.002 1.004 1.006 1.008 1.01
-410
-310
-210
-1101
10
210
310
410
510
610
710
810
910
abs((pow((1-x),6)-(1-6*x+15*pow(x,2)-20*pow(x,3)+15*pow(x,4)-6*pow(x,5)+pow(x,6)))/pow((1-x),6))
Abbildung 1.1.: Relative numerische Abweichung zwischen der naiven Berechnung undder zusammengefassten Berechnung des Polynoms. In der Umgebungder Nullstelle wird der Fehler beliebig groß.
1.2. Zahlen, Operationen und elementare Funktionen amComputer
Im Folgenden werden generische Eigenschaften von Maschinenzahlen, soweit fur die Da-tenanalyse relevant, diskutiert. Die genaue Speicherform hangt von der Wortbreite derverwendeten CPU, dem Betriebssystem und der Programmiersprache ab.2 Als Konse-quenz hangen auch die numerischen Eigenschaften einer Analyse oder Rechnung vondem Betriebssystem, der Programmiersprache und von Eigenschaften der Codierung ab.
1.2.1. Ganze Zahlen
Ganze Zahlen (Integer, ± xxxxx) werden als Kombination aus Vorzeichenbit und Ziffern-folge gespeichert. Die Codierung der Zahl erfolgt z.B. im 4-Bit-Zweierkomplement:
2Fur Details sei daher auf die entsprechenden Handbucher verwiesen.
7
1. Numerische Grundlagen
Bits Zahlenbereich Bezeichnung des Zahlentyps
8 -128 ... 127 Byte16 -32768 ... 32767 short integer32 −231...231 − 1 integer64 −263...263 − 1 long integer8 0 ... 255 unsigned byte
16 0 ... 65535 word
Tabelle 1.1.: Generische Deklarationen von ganzen Zahlen.
0 0000 -8 10001 0001 -7 10012 0010 -6 10103 0011 -5 10114 0100 -4 11005 0101 -3 11016 0110 -2 11107 0111 -1 1111
Der darstellbare Zahlenbereich hangt von der Zahl der zur Codierung verwendetenBits sowie der Entscheidung ab, ob ein Vorzeichenbit benotigt wird.
Es sei darauf hingewiesen, dass mit Analog-Digital-Convertern (ADC) (siehe z.B. 1.2)elektronisch registrierte Messdaten in einem ganzzahligen Maschinenformat gespeichertwerden. Sei ein ADC so konfiguriert, dass die Spannungsdifferenz ∆U gemessen werdenkann, ohne in den Uber- oder Unterlauf zu treten. Dann wird das Intervall bei einerVerwendung von n Speicherbits in 2n Intervalle unterteilt. Bei einer Auflosung von 8 Bitwird z.B. in dem Intervall ∆U eine relative Auflosung von 1/28 = 1/256 oder ≈ 0, 39%erreicht.
1.2.2. Gleitpunktzahlen
Gleitpunktzahlen bestehen aus folgenden zu codierenden Funktionseinheiten:
Ziffern Ziffern
±︷ ︸︸ ︷xxxxxx ·
︷ ︸︸ ︷xxxx 10 ± xxx
↑ ↑ ↑ ↑Vor- Dezimal- Vor- Expo-
zeichen punkt zeichen nent
Eine Gleitpunktzahl wird von links nach rechts wie folgt dargestellt: Vorzeichenbit,Ziffernfolge geteilt durch den Dezimalpunkt, Vorzeichenbit des Exponenten, Exponent.
8
1.2. Zahlen, Operationen und elementare Funktionen am Computer
Abbildung 1.2.: Kennlinie fur einen AD-Umsetzer mit einer Auflosung von 2 Bit. DieAuflosung betragt 1/22 = 25%
Die Genauigkeit der Zahl erhoht sich mit der Anzahl der dargestellten Ziffern. Der Ex-ponent, eine binar gespeicherte Zahl z.B. [-64,63], gibt an, mit welcher Potenz einerBasiszahl (i.d.R. 10) die vorliegende Zahl zu multiplizieren ist.
In der Regel (Norm IEEE 754) werden Gleitpunktzahlen im Computer codiert zurBasis zwei gespeichert. Die codierte Zahl z ergibt sich wie folgt aus den n gespeichertenBits:
z = m1 · 2−1 +m2 · 2−2 + ...+mn · 2−n (1.7)
Abbildung 1.3.: Speicherform einer Gleitkommazahl. WIKIPEDIA
In Abbildung 1.3 ist die Computer-Reprasentation einer Gleitkommazahl bestehendaus einem Vorzeichenbit, 8 Bits fur den Exponenten und 23 Bits fur die Mantisse dar-
9
1. Numerische Grundlagen
Bits Vorzeichen Exponent Mantisse Zahlenbereich Zahlentyp
32 1 8 23 1.4 · 10−45...3.40... · 1038 float64 1 11 52 4.9 · 10−324...1.79... · 10308 double
Tabelle 1.2.: Generische Deklarationen von Gleitpunktzahlen.
gestellt.Wie am Beispiel der Zahl 0.1 leicht gezeigt werden kann, geht die Umrechnung einer
beliebigen Dezimalzahl in dezimale Gleitpunktzahlen nicht immer auf. Fur 0.1 konnennur eine kleinste obere und eine großte untere Gleitpunktzahl (float) geschrieben werden.
Darstellungen der Zahl 0.1:
V E(8 Bit) M(23 Bit)V e1e2e3...e8 m1m2m3...............m21m22m23
0 0111.1011 1001.1001.1001.1001.1001.100 = 0, 09999999400 0111.1011 1001.1001.1001.1001.1001.101 = 0.1000000015
Bereits bei der Darstellung einer Zahl ergeben sich somit Rundungsfehler.
Uberlauf, Unterlauf
Das Verhalten eines Programmes bei Uber- oder Unterschreitung des erlaubten Zahlen-bereiches hangt von der Programmiersprache und Compilereinstellungen ab. Bei Unter-lauf werden Zahlen i.d.R. gleich Null gesetzt, bei Uberlauf wird NaN (Not a Number)oder auch Inf (Infinity) codiert. Beide Ereignisse sind fur eine exakte Berechnung vonErgebnissen naturlich unerwunscht. Das Problem muss daher vor Berechnung numerischin angemessener Weise skaliert werden.
Darstellung reeller Zahlen, Rundung
Betrachte Maschinenzahlen der Form:
z = 0.x1...xLBe
mit L als Mantissenlange, B als Basis und e als Exponent
Einer solchen Zahl z wird als Wert zugeordnet:
z =L∑i=1
xiBe−i = x1B
e−1 + x2Be−2 + ...+ xLB
e−L
= BL−e(x1BL−1 + x2B
L−2 + ...+ xL)
10
1.2. Zahlen, Operationen und elementare Funktionen am Computer
Zur Ermittlung der Wertes, muss ein Polynom an der Stelle B ausgewertet werden. DieKonvertierung zwischen Zahlen zur Basis B und Dezimalzahlen kann mit dem Horner-Schema erfolgen. Offenbar gibt es zwischen darstellbaren Maschinenzahlen immer nicht-darstellbare reelle Zahlen.
Als Rundung bezeichnet man die Suche nach einer nahegelegenen Maschinenzahl. Manunterscheidet zwei Typen von Rundungen:
• Die Optimale Rundung: nachstgelegene Maschinenzahl. Sind zwei Zahlen gleichweit entfernt, wird aus statistischen Grunden diejenige mit xL gerade genommen.
• Die Rundung durch Abschneiden: Die Ziffern nach der L-ten Stelle werden wegge-lassen.
Beispiel Optimal Abschneidenx x xx1 = .123456105 .123105 .123105x2 = .56789105 .568105 .567105x3 = .123500105 .124105 .123105x4 = .234500105 .234105 .234105
Definition:
Eine Rundung heißt korrekt, falls zwischen x und x keine Maschinenzahl liegt. BeideRundungstypen sind korrekt.
Der relative Fehler lasst sich wie folgt abschatzen:
Nimmt man an, dass der Exponentenbereich unbeschrankt ist, dann gilt fur einenkorrekt gerundeten Wert:
|x− x||x|
≤ ε = B1−L x 6= 0
Der Rundungsfehler ist durch die maschinenabhangige Zahl ε beschrankt.
1.2.3. Operationen und Funktionen
In Abhangigkeit von den Operationen, die fur eine Berechnung benutzt werden, wird dernumerische Fehler der Berechnung auf verschiedene Weise bestimmt:
11
1. Numerische Grundlagen
• Fur die zweistelligen Operationen ◦ ∈ {+,−, ∗, /} gilt bei korrekter Rundung:
|x ◦ y − ˜x ◦ y||x ◦ y|
≤ ε = B1−L x 6= 0
B sei hier die Basis, L sei die Mantissenlange, es liege kein Uber- oder Unterlaufvor.
• Bei Berechnung einer Potenz x ∗ ∗y oder x∧y gilt:Ist y klein und ganzzahlig (=2, 3, ...) kann der Wert durch Ausmultiplizierenbestimmt werden. Ansonsten wird
x∧y = exp(y · ln(x))
berechnet. Der relative Fehler ist i.a. großer als bei zweistelligen Operationen
∆f
f= c · ε mit c > 1.
• Bei anderen elementaren Funktionen werden Approximationsverfahren zur Be-stimmung des Funktionswertes eingesetzt. Diese Verfahren werden im Folgendenzunachst erlautert, dann wird der relative Fehler des Verfahrens diskutiert.
Die Approximationsverfahren werden nach folgendem Schema durchgefuhrt:
x→ Argumentreduktion → Approximantion → Ergebnisanpassung → f(x).
In der Argumentreduktion wird mit Hilfsformeln die Berechnung eines Funktions-wertes auf einen kleinen Argumentbereich zuruckgefuhrt. Dies wird bei der Ergebnis-anpassung nach der approximativen Berechnung des Funktionswertes wieder ruckgangiggemacht.
Die Approximation kann mit verschiedenen Verfahren erfolgen, von denen mit derKettenbruchdarstellung meist am schnellsten ein hinreichend genaues Ergebnis erzieltwird:
• Kettenbruchdarstellung
• Polynomapproximation
• Potenzreihenentwicklung
• Iterationsverfahren
12
1.2. Zahlen, Operationen und elementare Funktionen am Computer
In einigen Situationen ist auch eine Ruckfuhrung auf andere Funktionen moglich(z.B. cos(x) = sin(π/2− x)).
Als Beispiel betrachten wir die Wurzelfunktion:
√x = sqrt(x) fur x = mBe
mit der Basis B, dem Exponenten e und der Mantisse m ∈ [1/B, 1] .
• Zur Argumentreduktion und Ergebnisanpassung schreibt man
√x =√x0 ·BS mit
{x0 = m, S = e/2, fur e geradex0 = m
B , S = (e+ 1)/2, fur e ungerade
Dabei ist
x0 ∈[1/B2, 1
]Die Funktion muss nur noch fur Argumente aus dem Intervall bestimmt werden.Die Ergebnisse erhalt man aus der Multiplikation mit BS .
• Fur die Aproximation gibt es drei Moglichkeiten:
a) Entwicklung der Funktion in eine Potenzreihe, z.B. die Funktion√
1− z ineine Potenzreihe um 1:
√1− z = 1− 1
2z − 1
8z2 − 1
16z3 − 5
128z4 − ...
Der Konvergenzradius ist hier gleich 1. Die Konvergenz erfolgt aber nur furx ≈ 1(z ≈ 0) genugend schnell. Fur kleine x (z.B. x0 ≈ 0.01) wird derAlgorithmus langsam und ist somit fur praktische Zwecke nicht tauglich.
b) In dem Iterationsverfahren wird die Wurzel
w =√x
zu gegebenen x > 0 gesucht. Fur w gilt w2 = x oder w = xw .
Es sei w0 eine Anfangsnaherung fur√x. Man definiert dann
w0 =x
w0.
13
1. Numerische Grundlagen
Falls w0 = w0 gilt, ist die Iteration abgeschlossen3. Sonst schreibt man
w0 =w2
w0= w
w
w0< w
(fur den Fall w0 > w).Es folgt ein neuer Naherungschritt mit w1 = w0+w0
2 usw.
Das Verfahren wird abgebrochen, wenn die gewunschte Genauigkeit
|w0 + w0| < 10−x
erreicht ist.
Wir benutzen somit folgendes Iterationsverfahren:
w0 > 0 (1.8)
wi = x/wi (1.9)
wi+1 = (wi + wi)/2 (1.10)
Wir testen das Verfahren auf einem Taschenrechner mit B=10 und L=12x = 0.01, w0 = 1
i 0 1 2 3 4 ... 7
wi 1 0.505... 0.2624... 0.1502... 0.1084 ... 0.1
Bereits nach wenigen Iterationsschritten ist der Algorithmus quadratisch kon-vergiert.
c) Kettenbruchentwicklung mit optimalen Koeffizienten fur das Intervall [0.01,1]
Als Ansatz betrachtet man den Kettenbruch
w ∗ (x) = t2x+ t1 +t0
x+ s0
Nun bestimmt man die Koeffizienten so, dass
supx∈[0.01,1]
|√x− w∗(x)|
minimal wird.
Es ergibt sich t2 = 0.5881229, t1 = 0.467975327625; t0 = −0.0409162391674und s0 = 0.099998. Fur alle x0 ∈ [0.01, 1] ist der relative Fehler von w∗ kleinerals 0.02. Mit drei nachtraglichen Iterationen (Methode b) und 14 Rechenope-rationen werden relative Fehler von < 10−5 erreicht.
3Allein wegen der dargelegten Darstellungs- und Rundungsprobleme rechne man nicht damit, dassdieser Fall eintritt.
14
1.2. Zahlen, Operationen und elementare Funktionen am Computer
• In der Ergebnisanpassung wird die Argumentreduktion ruckgangig gemacht.
Außer in der Nahe von Nullstellen und Polen gilt fur die relative Abweichung:
f(x)− f(x)
f(x)≤ cf · ε
mit ε = B1−L, cf hangt von der Approximation und der Argumentreduktion ab.
15
1. Numerische Grundlagen
1.3. Stabilitat
Wie gezeigt, fuhrt die beschrankte Ziffernzahl im Rechner zu einer
• ungenauen Speicherung von Zahlen und einer
• ungenauen Ausfuhrung von Operationen.
Jedoch ist der relative Fehler in einer Operation beschrankt.Das gilt bei einer langerer Folge von Operationen nicht mehr.
Man betrachte zur Motivation die Funktion f(x) in den beiden Schreibweisen
• a) f(x) = (1− x)6
• b) f(x) = 1− 6x+ 15x2 − 20x3 + 15x4 − 6x5 + x6
In Fall (a) ist die Berechnung der Funktion stabil, in Fall (b) ist die Berechnung nume-risch instabil und fuhrt zu großen Fehlern.
1.3.1. Entwicklung von Kriterien fur die numerische Stabilitat
Im Folgenden werden anhand von Beispielen Kriterien fur die numerische Stabilitateines Algorithmus erarbeitet.
• a) f(x) = (x3 + 13)− (x3 − 1
3), ⇒ ∀x : f(x) = 23
x f(x)1 0, 666.666.666...
103 0, 666.666.663...109 0, 663...1011 0
Man beobachtet eine Abnahme der Genauigkeit fur großer werdendes x, also beiDifferenzbildung von großen Zahlen.
• b) f(x) = ((3 + x3
3 )− (3− x3
3 ))/x3, ⇒ ∀x : f(x) = 23
x f(x)
10−1
0, 666.666.667...10−3 0, 666.67...10−6 0
Man beobachtet eine Abnahme der Genauigkeit fur kleiner werdende x durchSummen- bzw. Differenzbildung zwischen zwei etwa gleich großen Zahlen. Diesfuhrt zu einer Ausloschung der fuhrenden Ziffern und einer Verstarkung der rela-tiven Fehler.
16
1.3. Stabilitat
• c) f(x) = sin2(x)1−cos2(x)
, ∀x : f(x) = 1.
Fur x → 0 erfolgt eine Abnahme der Genauigkeit wegen der Division durch einekleine Zahl, die aus Subtraktion von gleich großen Zahlen entstanden ist.
• d) f(x) = sin(x)√1−sin2(x)
, ∀x : f(x) = tg(x).
Die Berechnung wird in der Nahe des Pols bei 90◦ instabil.
• e) Betrachte f(x) = ex2
3 − 1.
Ein Vergleich mit Anfang der Reihenentwicklung (= x2
3 + x4
8 )... zeigt eine Abnahmeder Genauigkeit mit wachsendem x.
• f) Instruktiv ist auch die Progammierung der Formel fur die Standardabweichungeiner Gaußverteilung und deren Test mit konstanten Messwerten:
σn =
√√√√√ 1
n
n∑i=1
x2i −
1
n
(n∑i=1
xi
)2
Fur xi = x =const. gilt σn = 0.
x σ10 σ20
100/3 1.204 · 10−3 1.131 · 103
1000/29 8.062 · 10−4 0↑
Negative Wurzelwird 0 gesetzt
Folgende Schritte sollten den Beispielen nach bei der Codierung einer Formel vermie-den werden:
a) Die Subtraktion von etwa gleichgroßen Zahlen wegen der Ausloschung fuhrender Zif-fern und der Verstarkung der relativen Fehler.
b) Eine Division durch eine kleine Zahl.
c) Eine Multiplikation mit einer großen Zahl.Letztere ((b) und (c)) fuhren zu einer Verstarkung des absoluten Fehlers.
17
1. Numerische Grundlagen
Bemerkungen
1) Bei Ausloschung wird die Differenz i. allg. fehlerfrei berechnet, die Instabilitatkommt von der Vergroßerung von vorher akkumulierten Fehlern.
Beispiele
a) B=10, L=7
0.2789014 · 103
−0.2788876 · 103
Differenz: 0.0000138 · 103 keine Rundungsfehlernormalisiert: 1.38 · 10−2
b) B=10, L=6
Rel. Fehleropt. gerundet 0.278901 · 103 < 2 · 10−6
opt. gerundet −0.278888 · 103 < 2 · 10−6
Differenz: 0.000013 · 103
normalisiert: 1.3 · 10−2 > 5 · 10−2
In diesem Fall ist bei der Subtraktion keine Rundung notig, der relative Fehler desEndergebnisses ist jedoch 25000 mal so groß wie der Eingangsfehler.
2) Trotz Division durch eine kleine Zahl sind folgende Ausdrucke stabil:
f(x) =
{1, fur x = 0
sin(x)x , fur x 6= 0
bei x→ 0
und
f(x) =
{1, fur x = 1x−1ln(x) , fur x 6= 1
bei x→ 1
Sowohl im Zahler als auch im Nenner ist jeweils nur eine Operation auszufuhren. DieStabilitat folgt aus der Gestalt der Schranken fur den relativen Fehler.
Stabilisierung instabiler Ausdrucke
Zur Stabilisierung instabiler Ausdrucke mogen die folgenden Beispiele und Kochrezep-te hilfreich beitragen:
a)√x+ 1−
√x = 1√
x+1+√x
↑ ↑Instabil fur grosse x Stabil fur grosse x.
18
1.3. Stabilitat
√x+ 1−
√x =
(x+ 1)− x√x+ 1 +
√x
=1√
x+ 1 +√x
b) 1− cos(x) = sin2(x)1+cos(x) = 2 sin(x2 )
↑ ↑Instabil fur x→ 0 Stabil fur x→ 0.
In beiden Fallen wird die Differenz der Funktionswerte durch Erweitern so umgeformt,dass eine analytisch auswertbare Differenz im Nenner steht. Die Differenz kann dannrundungsfehlerfrei berechnet werden.
c) Die Funktion ex − 1 ist instabil fur x→ 0.
Man substituiert y = ex underhalt ex − 1 = y − 1, das weiterhin instabil ist.
Fur x 6= 0 gilt
ex − 1 = (y − 1)x
x=
(y − 1) · xln(y)
,
und dadurch
ex − 1 =
{x, fur y = 1
y−1ln(y)x, fur y 6= 1
mit y → ex.
d) Die Funktion ex − 1− x ist instabil fur x→ 0.
Zur Stabilisierung entwickelt man ex in eine Potenzreihe
ex =n∑i=1
xk
k!.
Dann ist
ex − 1− x =
{ex − 1− x, fur |x| ≥ cx2
2 + x3
6 + ..., fur |x| ≤ c
Die Konstante c ist durch Ausprobieren und Fehlerabschatzung geeignet zu wahlen.
19
1. Numerische Grundlagen
e) Mittelwert und Standardabweichung
Mittelwert sn und Standardabweichung σn fur die Messwerte x1, ..., xn sind wie folgtdefiniert:
sn =1
n
n∑i=1
xi
σ =√tn/n oder
√tn/(n− 1), tn =
n∑i=1
(xi − sn)2
Eine direkte Berechnung der Großen nach diesen Formeln ist unpraktisch, weil alle xibei der Berechnung gespeichert werden mussen.
Naive Umformung:
In den vielen Statistikbuchern wird, um das Zwischenspeichern zu vermeiden, dieuntenstehende Formel zur Berechnung der Standardabweichung empfohlen.
tn =n∑i=1
x2i − 2sn
n∑i=1
xi + s2n
n∑i=1
1
=n∑i=1
x2i − ns2
n
=n∑i=1
x2i − 1
n
(n∑i=1
xi
)2
Der Vorteil dieser Schreibweise besteht darin, dass tn ohne Speicherung der xi berech-net wird. Als Nachteil fuhrt die Rechenmethode zu Instabilitaten wegen Ausloschung,da i.a.
n∑i=1
x2i ≈ ns2
n
Zur Umformung in einen stabilen Ausdruck benutzt man folgende Idee: Beim Hin-zufugen eines neuen Messwertes andern sich sn und tn nur wenig. Deshalb betrachtenwir die Differenzen sn − sn1 und tn − tn1 :
sn − sn1 = (n−1)sn−1+xnn − sn−1
= xn−sn−1
n = δnn ,
δn : = xn − sn−1,
20
1.3. Stabilitat
tn − tn1 =
(n∑i=1
x2i − ns2
n
)−(n−1∑i=1
x2i − (n− 1)s2
n
)= x2
n − ns2n + (n− 1)s2
n−1
= (δn + sn−1)2 − n(sn−1 + δn
n
)2+ (n− 1)s2
n−1
= δn(δn + δn
n
)= δn [(xn − sn−1)− (sn − sn−1)]= δn(xn − sn).
Man erhalt so neue Rekursionsformeln fur die Berechnung von sn und tn:
s1 = x1, t1 = 0δi = xi − si−1, i ≥ 2
si = si−1 + δii , i ≥ 2
ti = ti−1 + δi(xi − si), i ≥ 2
mit tn =√tn/n bzw. tn =
√tn−1/n.
Die Differenzen
xi − si−1 und xi − si
sind harmlos, weil die mogliche Ausloschung keine große Verstarkung des relativenFehlers bewirken kann, da die Differenz mit einer kleinen Zahl δi multipliziert und dannzu einer i.a. großeren Zahl ti−1 addiert wird.
f) Die Losung einer quadratischen Gleichung
ax2 + bx+ c = 0
lautet:
x1,2 =−b±
√b2 − 4ac
2a
Ein so programmierter Zusammenhang ist instabil fur b2 � 4ac, wenn die Wurzel undb das gleiche Vorzeichen haben. Eine Umformung liefert:
x1,2 =2c
−b∓√b2 − 4ac
.
was jedoch dazu fuhrt, dass der Term instabil wird, wenn die Wurzel und b das ent-gegengesetzte Vorzeichen haben. Eine sinnvolle Kombination beider Schreibweisen ist
21
1. Numerische Grundlagen
q := −(b · sign(b)
√b2 − 4ac
)/2
mit
x1 =q
a, x2 =
c
q.
1.4. Fehlerfortpflanzung und Kondition
Die oben betrachtete Stabilitat machtAussagen uber den Einfluss von Run-dungsfehlern bei ungenauer Rechnung.
Im Gegensatz dazu beschreibt die Kondition die Fortpflanzung von Anfangsfeh-lern bei genauer Rechnung.
Betrachten wir dazu das Beispiel:
f(x) =1
1− xx = 0, 999 ⇒ f(x) = 1000.
Man fuhre nun eine analytische Fehleranalyse fur x = 0, 999 + ε mit ε klein aus.
f(x) =1000
1− 1000ε= 1000(1 + 103ε+ 106ε2 + ...)
Der relative Fehler ist dann:
|x− x|x
< 1.1ε und|f(x)− f(x)|
f(x)= 103ε+O(ε2).
Unabhangig von der numerische Methode wird der relative Fehler um einen FaktorO(1000) vergroßert.
Solche Probleme bezeichnet man als schlecht konditioniert.
Ein quantitatives Maß fur die Kondition von differenzierbaren Funktionen ist die Kon-ditionszahl K, die den Verstarkungsfaktor des relativen Fehlers angibt.
22
1.4. Fehlerfortpflanzung und Kondition
Sei x eine Naherung von x mit dem relativen Fehler
ε =x− xx
bzw. x = x(1 + ε).
Entwickelt man f(x) in einer Taylor-Reihe:
f(x) = f(x+ εx) = f(x) + εxf ′(x) +O(ε2))
f(x) =(
1 + xf′(x)f(x) ε+O(ε2)
)so gilt fur den relativen Fehler von f :
|f(x)− f(x)||f(x)|
=
∣∣∣∣xf ′(x)
f(x)
∣∣∣∣ · |ε|+O(ε2) = K · |ε|+O(ε2).
Im Gegensatz zur Stabilitat kann die Konditionszahl exakt bestimmt werden.
K :=
∣∣∣∣xf ′(x)
f(x)
∣∣∣∣Es gilt bei Vernachlassigung der hoheren Glieder:
|f(x)− f(x)||f(x)|
= K|x− x|x
,
wobei man folgende Falle fur K unterscheidet:
• K < 1 Fehlerdampfung;
• K > 1 Fehlerverstarkung;
• K � 1 Problem schlecht konditioniert.
In einer eindimensionalen Konditionsanalyse gelten folgende Zusammenhange:
a) Falls an einer Stelle f ′(x∗) 6= 0 die Funktion f(x) → 0 fur x → x∗ 6= 0 geht,dann strebt K → ∞ fur x → x∗. Mit anderen Worten: f ist in der Nahe von einfachenNullstellen 6= 0 schlecht konditioniert .
b) Sei f(x) = (x− x∗)mg(x) bei g(x∗) 6= 0 und m 6= 0. Dann ist fur m > 0 bei x∗ eineNullstelle m-ter Ordnung und fur m < 0 ein x∗ Pol m-ter Ordnung.
Es gilt weiter:
23
1. Numerische Grundlagen
f ′(x) = m(x− x∗)m−1g(x) + (x− x∗)mg′(x), und wir erhalten
K = x|f ′(x)||f(x)|
= |x| ·∣∣∣∣ m
x− x∗+g′(x)
g(x)
∣∣∣∣ = |m| ·∣∣∣∣x− x∗x
∣∣∣∣−1
+ ...
Fur x→ x∗ ist
K =
{∞ falls x∗ 6= 0|m| falls x∗ = 0
In der Nahe von Polstellen bzw. Nullstellen x∗ 6= 0 ist die Kondition schlecht, namlichgenau umgekehrt proportional zu
∣∣x−x∗x
∣∣. Harmlos sind dagegen Pol- bzw. Nullstellenx∗ = 0, da hier die Konditionszahl nur etwa die Ordnung des Pos bzw. der Nullstelleangibt.
c) Falls f ′(x) einen Pol bei x∗ hat, ist die Kondition bei x∗ ebenfalls schlecht.
Betrachte z.B. f(x) = 1 +√x− 1. Diese Funktion hat die Konditionszahl
K =
∣∣∣∣x2(
1 +1√x− 1
)∣∣∣∣ ,und K →∞ fur x→ 1.
1.4.1. Vergleich von Kondition und Stabilitat
Zur Illustration der Tatsache, dass Kondition und Stabilitat nicht mit einander kor-relieren mussen, werde folgende Funktion betrachtet:
f(x) =
√1
x− 1−
√1
x+ 1
fur
0 < x < 1.
Eine Untersuchung der Stabilitat liefert fur:
x→ 0 eine Ausloschung, die zur Instabilitat fuhrt.
x→ 1 ein stabiles Verhalten des Ausdrucks.
Berechnet man nun die Kondition, so erhalt man:
f ′(x) = −1/x2
2√
1x−1− −1/x2
2√
1x
+1=
√1x−1−
√1x
+1
2x2√
1x−1
√1x
+1
24
1.4. Fehlerfortpflanzung und Kondition
und
K = x|f ′(x)||f(x)|
=1
2√
1− x2
Im Ergebnis ist das Problem fur
x→ 0 gut konditioniert, da K = 12
und fur
x→ 1 schlecht konditioniert, da K =∞.
25
1. Numerische Grundlagen
Fragen zur Selbstkontrolle
• Was sind arithmetische Ausrucke? Aus welchen Bestandteilen bestehen sie?
• Wie werden welche Typen von Zahlen dargestellt? Konnen alle Zahlen dargestelltwerden? In welchem Typ werden Daten zunachst aufgezeichnet? Welche Formender Rundung gibt es?
• Wie wird ein ADC dimensioniert, wenn bekannt ist, mit welcher Genauigkeit ge-messen werden soll?
• Was sind zweistellige Operationen? Wie genau werden sie berechnet?
• Mit welchen Schritten konnen elementare Funktionen berechnet werden? Wie wirddie Genauigkeit dieser Berechnung bestimmt?
• Erlautere die Begriffe Stabilitat und Kondition und vergleiche sie.
• Welche Stabilisierungsregeln gibt es?
• Was ist eine Konditionsanalyse?
26
2. Wahrscheinlichkeit und Generatoren
2.1. Vorbemerkungen
Grundsatzlich konnen in der Physik zwei Typen von Experimenten durchgefuhrt werden.Entweder sollen Parameter mit moglichst hoher Genauigkeit gemessen werden (parame-ter determination/estimation) oder die Gultigkeit von Hypothesen soll getestet werden(hypothesis testing). Der Ubergang zwischen diesen Typen kann im Experiment fließendsein. (Messung eines Parameters bzw. Test, ob dieser Parameter mit einer Vorhersageubereinstimmt.)Die Wahrscheinlichkeit kann bei diesen Analysen in drei unterschiedlichen Formen rele-vant werden. Entweder ist sie durch mathematische Berechnung quasi a priori gegebenoder sie muss experimentell aus einer moglichst großen Zahl von Messungen ermitteltwerden - schließlich kann sie, wenn das Experiment etwa aus Kostengrunden bzw. we-gen der Einmaligkeit der Situation (z.B. Supernova-Explosion) nicht wiederholt werdenkann, abgeschatzt werden. Diese sogenannten Bayesischen Ansatze, die darauf beruhenunbekannte Wahrscheinlichkeiten abzuschatzen bzw.
”begrundet zu erraten“ werden in
Kapitel ?? diskutiert.
2.1.1. Definitionen
Der Begriff der Wahrscheinlichkeit kann abhangig davon, ob a priori Wissen uber denbetrachteten Vorgang zur Definition benutzt werden kann oder nicht, auf zwei unter-schiedliche Weisen eingefuhrt werden:
• 1) Falls ein Ereignis auf n verschiedene und gleich wahrscheinliche Arten eintretenkann und k davon die Eigenschaft A haben, so ist die Wahrscheinlichkeit von dasAuftreten von A
P (A) =k
n=
gunstige
moglicheFalle.
• 2) Besitzt man kein a priori Wissen uber die Eigenschaften des Zufallsexperi-ments, kann man wie folgt empirisch vorgehen: Die Eigenschaften A und nicht-Aeines Experimentes werden n-fach unabhangig beobachtet. Dabei trete k mal dieEigenschaft A auf. Dann ist die Wahrscheinlichkeit P (A) gegeben durch
P (A) = limn→∞
k
n.
27
2. Wahrscheinlichkeit und Generatoren
Beispiel: Bei einem Munzwurf ist der Ausgang Kopf oder Zahl moglich. (Wie hangenbeide Definitionen der Wahrscheinlichkeit zusammen?) Das Experiment bestehe aus zweiWurfen. Der erste Wurf ergebe Kopf. Wie groß ist die Wahrscheinlichkeit noch einmalKopf zu werfen?
2.1.2. Kombination von Wahrscheinlichkeiten
• Gegeben seien die Ereignistypen A und B mit den Wahrscheinlichkeiten P (A) undP (B), dann ist die Wahrscheinlichkeit fur A oder B
P (A ∨B) = P (A) + P (B)− P (A ∧B).
Wenn sich A und B ausschließen, gilt P (A ∧B) = 0 und
P (A ∨B) = P (A) + P (B).
Als Spezialfall sei: B = A (nicht A), dann ist
P (A ∨A) = P (A) + P (A) = 1.
• Gegeben seien die Ereignistypen A und B mit den Wahrscheinlichkeiten P (A) undP (B), dann ist die Wahrscheinlichkeit fur A und B
P (A ∧B) = P (A) · P (B|A).
P (B|A) ist dabei die bedingte Wahrscheinlichkeit, dafur dass B auftritt, wenn Aeingetreten ist.1
Sind A und B unabhangig, dann ist
P (B|A) = P (B)
und somitP (A ∧B) = P (A) · P (B).
1lies: P fur B gegeben A
28
2.2. Zufallsvariable und deren Verteilung
2.2. Zufallsvariable und deren Verteilung
Ziel ist zunachst die Klassifizierung von moglichen Endzustanden eines statistischenVorganges; dann sollen Methoden angegeben werden, mit denen beliebige Verteilungenbeschrieben werden konnen.
Beispiel zur Klassifizierung: Bei einem Munzwurf werden z.B. werden die Zuordnun-gen Kopf → 0 und Zahl → 1. vorgenommen.
Allgemein: Wird dem Ereignis Ai die ganze Zahl i zugewiesen, erhalt man eine diskreteZufallsvariable wie z.B. die Zahl der Teilchen in einem Detektor oder die Augen einesWurfels i.
Kontinuierliche Zufallsvariable werden genutzt, wenn es nicht moglich ist, die Ereignis-se ganzen Zahlen zuzuordnen, wie z.B. bei kontinuierlichen physikalischen Verteilungen(Winkelverteilung, Energiespektrum, usw.).
Wir suchen nun nach einer Beschreibung des moglichen Ausgangs von Zufallsexperi-menten. Die Zufallsvariable r moge den moglichen Ausgang des Experimentes angeben.Sie wird mit der reellen Zahl x verglichen, die jeden Wert zwischen −∞ und +∞ anneh-men kann. Gesucht ist die Wahrscheinlichkeit dafur, dass ein Ereignis eintritt, bei demdie Zufallsvariable r kleiner ist als ein vorher gewahltes x (r < x). Dazu bildet man dieVerteilungsfunktion:
F (x) = P (r < x),
die die Summe aller Ereignisse unterhalb von x normiert auf die Gesamtzahl der Ver-suche angibt.
Fur einen Wurfel, bei dem die Zahl der Augen r sechs diskrete Werte annehmen kann,ergibt sich fur F (r) eine sechsstufige Treppenfunktion, die monoton und nicht-fallendvon 0 auf 1 ansteigt.
Im Grenzfall einer kontinuierlichen Verteilung ist
limx→∞
F (x) = limx→∞
P (r < x) = 1.
Da die Summe aus P (A) + P (A) = 1 ist, gilt
P (r ≥ x) = 1− F (x) = 1− P (r < x).
Somit istlim
x→−∞F (x) = lim
x→−∞P (r < x) = 1− lim
x→−∞P (r ≥ x) = 0.
Wenn die Verteilungsfunktion stetig differenzierbar ist, gilt
dF (x)
dx= F ′(x) = f(x),
29
2. Wahrscheinlichkeit und Generatoren
0 1 2 3 4 5 6 70
0.2
0.4
0.6
0.8
1
x
F(x)
Abbildung 2.1.: Verteilungsfunktion eines Wurfels.
f(x) heißt dann Wahrscheinlichkeitsdichte von r und gibt ein Maß fur die Wahrschein-lichkeit in dem Intervall x ≤ r ≤ x+ dx an.
Die Wahrscheinlichkeit, dass r kleiner ist als ein vorgewahlter Wert a ist, ist gegebendurch:
P (r < a) =
a∫−∞
f(x)dx = F (a),
die Wahrscheinlichkeit, dass r in einem Intervall zwischen a und b liegt, ist:
P (a ≤ r ≤ b) =
b∫a
f(x)dx = F (b)− F (a).
Insbesondere gilt bei Integration uber den gesamten Bereich in x:
∞∫−∞
f(x)dx = 1.
30
2.3. Allgemeine Eigenschaften einer Zufallsvariablen: Erwartungswert, Streuung, Momente, etc.
0 0.2 0.4 0.6 0.8 10
0.2
0.4
0.6
0.8
1F(x)
x
Abbildung 2.2.: Gleichverteilung: Verteilungsfunktion
0 0.2 0.4 0.6 0.8 1
1
1.2
1.4
1.6
1.8
2F(x)
x
Abbildung 2.3.: Gleichverteilung: Wahrscheinlichkeitsdichte
2.3. Allgemeine Eigenschaften einer Zufallsvariablen:Erwartungswert, Streuung, Momente, etc.
Nachdem die Begriffe der Verteilungsfunktion und Wahrscheinlichkeitsdichte definiertworden sind, sollen nun allgemein (ohne Bezug auf eine spezielle Wahrscheinlichkeits-
31
2. Wahrscheinlichkeit und Generatoren
dichte zu nehmen) die Eigenschaften von Wahrscheinlichkeitsverteilungen beschriebenwerden. Betrachten wir dazu die Zufallsvariable r und die Funktion y = H(r), die selbstauch als Zufallsvariable betrachtet werden soll.
Der Mittelwert r oder Erwartungswert E(r) bei einer diskreten Verteilung von r ist:
x = E(r) =n∑i=1
(xi · P (r = xi)).
Der Erwartungswert einer Funktion von diskreten r ist
E[H(r)] =n∑i=1
(H(xi) · P (r = xi)).
Analog ist der Erwartungswert fur kontinuierlich verteilte r
E(r) = x =
∞∫−∞
x · f(x)dx,
und fur eine Funktion davon
E[H(r)] =
∞∫−∞
H(x) · f(x)dx.
Der Erwartungswert stellt bei einer Messung die beste Schatzung fur den wahren Wertund den Schwerpunkt der Verteilung dar.
Zu den Eigenschaften, durch die eine statistische Verteilung charakterisiert ist, gehorenneben ihrem Schwerpunkt auch ihre Breite und Symmetrie. Dazu betrachten wir als Spe-zialfall die Funktion: H(r) = (r − c)l.
Die Erwartungswerte dieser Funktion
al = E[(r − c)l]
heißen dann die l-ten Momente um den Punkt c.
Berechnen wir nun die Momente µl um den Mittelwert
µl = E[(r − x)l],
dann sind die Momente µ0 = 1 2 und µ1 = 0 3 trivial zu bestimmen.
2∞∫−∞
(x− x)0f(x) dx = 1
3siehe Definition des Mittelwertes
32
2.3. Allgemeine Eigenschaften einer Zufallsvariablen: Erwartungswert, Streuung, Momente, etc.
Fur das zweite Moment gilt
µ2 =
∞∫−∞
(x− x)2f(x) dx, ,
es ist das niedrigste Moment, das etwas uber die mittlere Breite der Verteilung der Ab-weichung von x von dem Mittelwert x aussagt.
Die so definierte Varianz
E[(r − x)2] = σ2(r) = var(r) =
∞∫−∞
(x− x)2f(x) dx,
ist ein Maß fur die Breite der Verteilung. Die Wurzel aus der Varianz σ =√σ2(r) heißt
Streuung, Standardabweichung oder Signifikanz. Da die Standardabweichung die gleicheDimension hat wie die (ggf. gemessene) Zufallsavriable r, wird die einfache Standardab-weichung mit dem Messfehler, σ(r) = ∆x identifiziert.4
Das dritte Momente um den Mittelwert, die Schiefe (= skewness), beschreibt die Sym-metrie der Verteilung. Dimensionslos ist die Schiefe definiert als
γ =µ3
σ3
und gibt die Asymmetrie in der Einheit der Streuung an.
Ist die Schiefe negativ, hat die Verteilung Auslaufer nach links; ist sie positiv, gibt esAuslaufer nach rechts. Je großer der Betrag der Schiefe, desto großer ist die Asymmetrie.Ist der Betrag der Schiefe gleich null, ist die Verteilung symmetrisch.Durch den Quotienten von dem vierten Moment um den Mittelwert und dem Quadratder Varianz wird die curtosis einer Verteilung definiert:
C = µ4/σ4.
C ist groß, wenn die Verteilung uber großere Auslaufer verfugt als die Gauß-Verteilung.Die Gauß-Verteilung selbst liefert c=3.
2.3.1. Regeln uber Mittelwerte und Varianzen
• Betrachten wir zunachst die Multiplikation jeder Zahl einer Verteilung mit (der-selben) Konstanten:
H(x) = cx, c = const.
4Vorsicht: Was dieser Fehler bedeutet, hangt von der speziellen Form der Verteilung ab und wird spaterdiskutiert.
33
2. Wahrscheinlichkeit und Generatoren
Es folgt, dassE(c · r) = c · E(r), und σ2(c · r) = c2 · σ2(r).
Daher istσ2(r) = E[(r − x)2] = E[r2 − 2rx+ x2] = E(r2)− x2.
• Die reduzierte Variable
u =r − xσ(r)
,
hat den Erwartungswert
E(u) =1
σ(r)E(r − x) =
1
σ(x)(x− x) = 0,
und die Varianz
σ2(u) =1
σ2(x)E[(r − x)2] =
σ2(x)
σ2(x)= 1.
• Der wahrscheinlichste Wert einer Verteilung ist jener Wert von x, bei dem
P (x = xm) = maximal.
Ist die Wahrscheinlichkeitsdichte differenzierbar, berechnet man
d
dxf(x) = 0; und testet, ob
d2
dx2f(x) < 0 ist
Besitzt die Verteilung ein Maximum, heißt sie unimodal, sonst heißt sie multimodal.
• Der Median ist derjenige Wert einer Verteilung, fur den die VerteilungsfunktionF = 1
2 ist
F (x1/2) = P (r < x1/2) =1
2.
Ist f(x) stetig, giltx1/2∫−∞
f(x)dx =1
2.
Ist die Verteilung unimodal, stetig und symmetrisch, dann ist Erwartungswertgleich dem wahrscheinlichsten Wert oder Median.
• Der quadratische Mittelwert root mean square = RMS ist definiert als
xrms =√E(x2) =
√σ2(x) + x2.
Ist der Erwartungswert gleich null (x = 0), dann ist xrms = σ(x).
34
2.3. Allgemeine Eigenschaften einer Zufallsvariablen: Erwartungswert, Streuung, Momente, etc.
0 2 4 6 8 100
0.5
1
1.5
2
2.5
3
3.5
4F(x)
x
Abbildung 2.4.: Asymmetrische Wahrscheinlichkeitsdichte (ahnlich Maxwellscher Ge-schwindigkeitsverteilung) Noch einzutragen: verschiedene Mittelwerte.
Betrachten wir zur Illustration eine Funktion, die der Maxwellsche Geschwindig-keitsverteilung von Teilchen eines idealen Gases:
f(ν) = N ·( m
2πkT
)3/2· e
mν2
2kT · 4πν2
ahnelt.
• Haufig interessiert man sich dafur, in welchem Intervall um den Mittelwert einerVerteilung ein bestimmter Prozentsatz der Zufallszahlen (zufallige Meßergebnisse)liegt. Oder man mochte berechnen, mit welcher Wahrscheinlichkeit eine Vorhersa-ge der theoretischen Physik durch eine (zufallsverteilte) Messung ausgeschlossenwerden kann. Dazu wird der Begriff des Quantils benotigt:
Das Quartil einer Verteilung ist analog zu x1/2 definiert als:
F (x1/4) = 0.25, F (x3/4) = 0.75,
unteres Quartil oberes Quartil
Entsprechend sind Dezile (q=10%) und Quantile (q=beliebige Prozentsatze) defi-niert als:
F (Xq) =
Xq∫−∞
f(x)dx = q.
35
2. Wahrscheinlichkeit und Generatoren
2.4. Gleichverteilung
Bisher haben wir uns das Werkzeug beschafft, um einfache statistische Zusammenhangezu verstehen und (z. B. gemessene) Verteilungen einfach zu beschreiben. In diesem Ab-schnitt werden die Eigenschaften der Gleichverteilung als spezieller Wahrscheinlichkeits-verteilung diskutiert. Die Gleichverteilung ist deshalb besonders wichtig, weil die nume-rische Erzeugung von gleichverteilten Zahlen die Grundlage fur die Erzeugung beliebigverteilter Zahlen ist.
Gegeben sei die Wahrscheinlichkeitsdichte mit folgenden Eigenschaften:
f(x) = c fur a ≤ x < b
f(x) = 0 fur x < a und x ≥ b.
Aus der Unitaritatsbedingung lasst sich die Konstante c bestimmen:
∞∫−∞
x · dx = c
b∫a
dx = c · (b− a) = 1
f(x) =1
b− afur a ≤ x < b
f(x) = 0 fur x < a und x ≥ b.
Dann ist die Verteilungsfunktion
F (x) =
x∫a
dx
b− a=x− ab− a
fur a ≤ x ≤ b
und
F (x) = 0 fur x < a sowie F (x) = 1 fur x ≥ b.
Der Erwartungswert ist durch
E(x) = x =1
b− a
b∫a
x · dx =1
2
1
b− a(b2 − a2) =
b+ a
2,
und die Varianz durch
σ2(x) =1
12(b− a)2
gegeben.
36
2.5. Erzeugung von gleich- und beliebig verteilten Zufallszahlen auf dem Computer
2.5. Erzeugung von gleich- und beliebig verteiltenZufallszahlen auf dem Computer
Es ist technisch moglich, echt zufallsverteilte Zahlen auf einem Computer zu generieren.Da man aber (nicht nur zur Suche nach Programmierfehlern) auf die Reproduzierbarkeitvon Rechenergebnissen angewiesen ist, muss auch die Folge von erzeugten Zufallszahlenreproduzierbar sein.Man sucht somit nach einer Reihe von Werte, deren Haufigkeit o.B.d.A. in einem Intervallvon 0 bis 1 gleichverteilt ist, deren Abfolge aber festliegt. Der Algorithmus soll also strengdeterministisch sein, die Abfolge der Zahlen ist pseudozufallig.
2.5.1. Lineare kongruente Generatoren (LCG)
In einem Generator von Zufallszahlen, wird in einer Sequenz aus allen j vorher gebildetenPseudo-Zufallszahlen eine neue Zufallszahl berechnet. Es ist also
xj+1 = f(x1, ..., xj).
Wir betrachten den Algorithmus:
xj+1 = ((a · xj + c)mod m)/m
Mit dem Multiplikator a, dem Inkrement c und dem Modulus m; alle drei Variablen sindganzzahlig. Das Ergebnis der Operation (e) mod m ist der Rest bei Division von e durchm (z.B. ist 7 mod 6 = 1).
Zufallsgeneratoren, die auf diesem Algorithmus basieren heißen linear kongruent. Wieman durch Einsetzen von kleinen Zahlen leicht sieht, liefert dieser Algorithmus offen-bar periodisch aber gleichverteilte Zufallszahlen im Abstand von jeweils 1/m-tel. Furpraktische Anwendungen will man bei gegebenem Zahlentyp im Rechner eine moglichstgroße Periodendauer erzielen. Bei der Suche nach geeigneten Werten hilft folgender Satz:
Satz uber die maximale Periode einer LCG mit c 6= 0.
Der LCG wird durch m, a, c, x0 definiert und hat dann die Periode m wenn
• c und m teilerfremd sind;
• b = a− 1 ein Vielfaches von p ist fur jede Primzahl p, die Teiler von m ist;
• b ein Vielfaches von 4 ist, falls m ein Vielfaches von 4 ist.
Nutzlich sind Teilfolgen, die kurz gegen die Periodenlange sind.
Beispiel: LCG-Ergebnisse fur c = 3, a = 5,m = 16, x0 = 0
0, 3, 2, 13, 4, 7, 6, 1, 8, 11, 10, 5, 12, 15, 14, 9, 0...
37
2. Wahrscheinlichkeit und Generatoren
2.5.2. Multiplikativ linear kongruente Generatoren MLCG
Sei in einem LCG gleich c = 0, dann ist
xj+1 = (a · xj) mod m.
Solche Generatoren werden als Multiplikativ Linear Kongruente Generatoren (MLCG)bezeichnet. Sie sind durch folgende Eigenschaften charakterisiert:
• Wegen der fehlenden Addition (im Vergleich zum LCG) ist pro Erzeugung einerZufallszahl eine Operation weniger auszufuhren. Die Rechnung ist schneller.
• Die Zahl 0 kann nicht mehr erzeugt werden
• Die maximale Periode ist kurzer.
Die Ordnung λ(m) eines primitiven Elements a modulus m ist wie folgt definiert:Seien a und m teilerfremd und λ ganzzahlig. Dann ist die Ordnung von m das kleinsteλ, fur das gilt: aλ mod m = 1. a heißt dann primitives Element zur Ordnung λ.
Sind z.B. a = 4,m = 7, dann ist
4λ mod 741 = 4 0 · 7 = 0 442 = 16 2 · 7 = 14 243 = 64 9 · 7 = 63 1
⇒ λ = 3.
Die maximale Periode eines MLCG, der durch m, a, c = 0 und x0 definiert ist, istgleich der Ordnung λ(m). Sie wird erreicht wenn
• der Multiplikator a ein primitives Element modulo m ist, und
• x0 und m teilerfremd sind.
In der Praxis sind zwei Konsequenzen relevant:
(a) Es sei m = 2l und m − 1 sei die großte auf dem Rechner darstellbare Zahl, dannist die maximal mogliche Periodenlange m
4 .
(b) Ist m = p = Primzahl, dann ist die maximale Periode m− 1.
38
2.5. Erzeugung von gleich- und beliebig verteilten Zufallszahlen auf dem Computer
2.5.3. Spektraltest
Betrachte die Folge 0, 1, 2, 3, ...,m; m sei die großte darstellbare Zahl.
Dann ist die Periodenlange groß und die generierten Zufallszahlen sind gleichervteilt.Aber: Ihre Abfolge ist nicht zufallig.
Spektraltest: Wir benotigen daher einen Test zur Aufspurung von nicht zufalligenAbhangigkeiten zwischen benachbarten Elementen in einer Folge.
Betrachte den MLCG mit a = 3,m = 7, c = 0, x0 = 1, dann ergibt sich sich die Folge
1, 3, 2, 6, 4, 5, 1, ...
Zum Test bildet man Paare aus benachbarten Zahlen (xj , xj+1) mit j = 0, 1, ..., n−1 DiePeriodenlange sei n (hier n = m− 1 = 6). Tragt man die Paare in ein zweidimensionalesKoordinaten-System ein, erkennt man ein spezifisches Muster.
1 2 3 4 5 6
1
2
3
4
5
6F(x)
x
Abbildung 2.5.: Spektraltest: zweidimensionales Diagramm fur den im Text genanntenGenerator.
Wir bemerken:
1. im Wertebereich 1 ≤ x ≤ n gibt es n2 mogliche Zahlenpaare, (n ganzzahlig)
2. davon sind jedoch nur n Moglichkeiten realisiert.
3. Der Gitterabstand ist gleich 1, wir gehen uber zu den transformierten Zahlen:uj = xj/m
39
2. Wahrscheinlichkeit und Generatoren
• Der Gitterabstand ist nun =1/m
• und die Kantenlange=1.
4. Durch die besetzten Punkte lassen sich endlich viele Familien von Geraden legen.
5. Betrachte den Abstand von benachbarten Linien einer Familie (Die Steigungendieser Geraden sind gleich);
6. Ist das Gitter gleichbesetzt ist der Abstand der Linienpaare der minimale realisierteAbstand d2 = m−1/2.
7. Ist das Gitter ungleichmaßig besetzt, dann ist der Abstand d2 � m−1/2
Betrachten wir zur Verallgemeinerung auf n Dimensionen das
n− tupel (uj , uj+1, ..., uj+n).
Die Familien von Geraden fur (n = 2) werden zu (n− 1)-dimensionalen Hyperebenen.Der Abstand von gleichbesetzten Gittern ist etwa dn ≈ m−1/n; von ungleich besetztenGittern ist er dn � m−1/n.
Als geeignete Moduli m und Multiplikation a ergeben sich:
32 Bit m a m a 16 Bit] 2147483647 39373 32749 162 ]] 2147483563 40014 32363 157 ]] 2147483399 40692 32143 160 ]] 2147482811 41546 32119 172 ]] 2147482801 42024 32722 146 ]] 2147482739 45742 32657 142 ]
Als Modulus m wurden Primzahlen nahe der großten darstellbare Zahl getestet.Als Multiplikator a wurden primitive Elemente modulo m mit a <
√m gewahlt.
Bemerkungen zur Praktische Implementation von Generatoren5
• Portabilitat: Zu testen ist, ob Zufallsgeneratoren mit dem gleichen Programmcodeauf anderen Rechnerarchitekturen exakt dieselben Ergebnisse liefern. In der Regeltun sie das nicht.
• Seed: Bei Nutzung von Zufallszahlengeneratoren aus Programmbibliotheken infor-miere man sich uber die Konsequenzen der Wahl der Startparameter (seed). Falschgewahlte Startwerte konnen die Periodendauer stark verkurzen und so artifizielle(nicht physikalisch bedingte) Muster im Output erzeugen.
• Ist die Periodendauer zu kurz, konnen mehrere MLCG miteinander kombiniertwerden.
5Weiterfuhrende Literatur z.B.: Siegmund Brand, Datenanalyse und William Pres, et al., NumericalRecipies
40
2.5. Erzeugung von gleich- und beliebig verteilten Zufallszahlen auf dem Computer
2.5.4. Erzeugung beliebig verteilter Zufallszahlen (Teil 1)
Mit Hilfe eines Generators, der linear verteilte Zufallszahlen erzeugt, konnen mit folgen-der Methode beliebig verteilte Zufallszahlen erzeugt werden.
Es sei xr eine gleichverteilte Zufallsvariable mit der Wahrscheinlichkeitsdichte
f(x) = 1 im Intervall 0 ≤ x < 1; sowie f(x) = 0 fur x < 0, oder x ≥ 1.
Es sei weiter yr eine Zufallsvariable mit der beliebigen (integrierbaren) Wahrschein-lichkeitsdichte g(y).
Beide Verteilungen konnen uber eine Substitution miteinander in Verbindung gebrachtwerden:
g(y) =
∣∣∣∣dxdy∣∣∣∣ · f(x).
Da f(x) = 1, folgt dass g(y)dy = dx. Die Verteilungsfunktion G(y) ergibt sich aus
dG(y)
dy= g(y) durch Umstellung erhalt man dG(y) = g(y)dy = dx
und die Integration liefert
x = G(y) =
y∫−∞
g(t)dt
Invertiert man die Funktion, erhalt man
y = G−1(x)
als Funktion von x. Generiert man eine gleichverteilte Zufallszahl x kann durch Anwen-dung von G−1 eine Folge von Zufallszahlen, die g(x) folgt, erzeugt werden.
Beispiel: Zu generieren seien Zufallszahlen, deren Verteilung im Bereich von 0 bis πder Dichtefunktion f(x) = sin(x) folgt.
Man berechnet zur Normierung zunachst die Flache unter der Kurve in dem Bereich:
A =
∫ π
0sin(x)dx = 2.
Dann die Flache bis zur Zufallsvariablen x:
A(x) =
∫ x
0sin(x)dx = 1− cos(x).
Die auf den Zahlenbereich zwischen 0 und 1 normierte relative Flache r(x) ist somit:
r(x) = A(x)/A = (1− cos(x))/2.
Der gesuchte Generator ergibt sich durch Invertieren der Funktion zu
x(r) = arccos(2 · (1− r)).
Die Zufallszahl r kann in einem Generator fur linear verteilte Zufallszahlen erzeugt wer-den.
41
2. Wahrscheinlichkeit und Generatoren
0 0.5 1 1.5 2 2.5 30
0.2
0.4
0.6
0.8
1F(x)
x
Abbildung 2.6.: Zu generierende Wahrscheinlichkeitsdichte: f(x)=sin(x)
0 0.5 1 1.5 2 2.5 30
0.2
0.4
0.6
0.8
1F(x)
x
Abbildung 2.7.: Normierte Verteilungsfunktion fur die Funktion sin(x) zwischen 0 undπ. Invertierung: die Ordinate (gleichverteilte Zufallszahl [0,1]) wird miteinem Generator erzeugt, die Abzisse ermittelt.
42
2.5. Erzeugung von gleich- und beliebig verteilten Zufallszahlen auf dem Computer
Fragen zur Selbstkontrolle
• Wie kann Wahrscheinlichkeit definiert werden?
• Welche Regeln gelten bei der Kombination von Wahrscheinlichkeiten?
• Erlautere die Begriffe Verteilungsfunktion und Wahrscheinlichkeitsdichte.
• Mit welchen einfachen Mitteln konnen Dichtefunktionen beschrieben werden?
• Wie unterscheiden sich die verschiedenen Definitionen fur den Mittelwert und ahn-liche Maße?
• Welche Eigenschaften hat die Gleichverteilung?
• Wie konnen gleichverteilte Zufallszahlen erzeugt werden?
• Wie konnen die Generatoren getestet werden?
• Wie konnen Zufallszahlen erzeugt werden, wenn die Dichtefunktion einer speziellenVerteilung folgen soll?
43
3. Spezielle Wahrscheinlichkeitsdichten
In diesem Kapitel werden die Eigenschaften von ausgewahlten speziellen Wahrschein-lichkeitsdichten in Zusammenhang mit Anwendungsbeispielen aus dem physikalischenUmfeld diskutiert.
3.1. Gleichverteilung
Die Gleichverteilung und ihre numerische Erzeugung sind im vorangegangenen Kapiteleingehend beschrieben worden. In der Physik tritt sie auf, wenn in dem beschriebenenProzess keinerlei Vorzugsrichtung vorliegt bzw. keine Selektion stattgefunden hat. Bei-spiele waren die Richtung der Emissionsachse beim Zweikorperzerfall im Ruhesystemdes Mutterteilchens, die Richtungen von Atomen oder Molekulen in einem Gas oder dieRichtungsverteilung geladener kosmischer Teilchen fern von lokalen Magnetfeldern.
3.2. Die Binomialverteilung
Die Binomialverteilung beschreibt multiple Ausfuhrungen von Versuchen, die diskretzwei Ergebnisse liefern konnen. Als Beispiele seien die Untersuchung der Ansprechwahr-scheinlichkeit einer Nachweiskammer, die Frage, mit welcher Wahrscheinlichkeit einebestimmte Ereigniszahl in einem bestimmten Bin eines Histogramms auftritt oder dieFrage, wieviele Teilchenemissionen bei einer Wechselwirkung in einen bestimmten Win-kelbereich erfolgen, genannt.
Solch ein Versuch hat die einfachsten Moglichkeit E = A + A, mit der Wahrschein-lichkeiten P (A) = p, und P (A) = 1− p = q.
Der Versuch werde n mal ausgefuhrt. Gesucht wird die Wahrscheinlichkeit fur k malA bei n Versuchen. Es folgt
P (k) = Cnk pkqn−q, Cnk =
n!
k!(n− k!).
Der Mittelwert der Binomialverteilung ist:
〈k〉 =n∑r=0
kP (k) = n · p
Die Varianz betragt:
45
3. Spezielle Wahrscheinlichkeitsdichten
V [k] =n∑r=0
(k − 〈k〉)2 · P (k) = n · p · (1− p)
Haufig soll die Effizienz eines Detektortyps bestimmt werden, der von n (z.B.=100)Teilchendurchgangen k (z.B.=80) korrekt anzeigt. Offenbar ist die Effizienz E dann E =k/n(= 80%). Der relative Fehler dieser Messung betragt ∆E =
√(k/n · (n− k)/n) =
1/n ·√
(k · (n− k))(= 4%).
Weiter ist die Schiefe skewness gegeben durch:
γ = ((1− p)− p)/√n · (1− p)p
und die curtosis durch:
C = (1− 6 · (1− p) · p)/(n · (1− p)p) + 3
Den Ubergang der Binomialverteilung in die Gaußverteilung fur große Versuchszah-len n und nicht verschwindende Wahrscheinlichkeiten p kann man zeigen, indem mandie Binomialverteilung um ihren Maximalwert, der etwa ihrem Mittelwert entspricht,entwickelt:
〈r〉 = n · p
Man betrachtet zunachst den Logarithmus von P(r):
lnP (r) = ln(pr · qn−r · n!
r!(n−r!)
)= r ln p+ (n− r) ln q + lnn!− ln r!− ln(n− r)!.
Fur den Ubergang zu großen n beuntzt man die Stirlingssche Formel
n! ≈(n+
1
2
)lnn− n+ ln 2
√(π).
Bildet man nach Einsetzen der Stirlingschen Formel die Ableitung nach r, ergibt sich:
d lnP (r)
dr= ln p− ln q − ln r + ln(n− r) = 0,
weil d lnx!dx ≈ lnx ist, flogt rmax = n · p = 〈r〉.
Nun kann die Taylorentwicklung um 〈r〉 durchgefuhrt werden:
lnP (r) = lnP (〈r〉) + (r − 〈r〉)d lnP (r)dr |r=〈r〉
+ 12!(r − 〈r〉)
2 d2 lnP (r)dr2
|r=〈r〉
46
3.3. Die Normal- oder Gauß-Verteilung
Der lineare Term ist gleich 0, da im Maximum die erste Ableitung 0 ist.
Die zweite Ableitung ist − 1npq = − 1
np(1−p) ,
So erhalten wir mit
p = k · e−r−〈r〉
2np(1−p)
die Form einer Gauß-Verteilung, mit µ = 〈r〉 und σ2 = np(1− p).
• Diese Naherung ist gut, wenn µ ≥ 10;
• Diese Naherung schlecht, wenn man sich weit außerhalb des Minimums befindet.
• Entsprechende Uberlegungen fur große Fallzahlen n aber kleine Wahrscheinlichkei-ten p fuhren auf die Poisson-Verteilung.
3.3. Die Normal- oder Gauß-Verteilung
Mit der Normal- oder Gauß-Verteilung werden sowohl zufallige Messfehler als auch Fer-tigungsfehler beschrieben. Weiter ergibt sich nach dem (spater zu diskutierenden) zen-tralen Grenzwertsatz, dass eine Uberlagerung einer großen Zahl kleinerer Einzelfehlerebenfalls durch die Gauß-Verteilung beschrieben wird. Die Wahrscheinlichkeitsdichte hatdie Form:
f(x) =1√2πσ
e−(x−µ)2
2σ2 , x ∈ [−∞,∞],
mit dem Mittelwert µ = E(x) und der Standardabweichung σ =√V (x).
Die Schiefe betragt γ = 0 und die curtosis C = 3.Die Wahrscheinlichkeitsdichte fur die Gauß-Verteilung wird auch geschrieben alsN(µ, σ2).
Ist die Abzisse so normiert, dass die Standardabweichung gleich eins ist, spricht man voneiner standardisierten Gauß-Verteilung: N(0, 1).
FIG6
47
3. Spezielle Wahrscheinlichkeitsdichten
Wegen der großen praktischen Bedeutung der Gauß-Verteilung sei zunachst auf ei-ne einfache Methode hingewiesen, die Standardabweichung aus der Darstellung einergemessenen Verteilung mit Auge und Lineal zu bestimmen. Dazu betrachtet man dieVolle Breite auf halber Hohe (FWHM). Man erhalt folgenden Zusammenhang zwischenFWHM und der Standardabweichung einer Gauß-Verteilung:
FWHM = 2σ√
2 ln 2 = 2.355σ.
Da die Form der Wahrscheinlichkeitsdichte bekannt ist, kann berechnet werden wie-viele Messungen einer großeren Messreihe inner- bzw. außerhalb von n Standardabwei-
chungen liegen. Dazu berechnet man das Integralb∫af(x)dx, das die Wahrscheinlichkeit
fur das Intervall [a = −n · σ, b = n · σ]. Man findet fur die wichtigsten Intervalle:
|x− µ| ≥ σ 31.74% |x− µ| ≤ σ 68.26%|x− µ| ≥ 2σ 4.55% |x− µ| ≤ 2σ 95.45%|x− µ| ≥ 3σ 0.27% |x− µ| ≤ 3σ 99.73%
Fuhrt man unabhangige Messungen durch und berucksichtigt man nur den statisti-schen Messfehler, erwartet man also, dass etwa ∼ 1
3 aller Messpunkte außerhalb desdurch 1 σ gegebenen Fehlerbereiches liegen muss. Um dies optisch deutlich zu machen,kann der systematische Fehler in graphischen Darstellungen separat gekennzeichnet wer-den.Liegt eine Messung außerhalb des 1 σ-Bereiches, z.B. unterhalb einer Vorhersage, kanndiese Vorhersage mit einer Wahrscheinlichkeit von (31.74/2)%+68.26% = 84.13% zuruck-gewiesen werden. (Von signifikanten Abweichungen von Vorhersagen wird erst bei Mes-sungen ausserhalb des 3 bis 5-σ-Bereiches gesprochen.)
Das Integral uber die Gaußfuntkion
φ(x) =1√2πσ
x∫−∞
e−(t−µ)2
2σ2 dt
lasst sich leider nicht in geschlossener Form berechnen. Eine Gaußverteilung kann dahernumerisch nicht durch das in Kapitel 2.5.4 vorgestellte Verfahren generiert werden. Aufdem Rechner lasst sich das Ergebnis der Integration durch Benutzung der Gauß schenFehlerfunktion bestimmen:
φ(x) =1
2
(1 + erf
(x− µ√
2σ
)).
48
3.4. Die Poisson-Verteilung
Die Gauß sche Fehlerfunktion erf(x) ist gegeben durch
erf(x) =2√π
x∫0
e−t2dt
miterf(0) = 0, erf(∞) = 1.
Sie ist in der Regel auf dem Rechner verfugbar.
3.4. Die Poisson-Verteilung
Die Poisson-Verteilung gibt die Wahrscheinlichkeit genau r Ereignisse zu erhalten,wenn die Zahl der Versuche n sehr groß ist und die Wahrscheinlichkeit fur das Auftre-ten von p in einem einzigen Versuch sehr klein ist. Die Poisson-Verteilung wird dannangewendet, wenn die Wahrscheinlichkeit p und somit der Mittelwert µ klein sind, dieVersuchszahlen n aber hoch sind.
Poisson-Verteilung: P (r) = µre−µ
r! , r = 0, 1, 2, ...
Mittelwert: 〈r〉 = µ = n · p
Ist die Wahrscheinlichkeit dafur bekannt, dass sich nichts ereignet (P (0)), ist
P (r) = P (0) = e−µ,
dann kann rekursiv weiter gerechnet werden:
P (r + 1) = P (r) · µ
r + 1.
Norm:∞∑r=0
P (r) = e−µ∞∑r=0
µr
r!= e−µ · eµ = 1.
Mittelwert:
E(r) =∞∑r=0
rµre−µ
r!=
∞∑r=0
re−µµr
r!= µ
∞∑r=0
re−µµr−1
(r − 1)!r= µe−µ
∞∑s=0
µs
s!= µ.
Varianz:V (r) = σ2 = E[(r − µ)2] = µ.
Beispiele fur die Anwendung der Poisson-verteilung sind:
49
3. Spezielle Wahrscheinlichkeitsdichten
• Die Anzahl von Teilchen, die von einem Zahler in einem Zeitintervall nachgewiesenwerden, wenn die Effizienz des Zahlers und der Fluss Φ zeitlich konstant sind unddas Produkt von Totzeit τ und Φ << 1 ist.
• Die Anzahl von Wechselwirkungen in einem dunnen Target, wenn ein intensiverTeilchenpuls einfallt.
• Der Zerfall eines kleinen Anteils von Materie in einem Zeitintervall, das einen nichtvernachlassigbaren Anteil der Lebensdauer der Quelle umfasst (Protonzefalle ineinem Jahr und 106t Material.
3.5. Die Gamma-Verteilung
Betrachten wir einen Poisson-verteilten Prozess mit dem Mittelwert µ = 1. Dann gibtdie Gamma-Verteilung die Verteilung der Wartezeiten vom ersten bis zum k−ten Ereig-nis an.
Sei k ganzzahlig und k ≥ 1, dann ist:
f(x|k) =xk−1e−x
Γ(k).
Mit (fur k ganzzahlig):Γ(k) = (k − 1)!
Eine Verallgemeinerung fur andere Werte von µ ergibt:
f(x|k, µ) =xk−1µke−µx
Γ(k).
Offenbar wird die Form der Verteilung durch k verandert, wahrend µ lediglich die Skalaverandert. Setzt man µ = 1/2, k = n/2, erhalt man eine χ2-Verteilung mit n Freiheits-graden.
Als Beispiel betrachten wir die Zeitdifferenz t zwischen zwei Ereignissen, die zufallig,mit der Rate λ aufeinander folgen. Gesucht sei die beschreibenden Wahrscheinlichkeits-dichte f(t).
Fur k − 1 Ereignisse erhalt man:
fk(t) =tk−1λke−λt
(k − 1)!.
Es kann sich dabei sowohl um die Zeitabstande (vom ersten bis k-ten Ereignis)nzwischen Wechselwirkungen in einem Detektor oder bei einem radioaktiven Zerfall alsauch zwischen dem Ausfall von elektronischen Bauelementen in einer Schaltung handeln.
50
3.6. Die χ2-Verteilung
3.6. Die χ2-Verteilung
Die Variablen x1, x2, ..., xn seien unabhangige Zufallsvariable, die alle einer Gauß-Verteilung folgen. Der Mittelwert der Verteilung sei 0, die Varianz sei 1.
Dann folgt die Summe der Quadrate
u = χ2 =
n∑i=1
x2i
einer χ2-Verteilung mit n Freiheitsgraden:
fn(u) = fn(χ2) =12
(u2
)n2−1e−
u2
Γ(n2 ).
Deren Maximum liegt bei (n− 2),
der Mittelwert bei 〈u〉 = 〈χ2〉 = n,
und die Varianz bei V [u] = V [χ2] = 2n.
Das Integral uber die Verteilung
F (u) =
∫ u
0fn(v)dv
fuhrt auf die Große
1− P = 1− F (χ2) = 1−∫ χ2
0fn(v)dv.
die die Wahrscheinlichkeit angibt, dass bei einer Stichprobe von n Messungen ein Wertvon mindestens χ2 auftritt.
Die χ2-Verteilung wird genutzt, um zu untersuchen, ob eine Messreihe durch einenhypothetisch angenommenen Zusammenhang (bezuglich dessen die Messpunkte Gauß-Verteilt sein sollten) beschrieben wird.
3.7. Die Cauchy-Verteilung
Die Wahrscheinlichkeitsdichte der Cauchy-Verteilung ist gegeben durch:
f(x) =1
π
1
1 + x2
Zu großen Werten von x fallt die Verteilung nur langsam ab, was zu einem etwasunfreundlichen Verhalten der Verteilung fuhrt. Insbesondere ist fur diese Verteilung der
51
3. Spezielle Wahrscheinlichkeitsdichten
Erwartungswert von x: nicht definiert. Alle Momente (Varianz, Schiefe, curtosis) sinddivergent (wenn zu ihrer Bestimmung bisher immer erfolgt in Grenzen bis ±∞ inte-griert wird). In der Praxis erzwingt man gelegentlich durch einen Abbruch bei endlichenArgumenten dass die Integrale nicht mehr divergieren, hat dann aber keine allgemeinvergleichbare Große mehr.
In physikalische Anwendungen taucht die Cauchy-Verteilung in der speziellen Formder Breit-Wigner-Verteilung auf:
f(x) =1
π
Γ/2
(x− x0)2 + Γ2/4.
und beschreibt die Energieverteilung von Zustanden nahe einer Resonanz, die exponenti-ell mit der Zeit zerfallen. Die Breit-Wigner-Verteilung ist symmetrisch um das Maximumbei x0. Das Gamma ist durch Γ=FWHM gegeben.
Die Cauchy-Verteilung ist mathematisch die Fouriertransformierte der Exponential-verteilung.
3.8. Die t-Verteilung
Die t-Verteilung (oder auch Studentsche-t-Verteilung) erlaubt den Test der Vertraglich-keit eines Stichproben-Mittelwertes x mit einem Erwartungswert µ oder der Vertraglich-keit der Mittelwerte zweier Stichproben.
Die Wahrscheinlichkeitsdichte der t-Verteilung ist
fn(t) =1√nπ
Γ((n+ 1)/2)
Γ(n/2)
(1 +
t2
n
)−n+12
.
Fur n = 1 ist die t-Verteilung identisch mit der Cauchy-Verteilung, und fur n→∞ gehtsie in die Gauß-Verteilung uber.
Die Herleitung dieser Verteilung basiert auf folgendem Gedankengang: Wir betrachtenzunachst eine Grundgesamtheit, die der standardisierten Normalverteilung folgt.
x sei darin das arithmetischer Mittel einer Stichprobe mit dem Umfang N (N seigroß), dann ist
σ2(x) =σ2(x)
N,
x ist normalverteilt (→ zentraler Grenzwertsatz) mit den Mittelwert x und der Streuungσ2(x). Also wird
y =x− xσx
,
52
3.8. Die t-Verteilung
durch eine standardisierte Normalverteilung beschrieben.
Statt σ(x) kennen wir jedoch nur die Schatzung
S2x =
1
N − 1
N∑j=1
(xj − x)2, fur σx
Daraus ergibt sich als Schatzung fur σ2(x):
σ2(x) =σ2(x)
N
Damit wird die Schatzung der Varianz bezuglich des Mittelwertes:
S2x =
1
N(N − 1)
N∑j=1
(xj − x)2.
Nun stellt sich die Frage, wie y = x−xσx von einer standardisierte Gauß-Verteilung ab-
weicht. Dazu wird eine Koordinatenverschiebung so duchgefuhrt, dass: x = 0.
Dann betrachteten die Verteilung der
t =x
Sx=
√Nx
Sx.
Der Ausdruck (N−1)S2x = nS2
x folgt einer χ2-Verteilung mit n = N−1 Freiheitsgraden.Einer aufwendigen Rechnung fuhrt auf die Wahrscheinlickeitsdichte
fn(t)
und die Verteilungsfunktion:
F (t) =1√nπ
Γ((n+ 1)/2)
Γ(n/2)
t∫−∞
(1 +
t2
n
)−n+12
dt.
53
3. Spezielle Wahrscheinlichkeitsdichten
3.9. Die F-Verteilung
Zwei Grundgesamtheiten haben den gleichen Erwartungswert (z.B. Messung der gleichenGroßen mit zwei verschiedenen Messgeraten, die keinen systematischen Fehler habenoder postiver Ausfall des Testes der t-Verteilung).
Mit Hilfe der F-Verteilung (Fisher-Verteilung) kann die Frage beantwortet werden, obsie auch die gleiche Streuung haben. Dazu betrachte man die normalverteilten Grund-gesamtheiten. (N1, S
21) und (N2, S
22).
Mit
S2 =1
n− 1
n∑i=1
(x− x)2,
definiert man
F =S2
1
S22
.
Sind die Streuungen gleich, wird F nahe an eins liegen.
Nach Konvention steht der großere Wert (S21 , S
22) im Zahler der Gleichung F =
S21
S22.
Daher ist F ist immer großer als eins.Die Wahrscheinlichkeitsdichte von F ist bei (n1, n2) Freiheitsgraden gegeben durch:
f(F ) =
(n1
n2
)n12 Γ((n1 + n2)/2)
Γ(n1/2) · Γ(n2/2)· F
n1−22
(1 +
n1
n2F
)−n1+n22
.
Fur wachsende und sehr große n1, und n2 konvergiert die F-Verteilung langsam gegeneine Gauß-Verteilung
Die Quantile der F-Verteilung berechnen sich nachQ∫1
f(F )dF = α, wobei α die Kon-
fidenz, d.h. die Wahrscheinlichkeit, einen Wert kleiner als Q zu erhalten, angibt
54
4. Beliebige verteilte Zufallszahlen (Teil 2)
In diesem Kapitel soll noch einmal die Erzeugung beliebig verteilter Zufallszahlen be-trachtet werden. Besondere Aufmerksamkeit wollen wir dabei der Simulation der wich-tigsten der oben vorgestellten speziellen Verteilungen widmen.
4.1. Transformation der Gleichverteilung
Dieses elegante Verfahren wurde bereits in Kapitel 2.5.4 vorgestellt. Wenn es anwendbarist, wird ohne Ereignisse verwerfen zu mussen - also ohne Verlust an Rechenzeit - furjede generierte Zufallszahl ein verwertbarer Eintrag zu der gesuchten Verteilung erzeugt.
Die Bedingungen fur die Anwendbarkeit dieser Methode sind:
• Die Verteilungsfunktion ist (x) = G(y) bekannt bzw. berechenbar. Die zugehorigeWahrscheinlichkeitsdichte muss daher integrierbar sein.
• Die Umkehrfunktion y = G−1(x) existiert.
Hat man es mit gemessenen Funktion zu tun, konnen diese durch Division durch dieNorm in Wahrscheinlichkeitsdichten transferiert werden. An diese Wahrscheinlichkeits-dichten konnen entweder geeignete Funktionen angepasst werden oder sie konnen durchHistogramme dargestellt werden und numerisch (wie in 2.5.4 diskutiert) weiterbehandeltwerden.
4.2. Das Neumannsche Ruckweisungsverfahren
Das Neumannsche Ruckweisungsverfahren kann unter schwacheren Anforderungen alsdas Transformationsverfahren angewandt werden. Hier muss nur die Wahrscheinlich-keitsdichte g(y) bekannt sein.
Gesucht wird nun eine Zufallszahl die im Bereich
a ≤ x ≤ b
entsprechend der Wahrscheinlichkeitsdichte g(y) verteilt ist.
Dazu betrachten wir in diesem Intervall
u = g(y)
55
4. Beliebige verteilte Zufallszahlen (Teil 2)
und suche eine Konstante, fur die gilt
d ≥ max{g(y)|y ∈ [a, b]}.
Nun werde Paare (yi, ui) von Zufallszahlen generiert, die Punktepaaren in den (y, u)Ebene entsprechen; yi und ui seien jeweils gleichverteilt.
Alle Paare, fur die gilt
ui ≥ g(yi)
werden verworfen. Nur die Punkten unterhalb der Kurve verbleiben.
Da die Punktedichte unterhalb der Kurve gleich ist, ist im Intervall ∆y um y die An-zahl der Punkte proportional zu g(y).
• Die Kurve g(y) muss nicht normiert sein (!).
• Das Verfahren ist ineffizient, wenn die Flache zwischen der Kurve g(y) und d = ugroß wird gegen die Flache zwischen u = 0 und der Kurve g(y).
Die Effizienz dieses Verfahrens ist:
E =
b∫ag(y)dy
(b− a)d.
Ist g(y) normiert (b∫ag(y)dy = 1), gilt
E =1
(b− a)d.
Das Effizienzproblem kann zwar nicht grundlegend behoben, wohl aber durch folgendeVerallgemeinerung abgemildert werden:
• Suche (statt der Konstanten d) eine Funktion s(y) die nahe bei g(y) liegt, g(y) je-doch nach oben beschrankt. (Beispiel: Schwarzkorperspektrum→Absorbtionsspektrumeines Stern). Dann ist:
g(y) ≤ c · s(y), a < y < b mit c < 1
• Erzeuge die beiden Zufallszahlen y gleichverteilt in a < y < b, und u gleichverteiltin 0 < y < 1.
56
4.3. Erzeugung normalverteilter Zufallszahlen
• Verwerfe y, falls
u ≥ g(y)
s(y).
• Die Effizienz ist in diesem Fall:
E =
b∫ag(y)dy
c ·b∫as(y)dy
.
• Istb∫ag(y)dy ≈
b∫as(y)dy, gilt
E =1
c.
4.3. Erzeugung normalverteilter Zufallszahlen
Die standardisierte Normalverteilung
f(x) =1√2πe−
x2
2
kann nicht nach dem Transformationsverfahren simuliert werden, da F (x)nur nume-risch mit Hilfe der Gaußschen Fehlerfunktion (erf(x)) bestimmt werden kann. Manbenutzt daher die Polarmethode, die auf einer Kombination von Transformation undRuckweisung beruht. Folgende Schritte werden dazu ausgefuhrt:
1. Erzeugung der gleichverteilten Zufallszahlen u1, u2 jeweils aus dem Intervall [0,1]. Anschließend Umformung durch die Transformationen v1 = 2u1 − 1 und v2 =2u2 − 1, sodass die Zahlenpaare v1, v2 ein Quadrat um den Ursprung mit denKantenlange 2v1, 2v2 gleichformig uberdecken.
2. Berechne s = v21+v2
2, und verwerfe (v1, v2) falls s ≥ 1 (In diesem Fall gehe zuruck zu(1)). Die verbleibenden Zahlenpaare (v1, v2) sind nun gleichformig uber die Flachedes Einheitskreises verteilt. Bilde die Polarkoordinaten des Punktes (v1, v2)
v1 = r cos θ, v2 = r sin θ
mit der Transformation
r =√s, θ = arctan (v2/v1)
.
57
4. Beliebige verteilte Zufallszahlen (Teil 2)
3. Bilde nun
x1 = v1
√−2
sln s,
x2 = v2
√−2
sln s,
x1 und x2 sind jetzt unabhangige Zufallszahlen, die der standardisierte Normal-verteilung folgen.
Zur Begrundung betrachten wir die Polarkoordinaten des Punktes (x1, x2)
x1 = cos θ√−2 ln s
x2 = sin θ√−2 ln s
Berchnet werden soll nun die Wahrscheinlichkeit, dass
√−2 ln s ≤ r =
√s
Die zugehorige Verteilungsfunktion ist:
F (r) = P (√−2 ln s ≤ r) = P (−2 ln s ≤ r2) = P (s ≥ e−
r2
2 ).
s = r2 ist gleichverteilt zwischen 0 und 1. Daher ist F (r) = 1− e−r2
2 .
Die zugehorige Wahrscheinlichkeitsdichte ergibt sich durch Differenzieren
f(r) =dF (r)
dr= re−
r2
2 .
Die gemeinsame Verteilungsfunktion von x1 und x2 ist:
F (x1, x2) = P (x1 ≤ k1, x2 ≤ k2)= P (r cos θ ≤ k1, r sin θ ≤ k2)
= 12π
∫x1<k1
∫x2<k2
re−r2
2 drdϕ
= 12π
∫x1<k1
∫x2<k2
e−x21+x
22
2 dxdy
=
(1√2π
k1∫−∞
e−x212 dx1
)·
(1√2π
k2∫−∞
e−x222 dx2
)
Dies ist das Produkt von 2 Verteilungsfunktionen der standardisierten Normalverteilung.
58
4.4. Erzeugung Poisson-verteilter Zufallszahlen
4.4. Erzeugung Poisson-verteilter Zufallszahlen
Wir betrachten eine Poisson-Verteilung mit dem Mittelwert µ:
P (r) =µre−µ
r!
Hier werden zwei Verfahren besprochen:
1. Vorgehen:
• Erzeuge exponentiell verteilte Zufallszahlen
• Summiere diese auf, bis die Summe großer ist als µ
• Dann ist die gesuchte Zufallszahl um eins kleiner als die Zahl der Summen-glieder.
Numerisch ist es gunstiger, statt exponentiell verteilten Zahlen zu addieren, derenLogarithmen, namlich gleichverteilte Zahlen zu multiplizieren und das produkt mite−µ zu vergleichen.
Exponentiell verteilte Zahlen sind gegeben durch
1
τe−
tτ
mit
t = −τ lnx.
Dann ist die Summe ∑ti = −τ
∑xi ⇒
∑tiτ
=∑
lnxi
Exponieren ergibt: e∑tiτ =
∏xi.
Zur Erlauterung betrachten wir folgenden Algorithmus. Gegeben seien:
µ, k = k,A0 = 1, n = 0
Man bildet zunachst x = e−µ, und erhoht n = n+ 1.
(*) un sei eine neue Zufallszahl. (Alle un seien im Intervall [0, 1] gleichverteilt.)Man berechnet
An = An−1 · un, n > 1
und vergleicht An mit x. Wenn An > x ist, gehe zu (*), sonst ist n die gesuchteZufallszahl.
59
4. Beliebige verteilte Zufallszahlen (Teil 2)
2. Fur große µ: kann die Gauß-Verteilung zur Naherung genutzt werden. Groß indiesem Sinn sind
µ > 10
.
Sei Z eine normalverteilte (standardisierte) Zufallszahl, dann ist
n = max(0, int(µ+ Z√µ+ 0.5)).
4.5. Erzeugung χ2-verteilter Zufallszahlen
Bei der Erzeugung χ2-verteilter Zufallszahlen unterscheidet man in Abhangigkeit vonder Zahl der Freiheitsgrade n drei Falle:
• Sei n gerade.
Bilde das Produkt von n/2 gleichformige verteilte Zahlen ui ∈ [0, 1]
x = −2 ln
n/2∏j=1
ui
.
Dann sind die x χ2-verteilte Zufallsvariablen.
• Sei n ungerade.
addiere in diesem Fall zu dem Produkt das Quadrat einer normalverteilten Zufalls-zahl Z
x = −2 ln
n/2∏j=1
ui
+ Z2.
Dann sind die x χ2-verteilte Zufallsvariablen.
• Fur große n (n >30)
erfolgt eine Naherung durch die Gauß-Verteilung. Die Funktion der χ2-verteiltenZufallsvariablen x
y =√
2χ2 − 2√
2n− 1,
ist N(0, 1)−-normalverteilt. Daher
– erzeugt man normalverteilte Zufallszahl Z
– invertiert x = 12(Z +
√2n− 1)2
– und verwirft alle Z <√
2n− 1.
60
5. Mehrdimensionale Verteilungen
5.1. Problemstellung
Betrachten wir zunachst zwei Zufallsvariablen xi und yi und fragen analog zu demeindimensionalen Fall nach der Wahrscheinlichkeit fur:
P ((X < x) ∧ (Y < y))
.
Die Verteilungsfunktion ist:
F (x, y) = P (X < x, Y < y).
Ist F differenzierbar, dann ist die Wahrscheinlichkeitsdichte
f(x, y) =∂
∂x
∂
∂yF (x, y),
und
P (a ≤ x < b, c ≤ y < d) =
b∫a
d∫c
f(x, y)dxdy.
Sei die Abhangigkeit von einer der Variablen irrelevant (x und y unabhangig), dannist
P (a ≤ x < b,−∞ ≤ y <∞) =
b∫a
∞∫−∞
f(x, y)dxdy =
b∫a
g(x)dx.
Die Wahrscheinlichkeitsdichten
g(x) =
∞∫−∞
f(x, y)dy
h(y) =
∞∫−∞
f(x, y)dx
sind die Randverteilungen der Variablen x und y.
61
5. Mehrdimensionale Verteilungen
Die Zufallsvariablen y und y sind unabhangig, wenn
f(x, y) = g(x) · h(y).
Die bedingte Wahrscheinlichkeit fur y, wenn x bekannt ist, ist gegeben durch dasIntervall
P (y ≤ Y ≤ y + dy|x ≤ X ≤ x+ dx)
Die entsprechende Wahrscheinlichkeitsdichte ist:
f(y|x) =f(x, y)
g(x),
mit der Wahrscheinlichkeit f(y|x)dy. Die Randverteilung von y bei gegebenem xist dann:
h(y) =
∞∫−∞
f(x, y)dx =
∞∫−∞
f(y|x)g(x)dx.
Wenn x und y unabhangig sind, gilt
f(y|x) =f(x, y)
g(x)=g(x) · h(y)
g(x)= h(y).
5.2. Erwartungswert, Varianz, Kovarianz und Korrelationbei zwei Variablen
Die Definition des Erwartungswertes der Funktion H(x, y) lautet jetzt
E[H(x, y)] =
∞∫−∞
∞∫−∞
H(x, y)f(x, y)dxdy.
und die Varianz von H(x, y) wird zu
σ2[H(x, y)] = E{[H(x, y)− E{H(x, y)}]2}.
Als Beispiel betrachten wir die Funktion H(x,y)=ux+by, fur die folgt:
E(ux+ by) = uE(x) + bE(y).
62
5.2. Erwartungswert, Varianz, Kovarianz und Korrelation bei zwei Variablen
Nun bilden wir mit H(x, y) = xlym, das lm-te Moment von x, y um den Ursprung:
λlm = E(xlym).
Momente um andere Punkte als den Ursprung werden dorthin verschoben:
H(x, y) = (x− a)l(y − b)m.
Der Erwartungswert des lm-ten Moments von x, y um den Punkt a, b ist dann:
µlm = E[(x− a)l(y − b)m].
Betrachtet man die Erwartungswerte in x und y
λ10 = E(x) = xλ01 = E(y) = y
Dann sind die Momente um λ10 und λ01 allgemein
µlm = E[(x− λ10)l(y − λ01)m].
Speziell ergibt sich fur die Momente
µ00 = λ00 = 1µ10 = λ10 = 0µ11 = E[(x− x) · (y − y)] = cov(x, y)µ20 = E[(x− x)2] = σ2(x)µ02 = E[(y − y)2] = σ2(y).
Wir betrachten zwei Beispiele:
1. Die Berechnung der Varianz von ax+ by
σ2(ax+ by) = E{[(ax+ by)− E(ax− by)]2}= E{[a(x− x)− b(y − y)]2}= E[a2(x− x)2 − b2(y − y)2 + 2ab(x− x)(y − y)]= a2σ2(x) + b2σ2(y) + 2ab·cov(x, y)
2. Die Berechnung H(x, y) = x, y, wenn x und y unabhangig sind.
E[xy] =
∞∫−∞
∞∫−∞
xyg(x)h(y)dxdy =
∞∫−∞
xg(x)dx ·∞∫−∞
yh(y)dy.
Es folgt, dass E(xy) = E(x) · E(y).
63
5. Mehrdimensionale Verteilungen
Die Kovarianz ist:
– positiv, wenn x > (<)x mit y > (<)y;
– negativ, wenn x > (<)x mit y < (>)y
– =0, wenn x, y unabhangig sind.
Der Korrelationskoeffizient
ρ(x, y) =cov(x, y)
σ(x) · σ(y),
ist ein grobes Maß fur die Abhangigkeit der Variablen x und y voneinander.
Rechnet man mit den reduzierten Variablen
u =x− xσ(x)
, v =y − yσ(y)
,
dann ist die Varianz der Summe:
σ2(u+ v) = σ2(u) + σ2(v) + 2ρ(u, v)σ(u)σ(v).
Nach der Definition von u und v ist
σ(u) = 1, σ(v) = 1.
Dadurch wirdσ2(u+ v) = 2 + 2ρ(u, v) = 2(1 + ρ(u, v))
desgleichen gilt fur die Varianz der Differenz
σ2(u− v) = 2 + 2ρ(u, v) = 2(1− ρ(u, v)).
Immer ist σ2 ≥ 0. Daraus folgt, dass −1 ≤ ρ(u, v) ≤ 1.
Man kann weiter zeigen, dass
ρ(u, v) = ρ(x, y)
daher ist auch−1 ≤ ρ(x, y) ≤ 1.
Betrachten wir eine Kovarianz, die exakt eins ist. aus ρ(u, v) = 1 folgt σ(u−v) = 0.Die Zufallsvariable (u− v) ist eine Konstante
u− v =x− xσ(x)
− y − yσ(y)
= const.
64
5.3. Mehrere Veranderliche
Dies ist immer erfullt, wenn
y = a+ bx, b > 0.
Bei exakt positiver Abhangigkeit zwischen x und y folgt ρ(x, y) = +1.
Bei exakt negativer Abhangigkeit gilt ρ(x, y) = −1.
Sind x und y unabhangig, dann ist die Kovarianz
cov(x, y) =∞∫−∞
∞∫−∞
(x− x)(y − y)g(x)h(y)dxdy
=∞∫−∞
(x− x)g(x)dx ·∞∫−∞
(y − y)h(y)dy
= 0.
5.3. Mehrere Veranderliche
Betrachtet man eine großere Zahl von Variablen x1, x2, ..., xn, ist die Verteilungs-funktion gegeben durch
F (x1, x2, ..., xn) = P (x1 < k1, x2 < k2, ..., xn < kn).
Ist F nach den xi differenzierbar, dann ist die gemeinsame Wahrscheinlichkeitsdichte
f(x1, x2, ..., xn) =∂n
∂x1∂x2...∂xnF (x1, x2, ..., xn).
Die Wahrscheinlichkeitsdichte der einzelnen Variablen xr, die Randverteilungen,sind
g(xr) =
∞∫−∞
...
∞∫−∞
f(x1, x2, ..., xn)dx1dx2...dxn.
︸︷︷︸n− 1–mal
Der Erwartungswert der Funktion H(x1, ..., xn) wird berechnet nach
E[H(x1, ..., xn)] =
∞∫−∞
...
∞∫−∞
H(x1, ..., xn)f(x1, ..., xn)dx1d...dxn.
65
5. Mehrdimensionale Verteilungen
Sei H(x) = xr, dann ist
E(xr) =
∞∫−∞
...
∞∫−∞
xrf(x1, x2, ..., xn)dx1dx2...dxn =
∞∫−∞
xrg(xr)dxr.
– Alle n Variablen sind unabhangig, wenn
f(x1, x2, ..., xn) = g1(x1) · g2(x2) · ... · gn(xn).
– Ein Teil von l < n Variablen sind unabhangig wenn
g(x1, x2, ..., xl) = g1(x1) · g2(x2) · ... · gl(xl).
– Die gemeinsame Randverteilung von l < n Variablen ist
g(x1, x2, ..., xl) =
∞∫−∞
...
∞∫−∞
f(x1, x2, ..., xn)dxl+1...dxn.
Die Momente der Ordnung l1, l2, ..., ln um dem Ursprung sind Erwartungswerteder Funktion
H = xl11 · xl22 · ... · x
lnn ,
d.h.
λl1,l2,...,ln = E[xl11 · xl22 · ... · x
lnn ].
Wieder ist
λ100...0 = E(x1) = x1, (5.1)
λ010...0 = E(x2) = x2, (5.2)
... = ..., (5.3)
λ000...1 = E(xn) = xn. (5.4)
Die Momente um die (x1, x2, ..., xn) sind
µl1,l2,...,ln = E[(x1 − x1)l1 · (x2 − x2)l2 · ... · (xn − xn)ln ].
Die Varianzen xi, i = 1...n sind
µ200...0 = E((x1 − x1)2) = σ2(x1), (5.5)
µ020...0 = E((x2 − x2)2) = σ2(x2), (5.6)
... = ..., (5.7)
µ000...2 = E((xn − xn)2) = σ2(xn). (5.8)
66
5.3. Mehrere Veranderliche
Die Kovarianz der Variablen xi und xj ist das Moment erster Ordnung (li = lj = 1)von den zwei Variablen (alle lk = 0) mit i 6= j 6= k
cij = cov(xi, xj) = E[(xi − xi)(xj − xj)].
Geht man uber zur vektoriellen Schreibweise (x1, x2, ...xn) →−→X , dann ist die
Verteilungsfunktion
F = F (−→X ),
und die Wahrscheinlickeitsdichte
f(−→X ) =
∂n
∂x1∂x2...∂xnF (−→X ).
Der Erwartungswert der Funktion H(−→X ) ist
E[H(−→X )] =
∫H(−→X )f(
−→X )d−→X
Varianzen und Kovarianzen konnen in einer Matrix zusammengefasst werden:
c =
c11 c12 ... c1n
c21 c22 ... c2n
... ... ... ...cn1 cn2 ... cnn
Diagonalelemente in dieser Matrix sind die Varianzen. Die Kovarianzmatrix istsymmetrisch da
cij = cji = E[(xi − xi)(xj − xj)].
Mit dem Vektor der Erwartungswerte E(−→X ) =
−→X , gilt dann
cij = E[(−→X −
−→X )(−→X −
−→X )>].
67
5. Mehrdimensionale Verteilungen
5.4. Die mehrdimensionale Gauß-Verteilung
Betrachten wir einen Vektor mit n Variablen−→X = (x1, x2, ..., xn).
Die Wahrscheinlichkeitsdichte der gemeinsamen Normalverteilung der xi ist danndefiniert als
φ(−→X ) = k · e−
12
(−→X−−→a )>B(
−→X−−→a ) = k · e−
12g(−→X ),
mit−→a : n−Komponenten Vektor,
B : n× n−Matrix, symmetrisch und positiv definit.
Es folgt, dass φ(−→X ) symmetrisch um
−→X = −→a ist, und der Erwartungswert
E(−→X −−→a ) =
∞∫−∞
...
∞∫−∞
φ(−→X )dx1...dxn = 0
Daher ist: −→E (−→X ) = −→a .
Durch das Differenzieren des Erwartungswert nach −→a erhalt man
∞∫−∞
...
∞∫−∞
[I − (−→X −−→a )(
−→X −−→a )>B]φ(
−→X )dx1...dxn = 0,
mit der Identitatsmatrix I. Der Erwartungswert in den eckigen Klammern ver-schwindet [Maximum], daher ist
E[(−→X −−→a )(
−→X −−→a )>]B = I
oder
C = E[(−→X −−→a )(
−→X −−→a )>] = B−1,
wobei C die Kovarianzmatrix ist.
Fur zwei Variablen ist
C = B−1 =
(σ2
1 cov(x1, x2)cov(x1, x2) σ2
2
)
Aus Inversion folgt
B =1
σ21σ
22 − cov(x1, x2)2
(σ2
2 cov(x1, x2)cov(x1, x2) σ2
1
)
68
5.4. Die mehrdimensionale Gauß-Verteilung
Wenn die Kovarianzen verschwinden, ist B diagonal
B = B0 =
(1σ21
0
0 1σ22
)
Einsetzen in die Wahrscheinlicheitsdichte liefert:
φ(x) = k · e−12
(−→X−−→a )>B0(
−→X−−→a ) = k · e
− 12
(x1−a1)2
σ21 e− 1
2(x2−a2)
2
σ22 ,
das ist das Produkt von zwei unabhangigen normalverteilten Variablen.Fur zwei Variablen ist die Normierung:
k = k0 =1
2πσ1σ2
und fur beliebige viele Variablen
kn =
(detB
(2π)n
)1/2
.
Seien die Variablen nicht unabhangig. Dann betrachten wir die reduzierten Varia-blen
ui =xi − aiσi
, i = 1, 2...,
und den zugehorigen Korrelationskoeffizienten
ρ =cov(x1, x2)
σ1σ2= cov(u1, u2),
dann erhalten wir
φ(u1, u2) = k · e−12−→u >B−→u = k · e−
12g(−→u ),
B =1
1− ρ2
(1 −ρ−ρ 1
)
In einer zweidimensionalen Darstellung konnen Linien gleicher Wahrscheinlich-keitsdichte als ’Hohenlinien’ gezeichnet werden. Dazu setzt man
φ(u1, u2) = const ⇒ −1
2g(−→u ) = const
⇒ −1
2
1
1− ρ2(u2
1 + u22 − 2u1u2ρ) = const.
69
5. Mehrdimensionale Verteilungen
Betrachten wir o.B.d.A. g(−→u ) = 1, und wahlen die ursprunglichen Variablen
(x1 − a1)2
σ21
− 2ρ(x1 − a1)
σ1· (x2 − a2)
σ2+
(x2 − a2)2
σ22
= 1− ρ2,
erhalten wir eine Ellipsengleichung mit den Mittelwert u1u2.Aus der Geometrie wissen wir, dass die Hauptachsen der Ellipse bezuglich der Ach-sen x1 bzw. x2 den Winkel α bilden und entlang der Hauptachsen die Halbmesserp1 und p2 aufweisen.
Die Kegelschnitte liefern
tan 2α =2ρσ1σ2
σ21 − σ2
2
p1 =σ2
1σ22(1− ρ2)
σ22 cos2 α− 2ρσ1σ2 sinα cosα+ σ2
1 sin2 α
p2 =σ2
1σ22(1− ρ2)
σ22 sin2 α− 2ρσ1σ2 sinα cosα+ σ2
1 cos2 α
als Kovarianzellipse.
Die Kovarianzellipse liegt in einen Rechteck, das durch den Punkt (u1, u2) und σ1
und σ2 bestimmt ist. Sie beruhrt das Rechteck in vier Punkten. Ist ρ = ±1, gehtdie Ellipse in eine Diagonale des Rechtecks uber.
Als 1σ-Bereich, wird analog zum eindimensionalen Fall die Region bezeichnet,innerhalb deren 68,39% der Ereignisse liegen.
70
6. Einfache statistische Methoden
6.1. Trennung von Datensatzen: Diskriminanzanalyse
– Fisher, Annals of Eugenics 7, 179-188, 1936.
Ein Datensatz bestehe aus zwei Unterklassen: A – Signal und B – Untergrund.Beide Klassen liegen als Ergebnis einer Monte-Carlo-Simulation vor. Die Methodeist auf experimentell bestimmte Datensatze nur dann anwendbar, wenn die gene-rierten und gemessenen Verteilungen ubereinstimmen.
Beide Datensatze konnen von den p ggf. korrelierten Zufallsvariablen x1, ...xp ab.Die Ereignisse in beiden Datensatzen sollen optimal voneinander getrennt werden.Ein einfacher Ansatz ware, dies durch Schnitte in den Einzelvariablen zu realisieren.In der Regel wird eine optimale Trennung aber durch eine geeignet gelegte p-1dimensionale Hyperebene erreicht.
Um diese Ebene zu berechnen definiert man zunachst die sog. Diskriminanzfunk-tion:
F =−→λ ·−→X
−→λ = (λ1, ..., λp) Parameter
−→X = (x1, ..., xp) Zufallsvariablen.
Mit diesem Werkzeug werden folgende Uberlegungen durchgefuhrt:
– Die Mittelwerte der Variablen werden berechnet
−→µ A =1
na·nA∑i=1
−→x A,i ,
−→µ B =1
nb·nB∑i=1
−→x B,i ,
nA, nB sind hier die Anzahl der Messungen.
Der Erwartungswert der Diskriminanzfunktion ist damit gegeben durch
FA =−→λ · −→µ A
71
6. Einfache statistische Methoden
FB =−→λ · −→µ B.
– Die Varianz V der F-Verteilung ist durch die Summe der Varianzen von Signalund Untergrund gegeben:
V =
nA∑k=1
(Fk − FA)2 +
nA+nB∑k=nA+1
(Fk − FB)2,
mit n = nA + nB gilt dann
V =
n∑k=1
λiλjSij .
Die Kombinierte Kovarianzmatrix Sij ist
Sij =(nA − 1) · SAij + (nB − 1) · SBij
nA + nB − 2,
wobei die SA,Bij =die Kovarianzmatritzen der einzelnen Klassen sind:
SA,Bij =1
nA,B − 1
nA,B∑k=1
(xi,k − xi) · (xj,k − xj)
es folgt
V =−→λ · S ·
−→λ .
– Mit der Distanzfunktion D2 wird der Abstand zwischen den Mittelwerts derKlassen A und B bestimmt:
D2 = (FA − FB)2 = [−→λ (−→µ A −−→µ B)]2
– Die Trennung zwischen den Klassen ist umso besser, je großer das Verhaltnisφ von Distanz D2 zur Varianz V ist.
Maximiere φ =D2
V=|−→λ · −→µ A −
−→λ · −→µ B|2
−→λ · S ·
−→λ
Gesucht ist somit das Maximum von φ an der Stelle, an der die erste Ableitungvon φ nach Z eine Nullstelle hat.
dφ
dλ=
2 · (−→µ A −−→µ B) ·−→λ · S ·
−→λ − 2S
−→λ (−→λ−→µ A −
−→λ−→µ B)
(−→λ · S ·
−→λ )2
.
72
6.1. Trennung von Datensatzen: Diskriminanzanalyse
dφ
dλ= 0 ⇒ −→µ A −−→µ B = S ·
−→λ
(−→λ−→µ A −
−→λ−→µ B)
−→λ · S ·
−→λ
Offenbar kann λ nur bis auf eine multiplikative Konstante, die sich in der Glei-chung fur φ herauskurzt, bestimmt werden. Eine Normierung erfolgt durch
(−→λ−→µ A −
−→λ−→µ B)
−→λ · S ·
−→λ
= 1 =D2
V.
Es folgt: −→λ = S−1(−→µ A −−→µ B).
– Variablen-Selektion: Aus Messungen stehen einem in der Regel mehr moglicheParameter zur Verfugung als sinnvoll fur die Diskriminanzanalyse einzusetzensind. Um geeignete Parameter zu selektieren, kann man fordern, dass sich dieMittelwerte der Klassen sollen mit moglichst großer Signifikanz unterscheiden.Dies kann mit einem t-Test gepruft werden:
t =xA − xB
S
√nA · nBnA + nB
mit der empirischen Streuung
S2 =1
nA + nB − 2
(nA∑i=1
(xA,i − xA)2 +
nB∑i=1
(xB,i − xB)2
)Dieser Test ist fur das Verfahren gunstig aber nicht notwendig.
– Um die Qualitat der Trennung beurteilen zu konnen werden folgende Histo-gramme gefullt und diskutiert:
1. Die Verteilungen der beiden Klassen in Abhangigkeit von λ (in einemHistogramm).
2. Die Reinheit (purity) in Abhangigkeit von λ:
Dazu plottet man in ein Histogramm fur A
RA =
k∑bin=0
nA,k
k∑bin=0
(nA,k + nB,k)
und fur B
RB =
k∑bin=kmax
nB,k
k∑bin=kmax
(nB,k + nA,k)
73
6. Einfache statistische Methoden
3. Ebenfalls in einem Histogramm wird die Effizienz bei einem Schnitt in λfur die beiden Klassen dargestellt:
εA =1
nA
k∑bin=0
nA,k
εB =1
nB
k∑bin=0
nB,k.
Mit Hilfe dieser Histogramme kann festgelegt werden, wieviele Ereignisse des Si-gnalsamples sich in dem selektierten Sample befinden und durch wieviele Unter-grundereignisse das selektierte Sample noch verunreinigt wird. Je nach untersuch-ter Fragestellungen konnen die Anforderungen unterschiedlich ausfallen.
6.2. Theoreme und Satze
6.2.1. Tschebyscheff-Ungleichung
Gesucht ist obere Schranke fur die Wahrscheinlichkeit, dass eine Zufallszahl mehrals kσ von Mittelwert abweicht. Die Varianz und der Mittelwert seien bekannt.
σ2 =
∞∫−∞
(x− 〈x〉)2 · f(x)dx =
〈x〉−kσ∫−∞
+
〈x〉+kσ∫〈x〉−kσ
+
∞∫〈x〉+kσ
(x− 〈x〉)2 · f(x)dx
Das Wegglassen des mittleren Terms fuhrt auf eine Ungleichung, da alle Teilinte-grale positiv sind:
σ2 ≥
〈x〉−kσ∫−∞
+
∞∫〈x〉+kσ
(x− 〈x〉)2 · f(x)dx
Fur das erste Integral gilt:x < 〈x〉 − kσ
⇒ x− 〈x〉 < −kσ
⇒ (x− 〈x〉)2 > k2σ2.
Fur das zweite Integral gilt:x > 〈x〉+ kσ
74
6.2. Theoreme und Satze
⇒ x− 〈x〉 > kσ
⇒ (x− 〈x〉)2 > k2σ2.
Einsetzen liefert die Ungleichung:
σ2 ≥ k2σ2
〈x〉−kσ∫−∞
f(x)dx+
∞∫〈x〉+kσ
f(x)dx.
Die Integrale geben die Wahrscheinlichkeit dafur an, dass die Zufallszahl aus demBereich|x− 〈x〉| ≥ kσ stammt. Man erhalt
〈x〉−kσ∫−∞
f(x)dx+
∞∫〈x〉+kσ
f(x)dx ≤ 1
k2.
Satz. Die Wahrscheinlichkeit dafur, dass die xi aus dem Intervall |x − 〈x〉| ≥ kσgezogen werden, ist kleiner gleich 1
k2. Dies gilt unabhangig von der Wahrschein-
lichkeitsdichte f(x).
Leider ist diese Bedingung nur schwach.
6.2.2. Gesetz der großen Zahl
Gegeben seien n unabhangige Experimente, in denen das Ereignis j genau nj malaufgetreten ist. Die nj seien binomialverteilt. Der Bruchteil hj =
njn sei die be-
trachtete Zufallsvariable.Dann ist der Erwartungswert von hj
E(hj) = E(nj/n) = pj ,
d.i. Wahrscheinlichkeit pj fur das Ereignis j.
Wie genau ist die Schatzung fur die unbekannte Wahrscheinlichkeit pj?
Da die nj = n · hj binomialverteilt sind, gilt fur die Varianz
V (hj) = σ(hj) = σ(nj/n) =1
n2· σ2(nj) =
1
n2· npj(1− pj).
Da pj(1− pj) ≤ 14 ist, gilt
σ2(hj) <1
n,
Das sog.Gesetz der großen Zahl besagt, dass der Fehler der Schatzung hj der Wahr-scheinlichkeit pj durch 1/
√n beschrankt ist.
75
6. Einfache statistische Methoden
6.2.3. Der Zentrale Grenzwertsatz
Die Wahrscheinlichkeitsdichte der Summe w =n∑i=1
xi einer Stichprobe aus n un-
abhangigen Zufallsvariablen xi mit einer beliebigen Wahrscheinlichkeitsdichte mitdem Mittelwert 〈x〉 und der Varianz σ2 geht im Grenzfall n → ∞ gegen eineGauß-Verteilung mit dem Mittelwert 〈w〉 = n · 〈x〉 und einer Varianz V (w) = nσ2.
Dieses gilt auch, wenn mehrere Wahrscheinlichkeitsdichten uberlagert werden.
Großen, die auf Summen von zufallsverteilten Ereignissen basieren, sollten Gauß-Verteilt sein.
Beweis
Betrachten wir zunachst den Mittelwert
E(w) = E
(n∑i=1
xi
)=
n∑i=1
E(xi) = nE(x) = n〈x〉
und dann die Varianz
V (w) = E[(w − 〈w〉)2]
= E
[(n∑i=1
xi −∑〈xi〉
)2]
= E
[(n∑i=1
(xi − 〈xi〉))2]
= E
[n∑i=1
(xi − 〈xi〉)2 +∑i 6=k
(xi − 〈xi〉)(xk − 〈xk〉)
]
Da die Variablen nicht korrelieren, ist∑i 6=k
(xi−〈xi〉)(xk−〈xk〉) = 0, und wir erhalten
V (w) = E
[n∑i=1
(xi − 〈xi〉)2
]= nV (x).
76
6.3. Methode der kleinsten Quadrate
6.3. Methode der kleinsten Quadrate
6.3.1. Vorbemerkungen
Legendre, Gauß, Laplace; Beginn 19. Jahrhundert
Situation: Wiederholte Messung einer Große yi(~xi). Interpretation als Summe einerwahren Große yw plus Meßfehler ξyi . Suche ein yw, sodass die Quadratsumme derFehler minimal wird: ∑
i
ξ2yi =
∑i
(yw − yi)2 = min. (6.1)
Die Abweichung zwischen den yi und den yw wird durch die Standardabweichungσ bzw. die Varianz σ2 beschrieben.yi: Stichprobe ↔ WahrscheinlichkeitsdichteAnnahme: Es existiert eine funktionelle Beziehung zwischen den yi und ~xi, das“Modell”.Dieses Modell kann von den zusatzlichen Variablen bzw. Parametern aj abhangen.Im allgemeinen Fall darstellbar durch eine oder mehrere Gleichungen/Beziehungender Form:
f(a1, a2, . . . ap, y1, . . . yn) = 0 (6.2)
Einfachster Fall:
- Daten unkorreliert
- alle haben das gleiche σ
⇒ S =
n∑i=1
∆y2i = minimal. (6.3)
Beispiel: n wiederholte Messungen yi der Variablen y.Bedingung S =
∑i (y − yi)2 = min wird erfullt durch den Mittelwert:
→ y =∑i
yin
, (6.4)
wobei y der “Schatzwert” von y ist.
Allgemeiner Fall:
- Daten beschrieben durch den Vektor y (der Dimension n)
- verschiedene Standartabweichungen σ
- Korrelationen durch die Kovarianzmatrix V gegeben
77
6. Einfache statistische Methoden
Matrixschreibweise:
S = ∆yTV −1∆y = min, (6.5)
wobei ∆y der Residuenvektor ist.
Haufige Anwendung:
Modell: y = f(a, a1, a2, . . . ap)
Meßergebnisse sollen uber eine Funktion f von den Parametern aj abhangen.
→ Minimierung von S
→ Bestimmung der Parameter
→ Suche nach funktionellen Zusammenhangen
→ Test, ob Form der Parametrisierung uberhaupt vertraglich mit den Meßdatenist
6.3.2. Kleinste Quadrate in linearen Modellen
Ziel: Bestimmung der Parameterwerte ~a = (a1 . . . ap).
yi = f(xi, ~a) (6.6)
f(x, ~a) hangt linear von den p Parametern aj ab.
y(x) = f(y, ~a) = a1f1(x) + a2f2(x) + . . .+ apfp(x) (6.7)
n Wertepaare xi, yi werden bestimmt Erwartungswert der Einzelmessung yi:
E[yi] = f(xi, ~a) = yi (6.8)
~a: wahre Werte von ~a.
Residuen ri:
ri = yi − f(xi, ~a) (6.9)
Fur ~a = ~a gilt
E[~ri] = 0 , E[r2i ] = V [ri] = σ2 (6.10)
- keine Aussage uber Wahrscheinlichkeitsdichte gemacht
- aber “unverzerrt”, Varianz endlich
- (hier) auch unkorreliert ⇒ Kovarianzen verschwinden
78
6.3. Methode der kleinsten Quadrate
∼Zu minimieren ist:
S =∑i
r2i =
∑i
[yi − a1f1(xi)− a2f2(xi)− . . .− apfp(xi)]2 (6.11)
→ alle partiellen Ableitungen nach den aj mussen verschwinden:
∂S
∂a1= 2
∑if1(xi) [a1f1(xi) + a2f2(xi) + . . .+ apfp(xi)− yi] = 0 (6.12)
∂S
∂a2= 2
∑if2(xi) [a1f1(xi) + a2f2(xi) + . . .+ apfp(xi)− yi] = 0 (6.13)
...
u.s.w.
→ Umschreiben in die sog. Normalgleichungen
a1∑if2
1 (xi) + . . .+ ap∑if1(xi)fp(xi) =
∑i
yif1(xi) (6.14)
a1∑if2(xi)f1(xi) + . . .+ ap
∑if2(xi)fp(xi) =
∑i
yif2(xi) (6.15)
...
u.s.w., p Gleichungen
→ Umschreiben in Matrixschreibweise
n× p Werte fj(xi) → Elemente in einer n× p Matrix A
A =
f1(x1) f2(x1) · · · fp(x1)f1(x2) f2(x2) · · · fp(x2)
......
. . ....
f1(xn) f2(xn) · · · fp(xn)
(6.16)
A heißt auch Design-Matrix.
~a =
a1
a2...ap
~y =
y1
y2...yp
(6.17)
79
6. Einfache statistische Methoden
A · ~a = Vektor der Erwartungswerte.
Residuen:~r = ~y −A~a (6.18)
Minimierungsbedingung
S = ~rT~r = (~y −A~a)T (~y −A~a) (6.19)
= ~yT~y − a~aTAT~y + ~aTATA~a (6.20)
Bed. ist−→ 2AT~y + 2ATA~a = 0 (6.21)
oder in Matrixform der Normalgleichungen:
(ATA)~a = AT~y (6.22)
~a = (ATA)−1AT~y (6.23)
mit der symmetrischen p× p-Matrix
ATA =
∑
i f21 (xi)
∑i f1(xi)f2(xi) · · ·
∑i f1(xi)fp(xi)∑
i f2(xi)f1(xi)∑
i f22 (xi) · · ·
∑i f2(xi)fp(xi)
......
. . ....∑
i fp(xi)f1(xi)∑
i fp(xi)f2(xi) · · ·∑
i f2p
(6.24)
und
AT~y =
∑
i yif1∑i yif1...∑i yifp
(6.25)
Die Kovarianzmatrix
Betrachte: ~a = (ATA)−1AT~y.~a geht aus einer linearen Transformation von ~y hervor. Kovarianzmatrix von ~a kann uber Fehlerfortpflanzung aus der Kovarianzmatrixvon ~y berechnet werden:
V [~y] =
var(y1) cov(y1, y2) · · · cov(y1, yn)
cov(y2, y1) var(y2) · · · cov(y2, yn)...
.... . .
...cov(yn, y1) cov(yn, y2) · · · var(yn)
(6.26)
var(yi) = σ2i , cov(yi, yk) = σik (6.27)
80
6.3. Methode der kleinsten Quadrate
Hier: ~yi unkorreliert, alle σi gleich
V [~y] = σ21 =
σ2 0 · · · 0
0. . .
. . ....
.... . .
. . . 00 · · · 0 σ2
(6.28)
Sei die Beziehung linear~a = B~y (6.29)
mit Fehlerfortpflanzung gilt:
V [~a] = BV [~y]BT (6.30)
B = (ATA)−1AT (6.31)
V [~a] = (ATA)−1ATV [~y]A(ATA)−1 (6.32)
mit V [~y] = σ21 (6.33)
V [~a] = σ2(ATA)−1 (6.34)
– V [~a] hangt nicht von den Residuen ab.
– σ2 hat keinen Eifluss auf die Parameterwerte, wenn σ2i = σ2.
Quadratsumme der Residuen
~a = (ATA)−1AT~y (6.35)
einsetzen inS = ~yT~y − 2~aTAT~y + ~aTATA~a (6.36)
ergibt
S = ~yT~y − 2~aTAT~y + ~aTATA(ATA)−1AT~y (6.37)
= ~yT~y − ~aTAT~y (6.38)
– Summe der Residuenkann direkt berechnet werden.
– aber:
∗ Differenz großer Zahlen
∗ Einzelbeitrage ggf. interessant
E[S] = σ2(n− p) (6.39)
Ist die Varianz unbekannt, kann abgeschatzt werden:
σ2 =S
n− p(6.40)
Fur große Werte ist das eine gute Schatzung.
81
6. Einfache statistische Methoden
Eigenschaften der Losung dieses Problems
Zur Erinnerung
- Daten sind unverzerrt (unbiased), “erwartungstreu”
- alle Varianzen gleich
d.h.
- E[~y −A~a ] = 0 oder E[~y] = A~a
- V [~y −A~a ] = σ21
Dann gilt:
1. Die Schatzwerte, die mit der Methode der kleinsten Quadrate ermittelt wer-den sind unverzerrt.
2. Sie haben den kleinsten Fehler von allen linearen Schatzwerten d.i. ”Gauß-Markoff-Theorem”.
3. Dies gilt unabhangig von der Wahrscheinlichkeitsdichte der Residuen.
6.3.3. Gaußverteilte Meßfehler
→ Wahrscheinlichkeitsdichte von S bekannt.Sσ2 folgt einer χ2-Verteilung mit (n− p) Freiheitsgraden.
χ2-Test zur Uberprufung, ob Daten und Modell vertraglich sind.→ Dieser wird auch bei (kleinen) Abweichungen von Gauß durchgefuhrt.
Unterschiedliche Fehler
– Einzelne Datenpunkte haben verschiedene Genauigkeiten:
Varianz V [yi] = σ2i (6.41)
– Datenpunkte statistisch unabhangig cov = 0
V [~y] =
σ2
1 0 · · · 0
0 σ22
. . ....
.... . .
. . . 00 · · · 0 σ2
n
(6.42)
Residuenquadrate
S =∑i
r2i
σ2i
?= min (6.43)
82
6.3. Methode der kleinsten Quadrate
Einfuhrung der Gewichtungsmatrix
W [~y] = V −1[~y] =
1/σ2
1 0 · · · 0
0 1/σ22
. . ....
.... . .
. . . 00 · · · 0 1/σ2
n
(6.44)
→ Diagonalelemente = inverse Varianzen.Ausdruck fur S:
(bisher: S = ~rT~r) (6.45)
nun: S = ~rTW~r (6.46)
= (~y −A~a)TW [~y](~y −A~a) (6.47)
Gilt auch fur korrelierte Datenpunkte!W [~y] = V −1[~y] ist dann nicht mehr diagonal.
Uberlegung:
– V [~y] ist symmetrisch.
– lineare Algebra: Zu jeder symmetrischen Matrix V [~y] gibt es eine orthogonaleMatrix U , die V [~y] in eine Diagonalmatrix V [~z] transformiert.
– dann ist
~z = UT~y (6.48)
und nach Fehlerfortpflanzung
V [~z] = UTV [~y]U (6.49)
S = (~z − UTA~a)T W [~z]︸ ︷︷ ︸diagonal
(~z − UTA~a) (6.50)
(vergleiche Formel 6.47.)
UW [~z]UT = V −1[~z] = W [~y] (6.51)
Losung fur allgemeine Kovarianzmatrix
S = ~rTW [~y]~r (6.52)
= (~y −A~a)TW [~y](~y −A~a) (6.53)
83
6. Einfache statistische Methoden
fuhrt fur den Schatzwert fur ~a auf
~a = (ATWA)−1ATW~y (6.54)
V [~a ] = (ATWA)−1 (6.55)
vergleiche:
~a = (ATWA)−1AT~y (6.56)
V [~a ] = (ATA)−1σ2 (6.57)
Summe der Residuenquadrate fur ~a = ~a:
S = ~yTW~y − ~aTATW~y (6.58)
E[S] = n− p”freie Parameter“ (6.59)
Praktische Hinweise
Geradenanpassungy = f(x, a1, a2) = a1 + a2x (6.60)
→ Designmatrix A
A =
1 x1
1 x2...
...1 xn
(6.61)
Messwerte seien unkorreliert V [~y] und W diagonal.
ATWA =
( ∑iWi
∑iWixi∑
iWixi∑
iWix2i
)=
(S1 SxSx Sxx
)(6.62)
ATW~y =
( ∑iWiyi∑iWixiyi
)=
(SySxy
)(6.63)
→ 6 Summen berechnen
(ATWA)−1 =1
D
(Sxx −Sx−Sx S1
)(6.64)
mit DeterminanteD = S1Sxx − S2
x (6.65)
a = (ATWA)−1ATWy (6.66)
Einsetzen:
a1 = (SxxSy − SxSxy)D (6.67)
a2 = (−SxSy + S1Sxy)D (6.68)
84
6.3. Methode der kleinsten Quadrate
Kovarianzmatrix:
V [a] =1
D
(Sxx −Sx−Sx S1
)(6.69)
Summe der Residuenquadrate:
S = Syy − a1Sy − a2Sxy (6.70)
2 Sonderfalle
1. Fehler in beiden Variablen?
y = a1 + a2x (6.71)
xi mit σxi (6.72)
yi mit σyi (6.73)
Zu minimieren ist dann die Summe der Quadrate des Abstands der Fehlerel-lipsen von der Geraden
S(a1, a2) =∑i
(yi − a1 − a2x)2
σ2yi + a2
2σ2xi
(6.74)
Gefordert∂S
∂a1= 0 und
∂S
∂a2= 0 (6.75)
Numerische Methoden gefordert, z.B. Variation von a2 und Berechnung von
a1 =
∑i
yiσ2yi
+a22σ2xi
−∑
ixi
σ2yi
+a22σ2xi∑
i 1/(σ3yi + a2
2σ2xi)
(6.76)
oder gleich Anwendung von Optimierungsmethoden.
2. Lineare Regression
Minimierung der Summe der Quadrate der senkrechten Abstande der Punktevon einer Geraden. Alle Standartabweichungen gleich σ,
d.h. S =∑i
(yi − a1 − a2xi)2
(1 + a22)σ2
, (6.77)
fuhrt auf
a1 = y − a2x (6.78)
a2 = q ±√q2 + q (6.79)
y =∑i
yin
, x =∑i
xin
(6.80)
q =
∑i (yi − y)2 −
∑i (xi − x)2
2∑
i (yi − y)(xi − x)(6.81)
Vorzeichen testen!
85
6. Einfache statistische Methoden
6.3.4. Nichtlineare kleinste Quadrate
f(x,~a) hangt nicht linear von den ai ab, z.B.
f(ax,~a) = a1 · exp(a2x). (6.82)
∂f∂a1
und ∂f∂a2
hangen von den Parametern ai ab.⇒
”nichtlinear“.
Allgemeine Losung durch Iteration!
vgl. Ausweg:
Linearisierung
Taylor-Entwicklung von f(x,~a):
f(x,~a) = f(x,~a∗) +
p∑j=1
∂f
∂aj(aj − a∗j )︸ ︷︷ ︸
Ableitungen an der Stelle a∗
(6.83)
Korrekturterm:∆~a = ~a− ~a∗ (6.84)
Die Residuenri = yi − f(xi,~a) (6.85)
sind in linearer Naherung~r = ~y −A∆~a− ~f (6.86)
mit der Jacobi-Matrix
A =
∂f(x1)/∂a1 ∂f(x1)/∂a2 · · · ∂f(x1)/∂ap
∂f(x2)/∂a1. . .
. . ....
.... . .
. . ....
∂f(xn)/∂a1 · · · · · · ∂f(xn)/∂ap
(6.87)
und dem Naherungsvektor ~f an der Stelle ~a∗.
S = ~rTW~r (6.88)
= (~y −A∆~a− ~f)TW (~y −A∆~a− ~f) = min (6.89)
Normalgleichung
(ATWA)∆~a = ATW (~y − ~f) (6.90)
mit ∆~a = (ATWA)−1ATW (~y − ~f) (6.91)
86
6.4. Nachtrag und Exkurs: Fehlerfortpflanzung
Addiere dann
~a∗ + ∆~a→ ~a∗ (6.92)
in der Regel bessere Naherung.Naherungsweise weiter gultig:
a = (ATWA)−1ATWy (6.93)
V [a] = (ATWA)−1 (6.94)
E[S] = n− p (6.95)
Konvergenz
Man hofft, dass
S(~a∗ + ∆~a) < S(~a∗) (6.96)
Bestimmt ist aber:
S(~a∗ + λ∆~a) < S(~a∗) (6.97)
→ Suche mit Salamitaktik nach geeignetem λ.
Konvergenzkriterium ?
Betrachte
∆S = ∆~aTW (~y − ~f) (6.98)
fur ∆S < 1 Abweichung innerhalb 1σ.Beenden, wenn z.B. ∆S < 0, 1, weil statistische Fehler dann dominieren.Stark nichtlineares Problem? allgemeine Optimierung.
6.4. Nachtrag und Exkurs: Fehlerfortpflanzung
6.4.1. Transformation einer Variablen
– Wahrscheinlichkeitsdichte f(x) bekannt
– Trannsformation von x in die Variable y(x)
→ Gesucht: Wahrscheinlichkeitsdichte f(y)
Betrachte:in x in y
Interval (x, x+ dx) (y, y + dy)
Flachentreue Abbildung
fx(x) dx = fy(y) dy (6.99)
87
6. Einfache statistische Methoden
fy(y) = fx(x(y))
∣∣∣∣dxdy∣∣∣∣ (6.100)
Falls Transformation nicht eindeutig
fy(y) =∑
Zweige
fx(x(y))
|dy/dx|(6.101)
(Keine Ausloschung! → Vorzeichen von dy/dx darf sich nicht andern.)Bsp.:
y = x2 → x =√y (6.102)
fy(y) =1
2√y
(fx(+√y) + fx(−√y)) (6.103)
Transformation von Mittelwert und Varianz
Nichtlinearitat → komplizierte Ausdrucke→ Entwicklung bis zur 2. Ordnung
y(x) ≈ y(〈x〉) + (x− 〈x〉) dydx
∣∣∣∣x=〈x〉
+1
2(x− 〈x〉)2 d
2y
dx2
∣∣∣∣x=〈x〉2
(6.104)
E[y] = y(〈x〉) + E[x− 〈x〉] dydx
∣∣∣∣x=〈x〉︸ ︷︷ ︸
=0 nach Def.
+1
2E[(x− 〈x〉)2]︸ ︷︷ ︸
=σ2x
d2y
dx2
∣∣∣∣x=〈x〉
(6.105)
〈y〉 = x(〈x〉) +1
2σ2x
d2y
dx2
∣∣∣∣x=〈x〉
(6.106)
Bei Vernachlassigung der 2. Ordnung
〈y〉 ≈ y(〈x〉) (6.107)
Varianz:
V [y] = E[(y − 〈y〉)2] = E
[dy
dx
∣∣∣∣〈x〉· (x− 〈x〉)2
](6.108)
=
(dy
dx
)2
· E[(x− 〈x〉)2] (6.109)
=
(dy
dx
)2
· V [x] (6.110)
d.i. Fehlerfortpflanzung fur eine Variable.
88
6.4. Nachtrag und Exkurs: Fehlerfortpflanzung
Bsp.: Lineare Transformation
y = ax+ b x =y − ba
(6.111)
fy(y) =1
|a|fx
(y − ba
)(6.112)
E[y] = a · E[x] + b (6.113)
〈V 〉 = a · 〈x〉+ b (6.114)
→ V [y] = E[(ax+ b− (a 〈x〉+ b))2] (6.115)
= a2 · E[(a− 〈x〉)2] = a2 · V [x] (6.116)
6.4.2. Transformation mehrerer Variablen
Lineare Transformation
~y = B~x (6.117)
~x: n-Vektor, Mittelwert ~µx, Kovarianmatrix V [~x]~y: m-Vektor, Mittelwert ~µy, Kovarianmatrix V [~y]
B: m× n-Matrix mit Bik = ∂yi∂xk
, m 6= n. Fehlerfortplanzunggesetz
〈~y〉 = ~µy = B~µx (6.118)
V [~y] = BV [~x]BT , (6.119)
da
〈~y〉 = E[~y] = E[B~x] = BE[~x] = B~µx (6.120)
und
V [~y] = E[(~y − ~µy)(~y − ~µy)T ] (6.121)
= E[(B~x−B~µx)(B~x−B~µx)T ] (6.122)
= E[B(~x− ~µx)(~x− ~µx)TBT ] (6.123)
= BE[(~x− ~µx)(~x− ~µx)T ]BT (6.124)
= BV [x]BT (6.125)
Spezialfall n = m→ B quadratisch
fy(~y) = fx(~x(~y)) · det∣∣B−1
∣∣ , (6.126)
mit
~x(~y) = B−1~y (6.127)
Allgemein (nicht unbedingt linear):
89
6. Einfache statistische Methoden
Transformationsgleichung:
yi = yi(~x) ~x = (x1, . . . , xn) (6.128)
Wenn dim(~y) = dim(~x) ist, kann fy(~y) ggf. aus fx(~x) berechnet werden. Wiederist
fx(~x) dx1 . . . dxn = fy(~y) dy1 . . . dyn. (6.129)
fx(~x), fy(~y) sind n-dimensionale Wahrscheinlichkeitsdichten.
fy(~y) = fx(~x) · J(x
y
), (6.130)
mit der Jacobideterminanten
J(x
y
)=
∣∣∣∣∂x1∂x2 · · · ∂xn∂y1∂y2 · · · ∂yn
∣∣∣∣ (6.131)
Wenn dim(~y) 6= dim(~x) ist:Naherungsrechnung fur die Kovarianzmatrix:Betrachte die Matrix der Ableitungen an der Stelle des Mittelwertes ~µy
B =
∂y1∂x1
∂y1∂x2
. . . ∂y1∂xn
.... . .
. . ....
∂ym∂x1
. . . . . . ∂ym∂xn
(6.132)
Allgemeines Gesetz der Fehlerfortpflanzung
V [~y] = BV [~x]BT |lin. Fall: B Transfermatrix (6.133)
mit den Kovarianzmatrizen V [~y] und V [~x].
Herleitung:Entwicklung von yi(~x) um den Mittelwert yi(~µx):
yi(~x)− µi =
n∑k=1
(xk − 〈xk〉)∂yi∂xk
+ · · · (6.134)
Element der Kovarianzmatrix Vpq:
Vpq[y] = E[(yp − µp)(yq − µq)] (6.135)
≈ E
(∑k
(xk − 〈xk〉)∂yp∂xk
)∑j
(xj − 〈xj〉)∂yq∂xj
(6.136)
= E
∑k,j
(xk 〈xk〉)(xj 〈xj〉)∂yp∂xk
∂yq∂xj
(6.137)
90
6.4. Nachtrag und Exkurs: Fehlerfortpflanzung
∂yi∂yk
: Element der Matrix B und daher
V [~y] = BV [~x]BT (6.138)
Spezialfall:Keine Korrelationen in ~x→ V [~x] ist diagonal, mit
Vii = σ2yi =
n∑k=1
(∂yi∂xk
)2
V [xk] (6.139)
=n∑k=1
(∂yi∂xk
)2
σ2xk
(6.140)
Merke: ~y sind im Allgemeinen auch dann korreliert, wenn ~x nicht korreliert sind.
Beispiel:Rotation eines Vektors um den Winkel θ:
~y = B~x , ~y =
(y1
y2
), ~x =
(x1
x2
)(6.141)
B =
(cos θ sin θ− sin θ cos θ
)(6.142)
BT =
(cos θ − sin θsin θ cos θ
)(6.143)
Kovarianzmatrix V [~x]: σ21 = V [x1], σ2
2 = V [x2]
V [~x] =
(σ2
1 cov(x1, x2)cov(x1, x2) σ2
2
)(6.144)
V [~y] = BV [~x]BT (6.145)
mit
V11 = σ2y1 = cos2 θσ2
1 + sin2 θσ22
+2 sin θ cos θcov(x1, x2) (6.146)
V22 = σ2y2 = cos2 θσ2
2 + sin2 θσ21
+2 sin θ cos θcov(x1, x2) (6.147)
cov(y1, y2) = V12 = V21 = sin θ cos θ(σ21 − σ2
2)
+(cos2 θ − sin2 θ)cov(x1, x2) (6.148)
cov(y1, y2) = 0, wenn
tan 2θ =2cov(x2, x2)
σ21 − σ2
2
(6.149)
91
6. Einfache statistische Methoden
→ Transformation in ein spezielles Koordianatensystem moglich, in dem die yiunkorreliert sind.
→ Auch bei unkorrelierten xi sind die yi im Allgemeinen korreliert.
6.5. Numerische Optimierung
6.5.1. Vorbemerkungen
Minimierung ohne Nebenbedingung
Minimiere F (~x), ~x ∈ R.F (~x) ist glatt ?, d.h. mindestens F ′ und F ′′ existieren.
~g(~x) = ~∇F (~x) =
∂F/∂x1
∂F/∂x2...
∂F/∂xn
(6.150)
~∇2F (~x) = H(~x) Hesse-Matrix (6.151)
=
∂2F∂x21
∂2F∂x1∂x2
· · · ∂2F∂x1∂xn
∂2F∂x2∂x1
∂2F∂x22
. . ....
.... . .
. . ....
∂2F∂xn∂x1
· · · · · · ∂2F∂x2n
(6.152)
Taylor-Entwicklung um xn:
F (~xk −∆~x) = F (~xk) + ~gTk ∆~x+1
2∆~xTHk∆~x+ . . . (6.153)
Notwendige Bedingung fur n = 1
d
dxF (x)
∣∣∣∣x=x∗
= 0 (6.154)
Maximum?, Wendepunkt?, Nebenminimum?
→ d2
dx2F (x)
∣∣∣∣x=x∗
> 0 (6.155)
(wenn F (x) an x∗ stetig ist.)
92
6.5. Numerische Optimierung
Notwendige Bedingung in n Dimensionen
d
dxiF (~x) = 0, fur alle i (6.156)
und H(~x∗) ist positiv definit, d.h. alle Eigenwerte sind großer als 0.
Spektralzerlegung von H:
H~u = λ~u, (6.157)
~u: Eigenvektor, λ: EigenwertH symmetrisch, n× n
n reelle Eigenwerte λi
n Eigenvektoren ~ui, orthogonal
Eigenvektoren seien normiert ~uTi · ~ui = 1
Bilde eine orthogonale Matrix U aus den Spaltenvektoren ~ui. Dann ist
D = UTHU =
λ1 0 . . . 0
0 λ2. . .
......
. . .. . . 0
0 . . . 0 λn
(6.158)
eine Diagonalmatrix mit den Eigenwerten λi.
– Es ist:
U−1 = UT (6.159)
H = UDUT =n∑j=1
λiUiUTi (6.160)
– Die Eigenwerte der inversen Hesse-Matrix H−1 sind invers zu denen von H,die Eigenvektoren sind identisch.
– Spektrale Norm einer symmetrischen Matrix: großter Eigenwert.
‖H‖ = λmax (6.161)
‖H−1‖ = λmin (6.162)
– Konditionszahl:
k = ‖H−1‖ · ‖H‖ =λmax
λmin(6.163)
(große Konditionszahl → numerisch problematisch)
93
6. Einfache statistische Methoden
Betrachte:
F (~x) = ~gT0 · ~x+1
2~xTH0~x (quadr. Funktion) (6.164)
und
~g(~x) = ~g0 +H0~x (Gradient.) (6.165)
Variiere die Funktion um den Punkt ~x∗, mit ~g(~x∗) = ~0 in Richtung von ~ui:
F (~x∗ + h~ui) = F (~x∗) +1
2λih
2 (6.166)
H: positiv definit λi > 0 Minimum um ~x∗.Da
~g(~x∗) = 0 (6.167)
wird aus
~g(~x) = ~g0 +H0~x (6.168)
H0~x∗ = −~g0 (6.169)
oder
~x∗ = −H−10 ~g0 (6.170)
→ Minimum kann in einem Schritt gefunden werden.
– Newton-Methode: quadratische Approximation der Funktion...
– Fur quadratische Funktionen f(x1, x2) (2Parameter):Hohenlinien konstanter Funktionswerte und Ellipsen, deren Achsen in Rich-tung der Eigenvektoren liegen.(Mehr als 2 Parameter: “Hyperellipsoide“.)
– Suche nach quadratisch approximierbarem Bereich:
→ lokales(?) Minimum
→ Sattelpunkt?Hesse-Matrix indefinit. Wahle zur Bestimmung der Suchrichtung eineLinearkombination von Eigenvektoren mit negativen Eigenwerten.
→ H positiv semidefinit?Mindestens ein λi = 0. Korrelationen? Problem schlecht konditioniert?
– Funktion nicht glatt?Wahl einer Suchmethode, die H nicht benutzt.
– log-Likelihood-Funktionen (spater!) sind in der Nahe der Losung meist qua-dratisch.
94
6.5. Numerische Optimierung
6.5.2. Eindimensionale Minimierung
Suchmethoden (genau ein Minimum!)
– Berechnung von Funktionswerten
– Einschließung des Minumums
– Verkleinerung des Intervalls
→ Methode des Goldenen SchnittsGegeben Intervall[x1, x3].Es wird durch den Goldenen Schnitt bei x2 geteilt, wenn
x3 − x2
x3 − x1=x2 − x1
x3 − x2. (6.171)
“Kleineres Intervall zu großerem, wie großeres zur ganzen Strecke.“
xt − x1
x3 − x1= τ = 0, 618 . . . Reduktionsfaktor (6.172)
Robuste Konvergenz, aber nicht ubermaßig schnell:
Abbildung 6.1.: Minimumsuche mittels Goldenem Schnitt
Nach n Schritten Reduktion um τn: z.B. τ10 ≈ 0, 00813 ≈ 1100 .
Newton-Methode
Suche nach Nullstellen in der Funktion:Iterationsvorschrift
xk+1 = Φ(xk) Φ(xk) = x− f(x)
f ′(x)(6.173)
Suche nach Nullstellen in der Ableitung
xk+1 = Φ(xk) Φ(xk) = x− f ′(x)
f ′′(x)(6.174)
→ Konvergenz gegen Fixpunkt x∗ (moglich)
– f ′′(x) = 0 → Schritt unendlich
– f ′′(x) < 0 → Konvergenz gegen Maximum
KonvergenzverhaltenIterationsmethode heißt lokal konvergent von der Ordnung p, wenn
|xk+1 − x∗||xk − x∗|p
< c, c > 0 (6.175)
95
6. Einfache statistische Methoden
mit Fixpunkt x∗. Wenn mindestens lineare Konvergenz
p = 1, c < 1 (6.176)
gilt, heißt die Methode global konvergent.
– lineare Konvergenz ist sehr langsam quadratische Konvergenz erwunscht.
– in der Newton-Methode ist
Φ(x) = x− f ′(x)
f ′′(x)(6.177)
und somit
Φ′(x) = 1− f ′′ 2(x)− f ′(x)f ′′′(x)
f ′′ 2(x)= −f
′(x)f ′′′(x)
f ′′ 2(x)(6.178)
an der Stelle x∗ ist f ′(x∗) = 0.
Φ′(x∗) = 0 (6.179)
zweite Ableitung:Φ′′(x∗) 6= 0 (6.180)
quadratisch konvergent (lokal).
Kombination der Methoden
– Suchmethoden:
∗ linear konvergent, robust
∗ keine Ableitung
– Newton:
∗ quadratische Konvergenz, kaum Fehlgehen
∗ Ableitungen
Kombination: Polynom-Interpolation:
– Beschreibung der Punkte durch ein Polynom 2. oder 3. Grades
– Mathematik des Newton-Verfahrens
xt = x2 −f ′(x2)
f ′′(x2)(6.181)
– symmetrische Teilung der Intervalle fuhrt auf
f ′(x2) ≈ f(x3)− f(x1)
2(x2 − x1)(6.182)
f ′′(x2) ≈ f(x3)− 2f(x2) + f(x1)
(x2 − x1)2(6.183)
96
6.5. Numerische Optimierung
→ Kann auf asymmetrische Teilung des Intervalls und damit zu Instabilitatfuhren.
→ Besser Teilung nach Goldenem Schnitt.
6.5.3. Mehrdimensionale Minimierung
Gittersuche
in n Dimensionen→ kn Funktionsberechnungen Vermehrung um Faktor 10, bei 10 Parametern Faktor 1010 Rechenoperationen mehr unpraktikabel
Monte Carlo Suche
a) Bestimmung der Unterteilung in einem Intervall in jeder der Dimensionenuber einen linearen Zufallsgenerator.
b) Wurfeln der Schrittrichtung vom kleinsten Funktionswert uber den Rich-tungskosinus
di =ri√∑ni=1 r
2i
(6.184)
Schrittweite ∆xi vorgegeben.
→ Gut fur die Erzeugung von Startwerten (manchmal).
Parametervariation
– Die n Parameter werden zyklisch jeweils eindimensional jeweils einem Opti-mierungsschritt unterworfen.
Abbildung 6.2.: Zyklische eindimensionale Minimierungen
→ Konvergiert (langsam) bei glatten Funktionen.
Simplex-Methode
– Schließe das Minimum durch einen k-dimensionalen Polyeder mit n Ecken ein.
– Sortiere die Funktionswerte nach Große.
– Ersetze den schlechtesten Punkt durch einen geeigneten in Richtung des Schwer-punktes variierten neuen Punkt.
97
6. Einfache statistische Methoden
→ Fallunterscheidung notig (Schwerpunkt 6= Minimum)
Newton-Methode
→ Wahle ~x0 als Anfangswert. Berechne F (~x0), setze Zahler k = 0.
→ Berechnung des Suchvektors ∆~xAnalog zum eindimensionalen Fall folgt aus der Bedingung, dass der Gradient= 0 sein soll:
~g(~xk + ∆~x) ≈ ~gk +Hk∆~xN = 0 (6.185)
→ Hk∆~xN = −~gk (6.186)
→ Eindimensionale MinimierungMinimiere F (~x) in der Richtung des Suchvektors.Betrachte dazu
f(z)︸︷︷︸1-dim.
= F (~xk + z∆~xN ) (6.187)
Minumum gegeben durch zmin!
→ IterationWahle das Minimum als Ausgangspunkt fur den nachsten Optimierungs-schritt:
~xk+1 = ~xk + zmin∆~xN (6.188)
Setze k = k + 1.
→ KonvergenztestIst die Approximation gut genug?
ja? → fertig
nein? → nachste Runde
zu viele Schritte? → aufgeben :-(
98
7. Spezielle Verfahren zur Datenanalyse
“Uber den Umgang mit a-priori-Wissen uber die Verteilungsfunktionen.“Probabilistenvs. Frequentisten
7.1. Die Maximum-Likelihood-Methode
7.1.1. Problemstellung
Vorgehen: n Messungen der Zufallsvariablen ~x, (dim ~x = 1 . . . k).Bekannt sei die Wahrscheinlichkeitsdichte fur die Messwerte ~x1 . . . ~xn bei gegebe-nem ~a: f(~x|~a). ~a: Parameter von denen die Wahrscheinlichkeitsdichte abhangt.
7.1.2. Aufgabe
Finde die beste Schatzung ~a aus den gegebenen Messdaten.
7.1.3. Ansatz
Bilde die Likelihood-Funktion
L(~a) = f(~x1|~a) · f(~x2|~a) · · · f(~xn|~a) =n∏i=1
f(~xi|~a) (7.1)
L(~a)
– ist ein Maß fur die Wahrscheinlichkeit, bei festliegenden Parametern ~a, denDatensatz ~x1 . . . ~xn zu nennen.
– ist keine Wahrscheinlichkeitsdichte.
7.1.4. Maximum-Likelihood-Prinzip
Die beste Schatzung von ~a maximiert die Funktion L(~a), d.h.
L(~a) = max, (7.2)
mit f(~x|~a) fur alle Werte von ~a auf 1 normiert:∫f(~x|~a) d~x = 1 ∀~a (7.3)
99
7. Spezielle Verfahren zur Datenanalyse
→ numerischer Aufwand wahrend aller Iterationsschritte!Maximum durch Differenzieren:
dL(~a)
d~a= 0
∂L(~a)
∂ai= 0 fur i = 1 . . . k (7.4)
7.1.5. Log-Likelihood-Funktion
l(~a) = ln(L(~a)) =n∑i=1
ln(f(~xi)|~a) (7.5)
ln ist monoton → Maximum von l(~a) und L(~a) an der gleichen Stelle.Ein Parameter:
l(a) = ln(L(a)) =
n∑i=1
ln(f(xi|a)) = max (7.6)
→ dl(a)
da= 0 (7.7)
k Parameter:∂l(~a)
∂ai= 0 i = 1 . . . k (7.8)
= Minimierung der negativen Log-Likelihood-Funktion:
F (~a) = −l(~a) = −n∑i=1
ln(f(~xi|~a)) (7.9)
Eigenschaften der Maximum-Likelihood-Methode
+ konsistent (unverzerrt) (siehe 7.3)
+ nicht immer erwartungstreu
+ beides aber fur n→∞+ effizient
- großer Rechenaufwand
- a priori Kenntnis der Form der Wahrscheinlichkeitsdichte notig
- Uberprufung notig, dass die Daten uberall mit den Parametern ~a vertraglichsind. . .bei mehreren Dimensionen in allen Teilintervallen. . .
7.1.6. Beispiele
1. Zerfallsintervallverteilung eines TeilchensGegeben durch
f(x|a) =1
2(1 + ax)︸ ︷︷ ︸
zufallig normiert! s.u.
mit x = cosϑ (7.10)
100
7.2. Fehlerbestimmung bei der Maximum-Likelihood-Methode
∫ 1
−1
1
2(1 + ax) dx =
1
2x+
1
4ax2
∣∣∣∣1−1
=1
21 +
1
21 +
1
4a− 1
4a = 1 (7.11)
n Werte fur xi gemessen, a gesucht, −1 ≤ x ≤ 1
F (a) = −n∑i=1
ln
(1
2(1 + axi)
)(7.12)
2. Gauß-Verteilung; Mittelwert a
f(xi|a) =1
√aπσi
e− (xi−a)
2
aσ2i (7.13)
F (a) = konst +1
2
n∑i=1
(xi − a)2
σ2i
(7.14)
dF (a)
da
!= 0 −
n∑i=1
xi − aσ2i
= 0 (7.15)
mit wi = 1/σ2i ist
a =
∑ni=1 xiwi∑ni=1wi
(7.16)
a ist das gewichtete Mittel aller xi
3. Vergleich zur Methode der kleinsten QuadrateAus Beispiel 2. folgt: Wenn die Daten Gaußverteilt sind, ist die Minimie-
rung von F (a) gleichbedeutend zur Minimierung von S(a) =∑n
i=1(xi−a)2
σ2i
, da
F (a) = konst + 12S(a).
7.2. Fehlerbestimmung bei derMaximum-Likelihood-Methode
7.2.1. Ein Parameter
Zeichne F (a) gegen a:
Bei 12
22
42
↓ ↓ ↓1σ 2σ 1σ
-Zunahme von F (a)
101
7. Spezielle Verfahren zur Datenanalyse
Da Likelihood-Funktion oft naherungsweise Gauß -verteilt:h→∞: Likelihood → Gauß und Varianz → 0.
F (a) = F (a) +1
2
d2F
da2(a− a)2 + . . . 1 (7.17)
L(a) ≈ const · e−12d2Fda2
(a−a)2 (7.18)
= const · e−(a−a)2
2σ2 (7.19)
in der Nahe des Maximums Gauß formig.
Koeffizientenvergleich: σ(a) =
(d2F
da2
∣∣∣∣a
)−1/2
(7.20)
→ neg. Log-Likelihood-Funktion: Parabel, 2.Ableitung = const.
→ F (a± rσ) = F (a) +1
2r2 (7.21)
7.2.2. Mehrere Parameter
Parameter a1, . . . , an = ~a.
L(~a) =n∏i=1
f(xi, a1, a2, . . . , an) (7.22)
Entwicklung der negativen log-Lokelihood-Funktion um ~a:
F (~a) = F (~a) +1
2
n∑i,k
∂2F
∂ai∂ak(ai − ai)(ak − ak) + . . . (7.23)
= F (~a) +1
2
n∑i,k
Gik(ai − ai)(ak − ak) + . . . (7.24)
Gik =∂2F
∂ai∂akV = G−1 = Kovarianzmatrix
– Gik=Hesse-Matrix, kann fur Optimierung genutzt werden
– Formal exakt fur große n, sonst Naherung
– 2 Parameter:
- Konturlinien als Linien gleicher Likelihood
F (~a) = F (~a+1
2r2 ↔ F (~a) +
1
2(7.25)
- Abweichung vom asymptotischen Verhalten→ asymmetrische Fehler.
1Wenn hohere Terme nicht verschwinden, aufwendigeres Verfahren: insbesondere sind Fehler nicht mehrsymmetrisch.
102
7.3. Die Maximum-Likelihood-Methode, Eigenschaften
7.3. Die Maximum-Likelihood-Methode, Eigenschaften
7.3.1. Konsistenz
a sei die beste Schatzung des Parameters a.a0 sei der wahre Wert dieses Parameters.Wenn
limn→∞
a = a0, (7.26)
heißt die Schatzung konsistent.Dies kann fur die Maximum-Likelihood-Methode gezeigt werden. (→ Blobel) X
7.3.2. Erwartungstreue?
Schatzung heißt erwartungstreu, wenn
E[a] = a0 (7.27)
Anschaulich: wenn der Mittelwert der Maximum-Likelihood-Funktion = dem Er-wartungswert, oder die Maximum-Likelihood-Funktion im Intervall [−pσ, pσ] sym-metrisch ist.Da die Maximum-Likelihood-Funktion durch Ableitungen g(a) verzerrt werdenkann (asymmetrische Fehler), gilt die Erwartungstreue nur fur n→∞.
7.3.3. Gaußahnlichkeit
Die Maximum-Likelihood-Funktion nahert sich asymptotisch der Gauß-Funktionan.
7.3.4. Varianz
Die Varianz
σ(a) =
(d2F
da2
∣∣∣∣a
)−1/2
(7.28)
ist die”kleinstmogliche“. Daher ist die Maximum-Likelihood-Methode asymptotisch
effizient.
7.4. Bayesische Statistik
Bisher:
– alle Aussage dienen dazu, Wahrscheinlichkeitsintervalle anzugeben,
– der Erwartungswert dient als Referenzpunkt in dem Intervall,
103
7. Spezielle Verfahren zur Datenanalyse
– a priori-Wissen uber die Form von Verteilungen fließt nur in der Maximum-Likelihood-Methode ein.
Vergleiche nun:
– gekoppelte Wahrscheinlichkeiten
P (A und B) = P (A) · P (B|A) (7.29)
= P (B) · P (A|B) (7.30)
→ Bayes Theorem
P (A|B) = P (B|A)P (A)
P (B)(7.31)
– mit dem vorliegenden Problem:
f(a|x) = f(x|a) · f(a)
g(x)(7.32)
a: Parameter
x: Daten
f(x|a): bedingte Wahrscheinlichkeit fur die Daten bei gegebenem a
”Monte-Carlo-Simulation“
f(a|x): bedingte Wahrscheinlichkeit fur die Parameter a bei gegebenen Daten
g(x): liegt nach Messung fest X
f(x|a): kann berechnet werden X
f(a): problematisch, da a eigentlich nur einen Wert haben.
”E inzug der Vorurteile in die Statistik“
z.B. bei Reflektionsproblemen auch berechenbar/motivierbar
f(a):”Prior“
f(a|x):”Posterior“
”Gedankenspiel“
Seien µ der Mittelwert einer Wahrscheinlichkeitsdichte und x der Mittelwert derdazugehorigen Stichprobe.Zentraler Grenzwertsatz:Wahrscheinlichkeitsdichte von x ist naherungsweise Gauß funktion
N(µ, σ) mit σ2: Varianz (7.33)
Prior: f(µ).
⇒ f(µ|x) = N(µ, σ) · f(µ) |ohne Normierung (7.34)
104
7.4. Bayesische Statistik
Maximum-Likelihood:
df(µ|x)
dµ= 0 =
dN(µ, σ)
dµ· f(µ) +N(µ, σ)
f(µ)
dµ(7.35)
d lnN(µ, σ)
dµ= − 1
f(µ)· df(µ)
dµ=
(x− µ)
σ2(7.36)
Sei der Prior f(a) gleichverteilt in a oder . . . (konnte z.B. auch gleichverteilt in a2
oder log a sein), dann ist µ = x.Allgemein gilt
µ = x+ σ2 · 1
f(µ)· df(µ)
dµ= x+ σ2d ln f(µ)
dµ(7.37)
→ Einfluss des Priors ∼ d ln f(µ)dµ .
Beispiel: Katastrophenwahrscheinlichkeit100 Raketenstarts waren erfolgreich.Gesucht: Wahrscheinlichkeit eines Fehlstarts beim 101. VersuchZur Beantwortung notig:
Wissen uber p=Wahrscheinlichkeit fur Fehlstart bei einem Versuch
Obere Grenze p0 mit 95% CL.
→ Binomialverteilung, Wahrscheinlichkeit fur”0“-Ereignisse = 100%-95%=5%.
P (100, 0) = p00 · (1− p0)100 = 5% (7.38)
→ 100 ln(1− p0) = ln(0, 05) (7.39)
→ p0 = 0, 03 (7.40)
”Ende der klassischen Statistik“.
Rechnung dazu:
Binomial:
Versuche=n︷ ︸︸ ︷(1000
)︸ ︷︷ ︸
Fehlstarts=k
p0(1− p)100−0 = P (7.41)
Poisson:(100p)0
0!︸ ︷︷ ︸=1
·e−100p = 0, 05 (7.42)
⇒ −100p = ln 0, 05 ≈ −3 (7.43)
⇒ p = 0, 03 (7.44)
Schatzwert statt Grenze → Bayesischer Ansatz→ Wahrscheinlichkeit, dass etwas passiert:
Poisson(1): p =λ1
1!e−λ, λ = n · p (7.45)
105
7. Spezielle Verfahren zur Datenanalyse
Prior gleichverteilt:
f(p, 1) = λ · e−λ · const (7.46)
df(p, 1)
dp= λ · e−λ − 1e−λ (7.47)
= e−λ(λ− 1)!
= 0 (7.48)
λ = 1 p = 1% (7.49)
7.5. Entfaltung
7.5.1. Problem
Gesucht: Verteilung in einer physikalisch relevanten Große x (z.B.”Energie“):
f(x).
Gemessen: Verteilung einer damit korrelierten Große y (”z.B. Zahl der Hits
im Detektor“): g(y). → Auflosung
– unerwunschte Ereignisse im Detektor, die durch Rauschen des Meßinstru-ments ausgelost werden → Untergrund (Background)
– nicht alle erwunschten Ereignisse werden vom Detektor registriert → Akzep-tanz
Betrachte die Fredholmsche Integralgleichung 1. Art :
g(y) =
∫ b
aA(y, x)f(x) dx (7.50)
x, f(x), a, b: gesuchte Verteilung in x im Intervall [a, b]
y, g(y): gemessene Verteilung
A(y, x): Kern; gibt die Wahrscheinlichkeit an, y zu messen, wenn der wahre Wert xist. →
”Messgerat-Eigenschaften“
Hinzu kommt:
ε(y): statistischer Fehler
b(y): Untergrund
g(y) =
∫ b
aA(y, x)f(x) dx+ b(y) + ε(y) (7.51)
Sehr gute Auflosung → A(y, x) = A′(x) · S(x, y) nur noch Akzeptanzkorrektur notig.
106
7.5. Entfaltung
Allgemeiner Fall
Auflosung ist schlecht
→ Rekonstruktion von f(x); g(y) =Entfaltung
→ Schlecht konditioniertes Problem
- oszillierende oder sinnlose Losungen
→ Ausweg”
regularisierte Entfaltung“
Bemerkungen
Ubereinstimmung von g(y)gemessen mit g(y)Monte Carlo ist nur ein Konsistenztest.Keine Aussagen uber Fehler; Losung nicht unbedingt eindeutig; kann ebenfallsoszillieren.
Ubergang zu gebinnten (diskreten) Variablen(notwendig fur numerische Problemlosung)
gi =m∑j=1
Aijfj + bi + εi (7.52)
~g = A~f +~b+ ~ε (7.53)
~x, ~f : m-dim. Vektoren~g, ~b, ~ε : n-dim. Vektoren
}Histogramme
A: Kern → n×m Transfermatrix
Es sei n = m A quadratisch, es folgt fur kleine statistische Fehler
~f = A−1( ~tildeg −~b) (7.54)
Untergrundfreie Messung:
E[~f ] = E[A−1~g] = A−1E[~g] = A−1A~f = ~f (7.55)
→ Schatzung ist konsistent.Fehlerfortpflanzung:
V [~f ] = A−1V [~g]A−1T (7.56)
Wegen der Oszillationen der Losungen (schlechte Kondition) Verfahren in dieserForm nicht befriedigend. → Regularisierung
107
7. Spezielle Verfahren zur Datenanalyse
7.5.2. Akzeptanzkorrektur
Akzeptanz P : Wahrscheinlichkeit, dass bei Vorliegen des wahren Wertes x die Variable ygemessen werden kann.
Messung sehr genau:
Transfermatrix ist diagonal mit Ajj = Pj
gemessene Verteilung: gj = Pjfj
Ursache:
– z.B. Geometrie
→ z.B. Raumwinkelkorrektur
– Detektoreigenschaften (Trigger)
→ Monte Carlo
Akzeptanz:
Pj =Nj rekonstruiert
Nj generiert=Nj r
Nj g(7.57)
Fehler
σj, rel =∆PjPj
=
√1
Nj g
1− PjPj Binomialverteilt (7.58)
σj =
√Pj(1− Pj)
Nj g=
1
Nj g
√Nj r ·Nj verworfen (7.59)
Korrigierte Werte
f =gjpj
(7.60)
σ(fj) =
√gj
pj
√1 + gjσ2
j, rel (7.61)
Verzicht in Regionen kleiner Akzeptanz! Fehler, atypische!
7.5.3. Diskretisierung
~g(y) =
∫ b
aA(y, k)f(x) dx+ b(y) + ε(y) (7.62)
→ ~g = A~f +~b+ ~y (7.63)
– Genauigkeitsverlust (i.d.R. unvermeidbar)
108
7.5. Entfaltung
Parametrisierung
Entwicklung in Pj
f(x) =
m∑j=1
ajPj(x), Pj(x): Basisfunktionen (7.64)
Dann ist ∫ b
aA(y, k)f(x) dx =
m∑j=1
aj
(∫ b
aA(y, x)Pj(x) dx
)(7.65)
=m∑j=1
ajAj(y) (7.66)
→ ~a: m-Vektor
mit
Aj(y) =
∫ b
aA(y, x)Pj(x) dx (7.67)
g(y) =
m∑j=1
ajAj(y) + b(y) + ε(y) (7.68)
Darstellung
Darstellung aller y-abhangigen Funktionen durch Histogramme mit den Grenzeny0, y1, · · · , yn.
gj =
∫ yi
yi−1
g(y) dy ~g : n-Vektor (7.69)
Aij =
∫ yi
yi−1
Aj(y) dy A : (n×m)-Matrix (7.70)
bi =
∫ yi
yi−1
b(y) dy ~b : n-Vektor (7.71)
~g = A~a+~b (7.72)
A:
– Spalte Aj : Histogramm in y fur f(x) = PJ(x)
– Elemente werden im Monte Carlo erzeugt
– Wahle beliebige Verteilung der wahren Werte f0(x)
– Generiere dazu ein y
109
7. Spezielle Verfahren zur Datenanalyse
– Addiere ein Gewicht proportional zu Pj(x) zu dem Histogramm Aj(y)
Basisfunktionen parametrisieren eine Dichte → Pj(x) ≥ 0
Wahl der Basisfunktionen:
- B-Splines mit den Eigenschaften
Pj(x) ≥ 0m∑j=1
Pj(x) = 1 (7.73)
Koeffizienten aj ≥ 0 werden bestimmt.
Berechnung der fk
fk =1
xk − xk−1
∫ xk
xk−1
f(x) dx =1
xk − xk−1
m∑j=1
aj
∫ xk
xk−1
Pj(x) dx (7.74)
f0(x) im Prinzip frei wahlbar, aber sinnvolle Wahl: nahe am erwarteten Ergebnis(Rechenzeit).
f(x) = f0(x) · f1(x) (7.75)
mit: f1(x) nur schwach x-abhangig.
g(y) =
∫ b
aA(y, x)f(x) dx (7.76)
=
∫ b
aA(y, x)f0(x) · f1(x) dx
=
∫ b
aA0(y, x)f1(x) dx (7.77)
7.5.4. Anwendung
Entfaltung in zwei Bins:
~g = A~f ; ~g =
(g1
g2
); A =
(1− ε11 ε12
ε21 1− ε22
)(7.78)
~f =
(f1
f2
)(7.79)
Fehler in g1, g2 unabhangig
V [~g] =
(g1 00 g2
)(7.80)
~f = B~g mit B = A−1 (7.81)
V [~f ] = BV [~g]BT (7.82)
110
7.5. Entfaltung
Seien nun die Akzeptanz 100% und alle ε = ε11 = ε22 = ε12 = ε21.Entfaltungsproblem hier: Ereignisklassifikation
A−2 = B =1
1− 2ε
(1− ε −ε−ε 1− ε
)(7.83)
f1 =1
1− 2ε[(2− ε)g1 − εg2] (7.84)
f2 =1
1− 2ε[−εg1 + (1− ε)g2] (7.85)
V [~f ] =1
(1− 2ε)2
((1− ε)2g1 + ε2g2 −ε(1− ε)(g1 + g2)−ε(1− ε)(g1 + g2) (1− ε)2g2 + ε2g1
)(7.86)
↪→ Divergenz fur ε→ 0, 5∧= keine Messung. (7.87)
Zahlenbeispiel:
g1 = 100± 10 g2 = 81± 9
Vermischwahrscheinlichkeit ε = 0, 1
f1 = 102, 4± 11, 3 f2 = 78, 6± 10, 2
Vermischwahrscheinlichkeit ε = 0, 4
f1 = 138, 0± 35, 0 f2 = 43, 0± 33, 6
7.5.5. Entfaltung (ohne Regularisierung)
Diskretisierung → ~g = A~a+~bAnzahl von Eintragen/Bin folgt normierter Poisson-Verteilung neg. Log-Likelihood-Funktion:
F (~a) =∑
(gi(~a)− gmi ln gi(~a)) (7.88)
m: Zahl der Bins in x Iterative Bestimmung von ~a
F (~a) = F (~a) + (~a− ~a)T · ~h+1
2(~a− ~a)TH(~a− ~a) (7.89)
h: Gradient; H: Hesse-Matrix
Korrektur: ∆~a = −H−1~h (7.90)
~a = ~a+ ∆~a bis Konvergenz, (7.91)
dann V [~a] = H−1 (7.92)
Problem ist immer noch instabiles Verhalten der Losungen Oszillationen imParameterraum.
111
7. Spezielle Verfahren zur Datenanalyse
7.5.6. Problemanalyse
Dazu: Diagonalisierung der Hesse-Matrix :
H = UHDUTH (7.93)
D = UTHH UH∣∣UHorthogonal mit UTHUH = 1 (7.94)
Eigenschaften von D:
– Diagonalelemente sind reelle Eigenwerte von H.
– Dii > 0, wenn H positiv definit.
– Dii konnen in absteigender Reihenfolge angeordnet werden. Typischer Be-reich: einige Großenordnungen
Definiere dann Digonalmatrix D1/2 mit
D1/2D1/2 = D, desgleichen fur D−1/2 (7.95)
Transformiere nun~a = UHD
1/2~b (7.96)
zu minimieren ist nun
F (~b) = F0︸︷︷︸2 konst. Beitr.
+~bT ~hb +1
2~bT ~b (7.97)
mit~hb = D−1/2 UTH
(−H ~a+ ~h
)(7.98)
Losung fur den transformierten Vektor
~b = −~hb (7.99)
mit Kovarianzmatrix
V [~b] = 1 (!) (7.100)
Darstellung des Ergebnisses zu der neuen Basis
f(x) =m∑j=1
bjp′j(x) (7.101)
p′j(x): passende Linearkombination aller pjRucktransformation des Ergebnisses
~a = UHD−1/2~b =
m∑j=1
(1
Djj
)1/2
bj ~uj︸︷︷︸Eigenvektor
(7.102)
Djj waren nach Große angeordnet,
112
7.5. Entfaltung
mit Null vertraglich
konnen aber nach Multiplikation mit bj großen Beitrag zu ~a liefern
”Oszillationen“ im Resultat
Abschneiden?
”Gibbsches Phanomen“ →
”Uberschwinger“ → Fluktuationen
Regularisierung
7.5.7. Regularisierung
Abbildung 7.1.:”Gibbsches Phanomen“
Ziel: vorsichtiges Abschneiden der insignifikanten Parameter.Dazu: vorsichtige Verzerrung des Problems.Betrachte die Krummung von f(x) als Maß fur die Fluktuationen:
r(~a) =
∫ (f ′′(x)
)2dx (7.103)
Parametrisierung durch kubische B-Splines
→ r(~a) = ~aT C~a, (7.104)
mit C =
2 −3 0 1 0 0 . . .
−3 8 −6 0 1 0. . .
0 −6 14 −9 0 1. . .
1 0 −9 16 −9 0. . .
0 1 0 −9 16 −9. . .
0 0 1 0 −9 16. . .
.... . .
. . .. . .
. . .. . .
. . .
(7.105)
Minimiere nun
R(~a) = F (~a) +1
2τ r(~a) (7.106)
= F (~a) +1
2τ~aT C~a (7.107)
τ → 0: Einfluss der Regularisierung verschwindet
τ →∞: lineare Losung
Wie ist τ zu wahlen, damit der Einfluss der Verzerrung klein bleibt?
113