Zusammenfassung und Konzept
Computergestützte Datenauswertung
Günter Quast, Andreas Poenicke SS 2016
Neue Veranstaltung
„Rechnergestützte Datenauswertung mit Computerübung“
mit 1 SWS Vorlesung und 1 SWS Übung erstmalig im SS16 (als verpflichtender Kurs Schlüsselkompetenz, 2. Semester im Bachelor Physik) G. Quast, A. Poenicke & viele Assistenten
Ziele: – einheitliche Arbeitsumgebung für alle Studierenden (für Windows, Linux & Mac OS, CIP-Pool, basierend auf python/matplotlib)
– grafische Darstellung von Funktionen und Daten
– Anpassung von Funktionen (=Modellen) an Messdaten & grundsätzliches Verständnis des theoretischen Hintergrunds
– Modellierung von experimentellen Daten („Monte Carlo“)
– Beherrschung des python-Pakets kafe zur Funktionsanpassung
1.1: Aufsetzen der Arbeitsumgebung
1.2: Kennenlernen der Arbeitsumgebung
1.3: Kennenlernen von Python
2.1: Kennenlernen von Python (2)
2.2: Arbeiten mit numpy
2.2: Arbeiten mit matplotlib
3.1: Darstellen von Funktionen
3.2: Berechnung statistischer Größen
3.3: Funktionen von Zufallszahlen
4.1: Würfelspiel
4.2: Zentraler Grenzwertsatz
4.3: Parameterschätzung: Resonanzkurve mit Messdaten
5.1 Signalverarbeitung: Frequenzbestimmung
5.2: Korrelation
5.3: Anpassung von Modellen an Daten mit kafe
6: Datenauswertung
Frequenzbestimmung beim Federpendel und Messung der Erdbeschleunigung
Empfohlene Software-Grundausstattung
Software-Pakete sind quelloffen (Open Source) d.h. frei unter Linux, MS Windows, Mac OSX und auf dem CIP-Pool verfügbar
(einfache) Entwicklungsumgebung für python IDLE Integrated Dvelopment and LEarning
interaktives Arbeiten mit python: IPython http://ipython.org/notebook.html
grafische Darstellung mit python: matplotlib
python-Bibliotheken für wissenschaftliches Arbeiten: numpy und scipy
Details s. CgDA, V01b_intro.pdf http://www.ekp.kit.edu/~quast/CgDA ( user: Students, password: only )
enthält auch Hinweise und ein vorkonfiguriertes Gesamtpaket für Win10, sowie eine „virtuelle Maschine“
http://www.ekp.kit.edu/~quast/CgDA/CgDA.html
Computergestützte Datenauswertung
Alle Dateien zu den gezeigten Beispielen (und noch viele mehr)
finden Sie unter dem Link zur Vorlesung
Spezielle Python-Beispiele zum Start ins Praktikum
http://www.ekp.kit.edu/~quast/CgDA/PhysPrakt
das python-Modul PhyPraKit.py enthält
1. Dateneingabe aus Text-Dateien
2. Signalprozessierung (Filter, Fourier-Transf., Autokorrelation)
3. Berechnung des gewichteten Mittelwerts
3. Histogramm-Tools barstat() statistische Information aus Balkendiagramm histstat() statistiche information aus 1d-Histogramm hist2dstat() statistische Information aus 2d-histogram profile2d() "profile plot" für 2d-Daten chi2p_indep2d() chi2-Test auf Unabhänggkeit zweier Datensätze
4. (lineare) Regression kRegression() mit kafe, x-, y- und and korrelierte Fehler
5. Erzeugung von simulierten Datensätzen smearData() generateData()
s. Beispiel-Scripts test_readColumnData.py
test_readCSV.py
test_readtxt.py
test_convolutionFilter.py
test_autoCorrelation.py
test_histogram.py
test_linRegression.py
test_kRegression.py
test_generateData.py
– Funktionsanpassung – Software-Paket kafe – Virtuelle Maschine
python ist eine
– moderne
– objektorientierte
– interpretierte
– intuitive
– leicht zu erlernende
– effiziente Skript- und Programmier- Sprache
mit
– Unterstützung praktisch aller Funktionen des Betriebssystems,
– eigenen Modulen insb. auch für wissenschaftliches Rechnen,
– Anbindung von Bibliotheken wichtiger Programmpakete,
– Einbindung von selbst-programmierter Funktionalität.
& friends:
- numpy
- scipy
- matplotlib
python ist verfügbar für Windows, Mac und Linux: Platform-übergreifend !
Wozu Statistische Datenanalyse ?
Die Aufgabe des Wissenschaftlers:
Vergleich von Modellen mit der Wirklichkeit
Frage: Passt das (theoretische) Modell zu den Beobachtungen (Messungen) ?
- wenn nein: Modell verwerfen oder verbessern - wenn ja: freie Modellparameter bestimmen
–– gutes Modell
–– schlechteres Modell
– Messungen sind immer mit statistischen Ungenauigkeiten behaftet (Rauschen, Ablese- oder Digitalisierungs- genauigteit, Eichung, ...)
– in der Quantenphysik sind die relevanten Modellvorhersagen und -Größen selbst Parameter von Zufallsverteilungen (mittlere Lebensdauer eines Zustands, Erwartungswert eines Operators, ...)
m = w + z
gemessener Wert
wahrer Wert
zufälliger
Beitrag
z kann viele Ursachen haben:
- zufälliger Beitrag zum Messwert („Rauschen“) → „statistische Unsicherheit“
- Genauigkeit des verwendeten Messinstruments → „systematische Unsicherheit“
- mitunter gibt es auch eine Unsicherheit auf den „wahren“ Wert, den man oft z zuschlägt → „theoretische Unsicherheit“
- Fehler im Messprozess – sollten nicht passieren !
z = zstat + zsyst + ztheo
Darstellung als
„Häufigkeitsverteilung“
oder „Histogramm“
(engl. frequency distribution“, „histogram“)
Visualisierung der Messreihe
(typisch wie z.B. mit MS Excel)
Mittelwert: 79.954
Standardabweichung: 0.522
Unsicherheit auf Mittelwert: 0.074
Statistische Information:
Anm.: Standardabweichung beschreibt „Breite der Verteilung“
Beispiel Strom-Spannungsmessung (1)
Messreihe:
I
U
+
–
Draht
?
Grafische Darstellung mit
(Millimeter -) Papier und Bleistift
oder Tabellenkalkulation beherrschen Sie alle mit ein wenig Python-Code wird es komfortabler ...
Beispiel Strom-Spannungsmessung (2)
# Measurements: Ohm's Law
# Prec.: U: +\-0.015V, I: +\-0.015A
# U[V] I[A]
0.15 0.14
0.25 0.22
0.40 0.36
0.50 0.48
0.60 0.53
0.70 0.64
0.75 0.65
0.80 0.74
0.90 0.81
Text-Datei
DataOhmsLaw.dat
import numpy as np
import matplotlib.pyplot as plt
#read columns from file
infile="DataOhmsLaw.dat"
U, I = np.loadtxt(infile,unpack=True)
# plot the results
fig, ax=plt.subplots(1,1)
ax.errorbar(U, I, fmt='o')
ax.set_xlabel('U')
ax.set_ylabel('I(U)')
plt.show()
Python Code
Grafik
Frage:
Gilt das Ohm'sche Gesetz ?
I=U/R = G R für
R bzw. G = const ?
Messreihe:
Häufigkeitsverteilung und Wahrscheinlichkeitsdichte
Beispiel Messen:
Häufigkeit von N Messergebnissen
in bestimmten Intervallen („Bins“)
Erwartungswert
der Verteilung
Intervallgrenzen („bin edges“)
Erwartete
Verteilung
Für eine große Anzahl von
Messungen nähert sich die
Häufigkeitsverteilung der
erwarteten Verteilung immer mehr an
Häufigkeitsdefinition der Wahrscheinlichkeit:
hi
bin i
Charakterisierung v. Verteilungen: Formeln
Stichprobe
diskrete Verteilung
kontinuierliche Vert.
Erwartungswert
Varianz
* „Bessel-Korrektur“: N → N-1 vermeidet Verzerrung
Plausibilitätsargument: für N=1 ist V nicht definiert !
*
Schiefe
Kurtosis
γ2 = 0 für Gauß-Vert.
+ höhere ...
nennt man das n-te Moment der Verteilung
Eine Verteilung ist über ihre Momente vollständig charakterisiert
Gauss ↔ Poisson ↔ Binomial - Verteilung
Binomial
n, p
Gauß
μ, σ
Poisson
μ
n→ ∞
p→ 0 np=μ
μ→ ∞
N(x; 0, 1)
N(x; 0, 2)
N(x; -1, 2)
Zentraler Grenzwertsatz oder warum sind Messfehler Gauß-verteilt ?
Im Grenzfall von großen N ist die Summe von N unabhängigen
Zufallszahlen eine Zufallszahl, die einer Gauß-Verteilung folgt.
xi aus beliebiger Verteilung mit
Mittelwert μi und endlicher Varianz σi
x ist gaußverteilt mit
Erwartungswert und Varianz
unter recht schwachen Annahmen (Lyapunov-Bedingung) gilt:
A
B
Projektion auf y-Achse
Projektion auf x-Achse
Randverteilungen
engl.: marginal distribution
kumulativ
Normierung
kumulativ
Definition Kovarianz:
Normiere so, dass Diagonalelemente = 1
⇒ Korrelationskoeffizenten :
Streudiagramme von zwei korrelierten Zufallsvariablen
Korrelationskoeffizienten ρij sind
anschaulicher als Kovarianzen !
x, y
- unkorreliert ≠ unabhängig
- korreliert ≠ kausal verknüpft
- Korrelation kann zufällig sein,
* x ⇒ y möglich
* y ⇒ x möglich
* z ⇒ (x, y) möglich
Kovarianzmatrix von korrelierten Messungen
Für n Messwerte mi erhält man die Kovarianz- Matrix, indem man die Varianz der gemeinsamen
Unsicherheit zg, (σg)², in die Nebendiagonale setzt. Die Diagonalelemente sind die Varianzen der Gesamtunsicherheiten, (σi)² + (σg)².
Die vollständige Kovarianzmatrix sieht so aus:
Mehrere (n) Messwerte mi , jede Messung hat - eine individuelle, zufällige Unsicherheit zi
- sowie eine gemeinsame Unsicherheit zg
→ mi = w + zi + zg
ein gemeinsamer wahrer Wert und n+1 Zufallszahlen (eine allen Messungen gemeinsame sowie n individuelle)
Multidimensionale Gaußverteilung
μi=0, σi=1, ρ=0.5
2d Gauß
allg. Kovarianzmatrix mit Korrelationskoeffizienten ρij
2-dimensional
Script: plot3dGauss.py
Kontur konstanter Wahrscheinlichkeitsdichte („Höhenlinie“) bei 2d-Gaußverteilung
ist eine Ellipsengleichung
gebräuchliche Darstellung bei Problemen mit vielen Variablen:
Kovarianz-Ellipsen von Variablen-Paaren
x2
x1
Kovarianz-Ellipse
Winkel zwischen x-Achse und Haupt-
achse der Ellipse hängt von ρ12 ab:
Fläche innerhalb der 1-σ Ellipse
entspricht ~39% Wahrscheinlichkeit
σ2
σ2
σ1
σ1
Monte Carlo-Methode zur Simulation
Die Monte-Carlo-Methode (auch „MC-Simulation“) - abgeleitet vom Namen der durch ihr Spielkasino berühmten Stadt Monte Carlo - ist eine auf (Pseudo-)Zufallszahlen basierende numerische Methode zur
Grundsätzliche Vorgehensweise:
Zufallskomponente z bei der Gewinnung empirischer Daten ebenfalls mit Zufallszahlen modellierbar:
m = w + z
gemessener Wert
wahrer Wert
zufälliger
Beitrag
Funktion (=Modell) des wahren Wertes, y = f(x), und
Modellierung der Messwerte mit MC:
my = f( wx + zx ) + zy
zx und zy entsprechen den Unsicher- heiten der Messgrößen x und y
Reale Daten: Digitale Abtastung („sampling“)
Digitale Messgeräte tasten Signale üblicherweise
regelmäßig zu festen Zeiten ti mit der Abtast- oder
Sampling-Rate fs ab.
Nummer der Messung
Script sampling.py
Datenexport
in Felder mit
- Zeitpunkten ti ( t[0, …,N-1] )
- Messwerten vi ( v[0, …,N-1])
Grafische Darstellung als Markierungen, evtl.
– verbunden durch gerade Linien „lineare Interpolation“
– verbunden durch Kurven, z.B. „Spline-Interpolation“ ( kubische Splines: stetig differenzierbar aneinander anschließende Polynomstücke 3. Grades )
Beispiel: Rohdaten einer Wellenform
Time,Channel A (ms),(V) -0.34927999,-0.00045778 -0.34799999,-0.00045778 -0.34671999,-0.00045778 -0.34543999,-0.00045778 -0.34415999,-0.00045778 -0.34287999,-0.00033570 -0.34159999,-0.00018311 -0.34031999,-0.00018311 -0.33903999,-0.00018311 -0.33775999,-0.00003052 -0.33647999,0.00006104 -0.33519999,0.00006104 -0.33391999,0.00006104 -0.33263999,0.00006104 -0.33135999,0.00021363 -0.33007999,0.00021363 . . . 9.63983995,0.02642903 9.64111995,0.02655110 9.64239995,0.02655110 9.64367995,0.02655110 9.64495995,0.02655110 9.64623995,0.02655110 9.64751995,0.02655110 9.64879995,0.02642903 9.65007995,0.02612384 9.65135995,0.02584918 9.65263995,0.02557451 9.65391995,0.02526933 9.65519995,0.02487259
7816 Wertepaare exportiert aus USB-Oszilloskop PicoScope im
csv-Format (comma separated values)
importiert in numpy-arrays t und v,
dargestellt mit plt.plot(t,v)
Datei
Welleenform.csv
Script:
test_readPicoscope.py
Messgeräte (auch „Datenlogger“) und einige Handy-Apps (z.B. phyphox)
nutzen einfache Datenformate in Text-Form:
Time,Channel A (ms),(V) -0.349,-0.000458 -0.348,-0.000458 -0.347,-0.000458 . . .
Beispiel PicoScope, „Comma Separated Values“ (CSV)
Kopfzeilen mit sog. „Meta-Daten“
die eigentlichen Daten als Dezimal-
zahlen, in Spalten durch „,“ getrennt
Üblich sind auch „Tabulator-getrennte“ Dezimalzahlen
und - bisweilen – auch Dezimalzahlen mit ‚, , “ statt „ .“
(dann müssen für die Verwendung in python-Programmen
Dezimalkommata durch Dezimalpunkte ersetzt werden!)
Durchgesetzt haben sich „beschreibende“ Datenformate, z.B. xml = „extensible markup language“ oder
json = „java script object notation“ (≙ python dictionary);
gut durch pyhton-Module unterstützt !
Rein „binäre“ Datenformate (also sehr kompakte Dar- stellungen in Maschinen-abhängigem digitalem Format)
werden wegen ihrer Plattformabhängigkeit heute kaum noch verwendet.
# Daten im csv-Format lesen
# Datei zum Lesen öffnen
f = open('AudioData.csv', 'r')
# Kopzeile(n) lesen
header=f.readline()
print "Kopfzeile:", header
# Daten in 2D-numpy-array einlesen data = np.loadtxt(f,
delimiter=',', unpack=True)
print "-> Anzahl Spalten", data.shape[0]
print "-> Datenzeilen",
data.shape[1]
# Daten in 1D-arrays speichern
t = data[0]
a = data[1]
l = len(a)
Beispiel-Code
s. auch – PhyPraKit.readCSV()
– PhyPraKit.readtxt()
– PhyPraKit.readColumnData()
– PhyPraKit.labxParser()
u. a.
Bearbeitung/Analyse abgetasteter Daten
s. Funktionen in PhyPraKit und die Beispiele – test_convolutionPeakFinder()
– test_AutoCorrelation.py
– test_Fourier.py
– UnivariteSpline.py
u. a.
Frequenz eine periodischen Wellenform
Faltung mit verschobener Version von sich selbst:
Wegen der Selbstähnlichkeit der jeweiligen Perioden nimmt ρ Maxima an, wenn die Verschiebung einer Periodenlänge entspricht
i=1, …, n-1
Autokorrelation
# einfache Implementierung
...
l=len(a)
rho=np.zeros(l)
for i in range(1, l):
rho[i]=np.inner(a[:-i],a[i:])
rho[0] = np.inner(a, a)
rho = rho/rho[0]
...
animate_autoCorrelation.py
np.inner(a,b) kennen wir in der Mathematik als das Skalarprodukt von a und b
Frequenz eine periodischen Wellenform (2)
Bestimmung der Maxima und Minima
der Autokorrelation und Differenzbildung
liefert sehr genaue Bestimmung der
Periodendauer
Zeitdifferenz T aus Maxima: (2.0668 +/- 0 [0,0060]*)) ms
bzw. Minima: (2.0668 +/- 0.0049) ms der Autokorrelation
Anm.: Unsicherheiten aus Unsicherheit des Mittelwerts („ “ ); *) Die Unsicherheit von 0 auf T aus Maxima bedeutet, dass Abweichungen kleiner sind als die Sampling-Zeit ts = 0.0208 ms; als Unsicherheit wurde daher die Standardabweichung einer Rechteckverteilung mit Breite 0.0208 ms angenommen
Beispiel-Skript: test_AutoCorrelation.py
Parameterschätzung mit der Methode der kleinsten Quadrate („χ2-Methode“)
Minimiere
„Summe der Residuen-Quadrate“ S
bzgl. der k Parameter p
σi2 sind Varianzen der N Messungen yi
Falls die Fehler korreliert sind, ersetze 1/σi2 → V -1 (Inverse der Kovarianzmatrix)
folgt
bei gaußverteilten Fehlern σi
einer χ2-Verteilung mit
nf = N - k Freiheitsgraden
Erwartungswert: < Smin> = nf bzw. < Smin / nf > = 1
χ2-Wahrscheinlichkeit:
dient zur Quantifizierung der Qualität einer Anpassung
Aussage, mit welcher Wahrscheinlichkeit ein größerer Wert von S
am Minimum als der tatsächlich beobachtete zu erwarten wäre.
als notwendige Bedingung für ein Minimum
lösbar in Spezialfällen, z. B. für lineare Probleme
z. B. kafe:
https://github.com/dsavoiu/kafe
http://www.ekp.kit.edu/~quast/kafe/html
oder
Je schärfer das Minimum von χ2(a),
desto kleiner die Parameterfehler:
p
p
S(p)
scharfes Minimum:
große Krümmung
Parameterunsicherheiten
flaches Minimum:
kleine Krümmung
S(p)
Numerische Bestimmung des Minimums
siehe Übung 4.3:
– 11 Messpunkte Ai
– Modell „Resonanzkurve“, freier Parameter Dämpfung D
– Abstandsmaß
Die Funktion S(D) ist (nahezu) parabelförmig mit klarem Minimum bei D = 0.414
Die Unsicherheiten können ganz allgemein bestimmt werden, indem man
die Anpassung für viele Stichproben wiederholt und die Verteilungen der
Parameterwerte analysiert.
Ergebnisse der Anpassung eines
Polynoms 2. Grades an 10'000
den beobachteten Messwerten
entsprechende Stichproben
Auf diese Weise lassen sich - Parameterunsicherheiten
- Korrelationen der Parameter
bestimmen.
Das Verfahren ist
aufwändig
aber
exakt
Im Zweifelsfall wird die
Exaktheit statistischer
Verfahren mit solchen „Toy – Monte Carlos“
untersucht.
toyMCErros.py
Likelihood-Funktion L :
Produkt der Werte der Wahrscheinlichkeitsdichte Pi für n unabhängige Messungen xi :
hängt nur noch von den Parametern p ab !
Maximum-Likelihood-Prinzip:
Der beste Schätzwert für den Parametervektor ist derjenige, der die Likelihood-Funktion maximiert
negativer Logarithmus der Likelihood-Funktion:
Bedingung für Optimum:
Beispiel: Likelihood der Gaußverteilung
= χ2 !
const. bzgl. a
Für Gauß-förmig um f(xi , a) verteilte Messungen yi ist
χ2 äquivalent zu – 2 lnL
Ist f(xi , a) zusätzlich eine lineare Funktion der Parameter a = (ai ),
dann ist – 2 lnL (ai | yi) ( und χ2 (ai , yi) ) eine Parabel
Näherungsweise kann – lnL (ai | yi) in der
Nähe des Minimums häufig durch eine
Parabel approximiert werden !
Maximum-Likelihood: Prameterunsicherheiten
Mathematisch exakt: die angegebenen
Fehlerabschätzungen sind Untergrenzen
Nur für Parabel-förmigen
Verlauf von F(a) sind die
beiden Fehlerdefinitionen
äquivalent
Varianz ≈ 1 / Krümmung 1/σ2 ≈ ∂2F / ∂a2
bei mehreren Parametern ai: (cov-1)ij ≈ ∂2F / ∂ai ∂aj
±1σ - Intervall (=68%) aus ΔF = 0.5
Typischer Verlauf einer negativen log-Likelihood Funktion und ihrer 1. und 2. Ableitungen
F(a) näherungsweise quadratisch
um das Minimum;
1. Ableitung näherungsweise
linear, 0 am Minimum
2. Ableitung ~ konstant
Parabel aus Krümmung
am Minimum
±1σ
parameter a
Übersicht: Bestimmung der Parameter-Unsicherheiten
Lineare Probleme
mit gaußförmigen Unsicherheiten:
- analytische Lösung möglich (χ2 -Minimierung)
- Position des Minimums gegeben durch Linearkombination der Messwerte:
- Varianz der Parameterschätzung durch Fehlerfortpflanzung:
Nicht- lineare Probleme
oder andere als gauß- förmige Unsicherheiten:
- Likelihood-Analyse:
Ausnutzen der Cramer-Rao-Frechet-Grenze:
2. Ableitungen am Minimum:
Nicht-lineare Probleme
mit nicht-parabolischer
Likelihood am Minimum:
(für Grenzfall hoher Statistik)
- Scan der (Profil-) Likelihood in der Nähe des Minimums
Bei Unklarheit oder
sehr kleiner Statistik:
Monte-Carlo-Studie: Anpassung an viele der Verteilung der Daten entsprechende Stichproben, Verteilungen der Parameter auswerten.
Achtung: nur die letzten beiden Methoden liefern „Konfidenz-Intervalle“ (z. B. 1σ ≙ 68%)
zu kafe gibt es eine Anzahl von gut dokumentierten Beispielen, s.
https://github.com/dsavoiu/kafe
http://www.ekp.kit.edu/~quast/kafe/html
oder
Anmerkung:
die Berücksichtigung von Unsicherheiten in Ordinaten-
und Abszissenrichtung oder die Anpassung von Modellen,
die nicht linear von den Parametern abhängen, sind
analytisch nicht möglich.
Die Empfehlung ist daher, grundsätzlich Anpassungen mit numerischen Werkzeugen durchzuführen.
Software-Paket zur Modellanpassung
In der alltäglichen Praxis verwendet man heute Programme,
(bzw. besser Softwarepakete) zur Anpassung von Modellen an Parameter:
Geradenanpassung: verschiedene Parametrisierungen
Ergebnisse für das angepasste Modell
sind identisch,
aber im 1. Fall sind
die Parameter sehr
stark korreliert.
s. Übungsaufgabe 5.3
Parametrisierungen:
1.
2.
1.
2.
Daten erzeugt mit Script Generate_OhmLawExperiment.py
Auswertung: Ohm'sches Gesetz (2)
s. Übungsaufgabe 6
Ergebnis (Methode 2) :
Hilfsanpassung: Parametrisierung des
Temperaturverlaufs als
Funktion der Spannung
Auswertung: Ohm'sches Gesetz (3)
s. Übungsaufgabe 6
Ergebnis (Methode 3) :
Analytisches Modell
λw aus Hilfsanpassung
eingeschränkt
Bei Anwendung der χ2 Methode …
nie die Messwerte transformieren Unsicherheiten sind dann nicht mehr gaußförmig und damit die Voraussetzungen für die Methode nicht mehr gegeben
Also:
ggf. Unsicherheiten der Therorievorhersage quadratisch addieren
Also:
ggf. “Constraints“ und „Fixieren“ von Parametern benutzen, um externe Abhängigkeiten zu berücksichtigen
ggf. gemeinsame (systematische) Unsicherheiten mit Hilfe der Kovarianzmatrix berücksichtigen Also:
χ2 - Methode: Do's und Don'ts (2)
bei nicht-linearen Problemen unbedingt den Verlauf der χ2-Kurve am Minimum überprüfen
kafe: Methode plot_profile()
ggf. asymmetrischen Fehler angeben !
Konturen von Paaren von Parametern überprüfen
kafe: Methode plot_contour()
falls nicht elliptisch, Grafik veröffentlichen
immer die χ2-Wahrscheinlichkeit überprüfen kafe: Variable fit.chi2prob
bei Anpassung mehrerer Parameter
immer die Kovarianzen / Korrelationen überprüfen und angeben Kovarianzmatrix in kafe: Variable fit.par_cov_mat
bei großen Korrelationen ggf. Parametrisierung ändern, Korrelationen angeben (sonst keine Interpretation als Wahrscheinlichkeit möglich)
x=np.linspace(0.,10, 200)
y=np.sin(np.pi*x)/x
plt.plot(x,y,label='example function')
plt.legend()
plt.show()
i.w. das Ergebnis dieses
kurzen Programms:
siehe simplePlot.py
Darstellung einer
Funktion
Datenpunkte mit Fehlerbalken
und einer Geraden
ym=[2.31,1.27,4.40,4.98,5.40,6.98,7.53,8.54,9.65,10.15]
xm=range(1,11)
plt.errorbar(xm,ym,yerr=0.5,fmt='ro')
x=np.linspace(0.,11.,100)
y=x+0.5
plt.plot(x,y)
plt.show()
simpleData.py
Allgemeine Vorlage zur Funktkionsdarstellung
# -*- coding: utf-8 -*-
# diese Zeile legt die Codierung von Umlauten fest
########################################################
# script PlotBeispiel.py
'''
Beispiel fuer Funktionsdarstellung mit matplotlib
HILFE:
Ausfuehren dieses Programms: python PROGRAMMNAME
Dieses Programm ist in der Sprachversion 2 von python geschrieben.
numpy-Objekte werden mit prefix "np." aufgerufen
matplotlib-Objekte werden mit prefix "plt." aufgerufen
.. author:: Günter Quast <g.quast@kit.edu>
für den Kurs Computergestützte Datenauswertung
'''
#--------------------------------------------------------------
. . .
PlotBeispiel.py
noch eine (flexiblere) Vorlage
# -*- coding: utf-8 -*-
# diese Zeile legt die Codierung von Umlauten fest
###########################################
# script FunctionPlotter.py
'''
Funktionsdarstellung mit matplotlib
.. author:: Günter Quast <g.quast@kit.edu>
für den Kurs Computergestützte Datenauswertung
'''
#--------------------------------------------------------------
. . .
siehe FunctionPlotter.py
nutzt myPlotstyle.py
animate_wave.py
animated_Gauss.py
throwCoin.py
Beispiele unter http://matplotlib.org
Auswertung einer Längenmessung
Daten in Text-Editor eingeben:
# Messung 28. Apr. 2016
# h[cm]
xx.xx
yy.yy
….
liveData.dat
import numpy as np
import matplotlib.pyplot as plt
#read columns from file:
infile="Data.dat"
A=np.loadtxt(infile,unpack=True)
print A # check input
# visualize as histogram (with matplotlib)
fig=plt.figure(figsize=(10,10))
plt.hist(A,100)
plt.title('Distribution of Data')
print "Mittelwert: %.3f" % np.mean(A)
print "Standardabweichung: %.3f" % np.std(A)
print "Unsicherheit auf Mittelwert: %.3f" \
% (np.std(A)/np.sqrt(len(A)))
plt.show()
read_array.py
Mit kleinem Python-script darstellen und auswerten →
s. auch liveDisplay.py
(in der Vorlesung verwendet)
Wahrscheinlichkeit: Kopf oder Zahl ?
Beispiel Münzwurf:
Wahrscheinlichkeit für Kopf: pK=0.5
Wahrscheinlichkeit für Zahl: pZ=1-pK=0.5
Führe N=1, .., N Computer-Experimnte durch, berechne jeweils die Häufigkeit hn= NK(n) / n
i
# throw a coin N times
import numpy as np
N=500
f=[ ]
Nh=0
for n in range(N):
if np.random.rand()>0.5:
Nh+=1
f.append(float(Nh)/(n+1.))
throwCoin.py
Häufigkeit nähert sich der Wahrscheinlichkeit an:
hn → pK(n)
n
hn
Häufigkeitsverteilung: diskret
Beispiel Würfeln:
Häufigkeit der Augenzahlen
1, 2, …, 6 bei 30 Würfen
# example to produce random numbers
# as obtained by throwing dice
import numpy as np
ndice=1 # number of dice
nthrow=30 # how often to throw
numbers=[ ]
for i in range(nthrow):
r=0
for j in range(ndice):
r += int(6*np.random.rand())+1
numbers.append(r)
throwDice.py
Erwartete Verteilung
Häufigkeitsverteilung mit Verteilungsdichte
Gauss.py
Häufigkeitsverteilung (Histogramm) mit Verteilungsdichte (rote Kurve) und Beschriftung
Binomial ↔ Poisson ↔ Gauss - Verteilungen
import numpy as np, matplotlib.pyplot as plt, scipy.special as sp
def fGauss(x,mu=0.,sigma=1.): # Gauss distribution
return (np.exp(-(x-mu)**2/2/sigma**2)/np.sqrt(2*np.pi)/sigma)
def fPoisson(x,mu): # Poisson distribution
k=np.around(x)
return (mu**k)/np.exp(mu)/sp.gamma(k+1.)
def fBinomial(x,n,p): # Binomial distribution
k=np.around(x)
return sp.binom(n,k) *p**k *(1.-p)**(n-k)
#--------------------------------------------------------------
# plot basic distributions
mu, sig = 50., np.sqrt(mu)
x = np.arange(20, 80., 0.1)
plt.plot(x, fGauss(x,mu,sig), 'r-')
plt.plot(x, fPoisson(x,mu), 'b-')
p, n = 0.1, mu/p
plt.plot(x, fBinomial(x,n,p), 'g-')
# ( … ) # produce nice text
plt.show()
← einige spezielle Funktionen in scipy.special [ n! = gamma(n+1) ]
basicDistributions.py
python / numpy / matplotlib:
numpy.histogram(),
numpy.histogram2d()
matplotlib.pyplot.bar(),
matplotlib.pyplot.pcolormesh()
siehe auch:
matplotlib.pyplot.hist() matplotlib.pyplot.hist2d()
Beispielscript Histogram.py
hist2d statistics: <x> = 10.0093 <y> = 19.9881
var_x = 4
var_y = 9.1
cov = 4.5, cor = 0.75
Mittelwert von 10 Messungen yi mit Unsicherheiten σ entspricht der Anpassung einer konstanten Funktion f(x;c)=c
analytisch:
identisch zum „Mittelwert“
berechnen und grafisch darstellen
„numerisch“:
Script PlotAverage-withChi2.py
Beispielcode: numerische Anpassung mit kafe
In der Praxis setzt man Programmpakete ein, die - Strukturen zur Verwaltung von Daten und deren Fehlern - Definition von Modellen - Anpassung mittels numerischer Optimierung - grafische Darstellung - Ausgabe der Ergebnisse bereit stellen. Zuverlässige Konvergenz zum globalen Minimum bei nicht-linearen Problemstellungen erfordert meist das Setzen geeigneter Startwerte!
Beispiel mit dem
python-Framework
„kafe“:
Anpassung der drei Parameter einer qua-
dratischen Funktion an Datenpunkte mit Unsicherheiten nur
in Ordinatenrichtung.
siehe Script fitexample_kafe.py
68%-Konfidenzband der Modellanpassung
Beispielcode: numerische Anpassung (scipy)
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
# --------------------------------------------------------------------------------
def chi2_fit():
#-- set data points
xm = np.array([.05,0.36,0.68,0.80,1.09,1.46,1.71,1.83,2.44,2.09,3.72,4.36,4.60])
ym = np.array([0.35,0.26,0.52,0.44,0.48,0.55,0.66,0.48,0.75,0.70,0.75,0.80,0.90])
ye = np.array([0.06,0.07,0.05,0.05,0.07,0.07,0.09,0.1,0.11,0.1,0.11,0.12,0.1])
#-- least-squares fit with scipy.optimize.curve_fit
p, cov = curve_fit( poly2, xm, ym, sigma=ye, absolute_sigma=True )
print "Fit parameters:\n", par
print "Covariance matrix:\n", cov
#-- plot data and fit result
xp = np.linspace( 0., 5., 100 )
ffit = poly2(xp, p[0], p[1], p[2])
plt.errorbar(xm, ym, yerr=ye, fmt='o')
plt.plot( xp, ffit, '-' )
plt.xlim( 0, 5 )
plt.ylim( 0, 1 )
plt.show()
#-- define fit function
def poly2(x, a=1., b=0., c=0.):
return a * x**2 + b * x + c
if __name__ == '__main__': # –-----------
chi2_fit()
curvefit_example.py
Tool curve_fit aus scipy.optimize
zur χ2-Anpassung und numerischen Optimierung
Die Unsicherheiten können ganz allgemein bestimmt werden, indem man
die Anpassung für viele Stichproben wiederholt und die Verteilungen der
Parameterwerte analysiert.
Ergebnisse der Anpassung eines
Polynoms 2. Grades an 10'000
den beobachteten Messwerten
entsprechende Stichproben
Auf diese Weise lassen sich - Parameterunsicherheiten
- Korrelationen der Parameter
bestimmen.
Das Verfahren ist sehr
aufwändig, aber exakt !
Im Zweifelsfall wird die
Exaktheit statistischer
Verfahren mit solchen „Toy – Monte Carlos“
untersucht.
toyMCErros.py
Beispiel: Likelihood beim Münzwurf
Erinnerung: Binomialverteilung beim Wurf einer Münze:
relative Häufigkeit des Auftretens
des Ergebnisses „Kopf“
Für einige der Ergebnisse aus
der Reihe von Münzwürfen
oben ist nebenan die jeweilige
Likelihood-Funktion gezeigt:
Mit zunehmender Zahl an Würfen wird der Parameter p
durch die Likelihood-Funktion immer genauer eingegrenzt
ist eine Funktion
des Parameters p für gegebene
Beobachtung (N, k )
Skript nlLCoin.py
[** Beispiel: Likelihood für Poisson-Prozess ]
Stichprobe:
Poisson-verteilte Zufallszahlen
{ki }= { 0, 2 , 4 , 4 , 2 ,1, 1, 1, 1 }
Darstellungen von –ln L (λ | {ki , i ≤ n })
für die jeweils n ersten Elemente der Stichprobe
Ein typisches Beispiel für die Anwendung der negLogL-Methode sind kleine
Stichproben von Poission-verteilten Beobachtungen mit kleinem Erwartungswert:
allg. Implementierung ist
etwas anders, s. Codebeispiel
nlLExample.py
kafe example1 : Vergleich von zwei Modellen
# import everything we need from kafe
from kafe import *
# additionally, import the two model functions we want to fit:
from kafe.function_library import linear_2par, exp_2par
############
# Load the Dataset from the file
my_dataset = Dataset(input_file='dataset.dat',
title="Example Dataset")
### Create the Fits
my_fits = [ Fit(my_dataset, exp_2par),
Fit(my_dataset, linear_2par) ]
### Do the Fits
for fit in my_fits:
fit.do_fit()
### Create the plots, save and show output
my_plot = Plot(my_fits[0], my_fits[1])
# show data only once (it's the same data
my_plot.plot_all(show_data_for=0) )
my_plot.save('plot.pdf')
my_plot.show()
linear_2par
chi2prob 0.052
HYPTEST accepted (CL 5%)
exp_2par
chi2prob 0.96
HYPTEST accepted (CL 5%)
Die Methode kafe.file_tools.build_fit_from_file()
ermöglicht es, einfache Anpassungen
vollständig über eine Datei zu steuern:
import sys, matplotlib.pyplot as plt, kafe
from kafe.file_tools import buildFit_fromFile
# check for / read command line arguments
if len(sys.argv)==2:
fname=sys.argv[1]
else:
fname='data.fit'
print '*==* script ' + sys.argv[0]+ ' executing \n',\
' processing file ' + fname
# initialize fit object from file and run fit
theFit = buildFit_fromFile(fname)
theFit.do_fit()
thePlot = kafe.Plot(theFit)
thePlot.plot_all( show_info_for=None)
#thePlot.save(fname.split('.')[0]+'.pdf') #
#theFit.plot_correlations() # eventually contours
# show everything on screen
thePlot.show()
Skript kafe_fit-from-file.py
# example showing fit driven by input file
# ----------------------------------------
# this file is to be parsed with
# kafe.file_tools.buildFit_fromFile()
# Meta data for plotting
*BASENAME quadraticFitExample
*TITLE some test-data
*xLabel $x$
*xUnit
*yLabel $f(x)$
*yUnit
# the data
*xData
0.50 0.1
1.00 0.1
- - - -
9.50 0.1
10.00 0.1
*yData
0.45 0.1
0.82 0.1
- - - -
7.07 0.1
7.29 0.1
*FITLABEL quadratic function of x
*FitFunction
def fitf(x,a,b,c):
~~return a*x**2 + b*x +c
*InitialParameters # set initial values and range
0. 0.5
1. 0.5
0. 0.5
Datei data.fit
gleiches Skript, anderer Input: Mittellung korrelierter Messungen
Gleiches Skript (kafe_fit-from-file.py)
mit anderen Eingabedaten
(Mess-)Werte mit Unsicherheit und
Kovarianzmatrix (angegeben als
Wurzel aus den Nebendiagonalelementen)
aveCorrDat.fit
# Simulated correlated measurements
# - common error between pairs of measurements
# - common error of all measurements
# ----------------------------------------------------------------------
*BASENAME averageCorrDat
*TITLE simulated correlated measurements
*xLabel number of measurement
*yLabel values
#*xData # commented out, as not needed for averaging
*yData_SCOV
# val err syst sqrt. of cov. mat. elements
0.82 0.10 0.15
0.81 0.10 0.15 0.15
1.32 0.10 0.15 0. 0.
1.44 0.10 0.15 0. 0. 0.15
0.93 0.10 0.15 0. 0. 0. 0.
0.99 0.10 0.15 0. 0. 0. 0. 0.15
# common (correlated) error for all
*yAbsCor 0.05
# python code of fit function
*FITLABEL Average
*FitFunction
@ASCII(expression='average')
@LaTeX(name='f', parameter_names=('m',), expression='m')
@FitFunction
def fitf(x, m=1.): # fit an average
~~~~return m
*InitialParameters
1. 1.