Fakultät für PhysikInstitut für Experimentelle Kernphysik

www.kit.eduKIT – Die Forschungsuniversität in der Helmholtz-Gemeinschaft

Günter Quast

SS '19

Vorlesung:

Computergestützte Datenauswertung

Parameteranpassung in der Praxis – kafe

Beispiel kafe

Vom Prinzip zur (wissenschaftlichen) Praxis

Bisher nur einfachste Grundlagen besprochen: an festen Stützstellen xi ausgewertete, parameter-

behaftete Modellfunktion f(xi ; {p} ) zur Beschreibung von Daten yi mit statistischen Unsicherheiten ∆i

Dies beherrschen zahlreiche (kommerzielle und freie) Programme: Origin, QtiPlot, GnuPlot, scipy curve_fit, MS Excel, ...

in der Praxis ist es aber meist komplizierter: auch die Stützstellen xi haben Unsicherheiten ∆xi

Unsicherheiten ∆i (und auch ∆xi) evtl. korreliert

manche der Parameter {p} sind evtl. „eingeschränkt“ (als externe Messwerte oder durch „Hilfsmessungen“)

angepasste Parameter häufig korreliert, deren Unsicherheiten nicht gaußförmig

Parameteranpassung: Zusatzfunktionen

Bisher nur Unsicherheiten in Richtung der abhängigen Größe betrachtet.

Wie lassen sich zusätzlich Unsicherheiten in Richtung der unabhängigen Größe behandeln ?

Vergleich der Übereinstimmung des Modells mit den Daten – wir kennen die χ2-Wahrscheinlichkeit; Unsicherheit des angepassten Modells durch Parameterunsicherheiten

Es hat sich eingebürgert, die Modellunsicherheit zusätzlich zur angepassten Kurve als „Konfidenzband“ anzuzeigen.

nur ganz wenige der verfügbaren Anpassungspakete beherrschen solche Zusatzfunktionen !

Wie kann man (externe) Einschränkungen an die Modellparameter berücksichtigen?

Nichtlineare Anpassungen erfordern zusätzliche Überprüfungen: – ist die parabolische Näherung am Minimum gerechtfertigt ?

– sind die Korrelationen linear, d.h. ist die 1σ-Kontour elliptisch?

Behandlung von Fehlern in x- und y-Richtung

Implementierung: Methode Fit() der Root-Klasse TGraphErrors, oder kafe

1. Ableitung nutzen, um x-Fehlerin y-Fehler umzurechnen undQuadratisch zu den y-Fehlen addierenIteratives Verfahren:

1. Anpassung ohne x-Fehler

2. f '(xi) bilden und neue Fehler bestimmen:

3. Schritt analog 2 zur Kontrolle wiederholen;

χ2 darf sich nicht mehr stark ändern!

Geometrische Interpretation: Minimierung des auf projizierte Messfehler normierten Abstands d der Punkte von der Tangente an die Funktion:

Bisher nur Unsicherheiten in Richtung der abhängigen Größe betrachtet

Konfidenzband um das angepasste ModellEine Anpassung liefert Schätzwerte für die Parameter

sowie ihre Unsicherheiten in Form der Kovarianzmatrix

Mit Hilfe der linearen Fehlerfortpflanzung lässt sich die Unsicherheit auf das angepasste Modell bestimmen:

Übliche Darstellung: helles ±1σ-Band um die Zentralwerte des Modells

±1σ-Band der Modellanpassung

Beispiel kafe

Berücksichtigung von EinschränkungenBisweilen tauchen in Anpassungsproblemen Parameter auf, die auf externem Wissen bzw. externen Messungen beruhen: - Naturkonstanten - Materialeigenschaften - Parameter der Messapparatur aus dem Datenblatt - …

Es gibt mehrere Möglichkeiten zur Behandlung:

1. Parameter im Theorie-Modell als Konstante enthalten, für Parameter ohne bzw. mit vernachlässigbar kleiner Unsicherheit

2. gut bekannte Parameter, deren Einfluss studiert werden soll, sind im Theorie-Modell als freie Parameter enthalten (analog zu den anzupassenden Parametern), können in der Anpassung festgehalten, aber auf verschiedene Werte gesetzt werden (Bsp. kafe: Methode fix_parameters() )

3. nur ungenau bekannt Parameter sind im Theorie-Modell enthalten und werden auch in der Anpassung frei gelassen, aber als weiterer „einschränkender Messwert“ (engl. constraint) hinzugefügt ( kafe: constrain_parameters() )

also: für einen auf c = c0 ± σc eingeschränkten Parameter

Profil-Likelihood und (Konfidenz-) Konturen

Bei nichtlinearen Anpassungen ist nicht garantiert, dass die (üblicherweise angenommene) parabolische Näherung genau genug ist.

Liefert auch den Konfidenzbereichfür die Parameterwerte.

Beispiel einer nichtlinearen Anpassung (mit kafe )

39,3% CL ≙ 1σ (Δχ2=1)

Zur Überprüfung dient ein Scan der χ2 – Funktion um das Minimum für jeden einzelnen bzw. für Paare von Parametern. Dabei wird für jeden Punkt auf den Kurven das Minimum bzgl. aller anderen Parameter bestimmt.

(„Profil-Likelihood“).

Parameter-abhängige Unsicherheiten

Oft hängen die Unsicherheiten von den Messwerten ab,

z.B. relative Fehler, ,

oder Poisson-Fehler,

Dann resultiert häufig eine „verzerrte“ Anpassung

Abwärtsfluktuationen der Messung

→ zu kleine Unsicherheit

→ überhöhtes Gewicht in der Anpassung

bzw.

Aufwärtsfluktuationen der Messung

→ zu große Unsicherheit

→ zu kleines Gewicht in der Anpassung

Unsicherheiten aus der Modellfunktion bestimmen,

dazu Anpassung iterieren

1.) mit den aus der Beobachtung geschätzten Messungen anpassen

2.) verbesserte Schätzung der Unsicherheit aus den nun bekannten

Funktionswerten verwenden

[ 3.) zur Konvergenzkontrolle kann man ein drittes mal iterieren ]

Abhilfe:

Illustration mit 20% rel. Unsicherheit

Einfluss einzelner Fehlerkomponenten

Häufige Frage:

Wie groß ist der Einfluss einer bestimmten Fehlerkomponente, k ,

auf die Unsicherheit , für einen Parameter a ?

Methodik zur Beantwortung:

Außer der Anpassung mit allen n Fehlerkomponenten für die Messwerte i =1,N eine weitere ausführen, bei der eine Komponente (m) weggelassen wird:

Die auf die Komponente m zurückzuführende Unsicherheit ergibt sich dann als

Diese Beziehung gilt, wenn die Fehlerkomponenten voneinander unabhängig sind !

Achtung: wenn sich der Zentralwert in der zweiten Anpassung signifikant verschiebt, Fehlerkomponente nur skalieren statt weglassen !

Parameterfehler ohne Datenfehler ?

Manche Programme zur Anpassung (QTIplot, Origin usw.) geben bei χ2-Anpassungen Parameterfehler aus, ohne dass Fehler auf die Datenpunkte angegeben wurden. Wie geht das?

Annahmen:

Modell beschreibt die Daten perfekt: χ2 := nf statt < χ2 >= nf

Alle Datenpunkte haben den gleichen Fehler: σi :=σ (nur selten so !)

d.h. alle

dann setzen und Parameterfehler bestimmen

– Verlust der χ2-Wahrscheinlichkeit – angenommenes Fehlermodell oft falsch !

Beispiele mit kafe

Software-Paket zur Modellanpassung

In der alltäglichen Praxis verwendet man heute Programme, (bzw. besser Softwarepakete) zur Anpassung von Modellen an Parameter:

➢ Einlesen der Eingabedaten, evtl. Unterstützung verschiedener Formate

➢ Verwaltung der Unsicherheiten, evtl. für Abszisse und Ordinate, evtl. Kovarianzmatrizen bzw. Korrelationen

➢ kafe: Berücksichtigung von Einschränkungen an Modellparameter

➢ Bereitstellung verschiedener Fit-Funktionen, evtl. frei programmierbar

➢ Durchführung der Anpassung mittels numerischer Optimierung, evtl. anpassbar

➢ Ausgabe des Ergebnisses, als Textdatei und Grafik

➢ Bewertung der Qualität der Anpassung (außer χ2 bisweilen auch andere gebräuchlich)

➢ evtl. Ausgabe der Korrelationen bzw. Kovarianz-Matrix der Parameter

➢ kafe: graphische Darstellung der Modellunsicherheit

➢ kafe: asymmetrische Fehler mittels Scan der χ2 -Kurve am Minimum

➢ kafe: Ausgabe von Kovarianz-Ellipsen (d. h. Wahrscheinlichkeitskonturen)

Struktur von kafekafe = Karlsruhe fit environment

– implementiert als pyhton-Klassen: Dataset, Fit und Plot sind die wichtigsten

– Fit-Funktion als python Code – zahlreiche Hilfsfunktionen zur Dateneingabe– ansprechende (konfigurierbare) grafische Ausgabe

kafe: Beispielezu kafe gibt es eine Anzahl von gut dokumentierten Beispielen, s.

https://github.com/dsavoiu/kafehttp://www.ekp.kit.edu/~quast/kafe/html oder

kafe example0 ''' general example for fitting with kafe - construct dataset - specify the error model using the add_error_source() method - perform fit (2nd order polynomial) - show and save output'''import kafefrom kafe.function_library import quadratic_3par

#### create a Dataset instance:my_dataset = kafe.Dataset( data = ([0.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], [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]), title='some data', axis_labels=['$x$', '$y=f(x)$'])

#### specify the error modelmy_dataset.add_error_source('y', 'simple', [0.06,0.07,0.05,0.05,0.07,0.07,0.09,0.10,0.11,0.10,0.11,0.12,0.10])

#### Create the Fit objectmy_fit = kafe.Fit(my_dataset, quadratic_3par)# Set initial values and error estimatesmy_fit.set_parameters((0., 1., 0.2), (0.5, 0.5, 0.5))# Do the Fitmy_fit.do_fit()

#### Create result plots and output themmy_plot = kafe.Plot(my_fit)my_plot.plot_all()my_plot.show()

Auszug aus example0_intro.py

Ausgabe(n) in Verzeichnis .kafe/ <name>.log Ergebnis fit.log fit log file minuit.log minuit log file

# 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()

kafe example1 : Vergleich von zwei Modellen

linear_2parchi2prob 0.052HYPTEST accepted (CL 5%)exp_2parchi2prob 0.96HYPTEST accepted (CL 5%)

kafe example1 : Vergleich von zwei Modellen# check contourscontour1 = my_fits[0].plot_contour(0, 1, dchi2=[1., 2.3])profile00=my_fits[0].plot_profile(0)profile01=my_fits[0].plot_profile(1)contour2 = my_fits[1].plot_contour(0, 1, dchi2=[1., 2.3])

#contour1.savefig('kafe_example1_contour1.pdf')#contour2.savefig('kafe_example1_contour2.pdf')#profile00.savefig('kafe_example1_profile00.pdf')#profile01.savefig('kafe_example1_profile01.pdf')

Geradenanpassung: verschiedene Parametrisierungen

Ergebnisse für das angepasste Modell sind identisch, aber im 1. Fall sinddie Parameter sehrstark korreliert.

s. Übungsaufgabe 5.3

Parametrisierungen:

1.

2.

1.

2.

Daten-getriebene Anpassung

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, kafefrom kafe.file_tools import buildFit_fromFile

# check for / read command line argumentsif 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 fittheFit = 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 screenthePlot.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 *FitFunctiondef fitf(x,a,b,c):~~return a*x**2 + b*x +c

*InitialParameters # set initial values and range0. 0.51. 0.5 0. 0.5

Datei data.fit

Anpassung (hier Mittelung) korrelierter Messwerte

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')@FitFunctiondef fitf(x, m=1.): # fit an average~~~~return m

*InitialParameters1. 1.

Mittelung korrelierter Messwerte (2)Bei den Daten handelt es sich um unser Beispiel aus Vorlesung 4:

6 Studenten in 3 Gruppen messen mit jeweils eigenem Messgerät vom gleichen Typ, von allen angewandte „Theorie-Korrektur“ mit Unsicherheit

Δf_i = ±0.10 unabhängiger Einzelmessfehler der 6 Einzelmessungen

Δs = ±0.15 korreliert zwischen den 3 Gruppen Δt = ±0.05 allen gemeinsamer „Theoriefehler“

Den Wert der Unsicherheit auf denMittelwert kann man einfach aus der Zahl der unabhängigen Messungen für jede Komponente der Unsicherheit verstehen:

Ergebnis der χ2 Anpassung mit Kovarianzmatrix:

m = 1.052 ± 0.108

Kovarianzmatrix V besteht aus 3 Blöcken :

6 x 0.10 → 0.10 / √6 3 x 0.15 → 0.15 / √3 1 x 0.05 → 0.05 ---------------------Wurzel aus der quadr. Summe: 0.108

Kovarianzmatrix tut genau, was man intuitiv erwartet !

Last words …

Das war die letzte Vorlesung aus der Reihe „Einführung in die Computergestützte Datenanalyse“

es war viel Material, konzipiert als „Learning by Doing“ → Beispiele aus Vorlesung und Übungen wichtig hoffentlich hilfreich bei zukünftigen Aufgaben (Praktika, Bachelorarbeit, ...)

die Vorlesungsfolien – z.T. mit „optionalem“ mathematischem Hintergrund, den ich in der Veranstaltung selbst nicht angesprochen habe – bleiben als Referenz auf dem Web-Server ( www.ekp.kit.edu/~quast )

→ kehren Sei bei Bedarf zu diesen Folien zurück

Wenn Sie es noch nicht getan haben: bearbeiten Sie die Übungsaufgaben ! „Programmieren lernt man durch Programmieren“

Bleiben Sie kreativ und neugierig !

Datenanalyse xxL

In dieser einführenden Veranstaltung konnte bei weitem nicht alles und schon gar nicht in der angebrachten Tiefe angesprochen werden; es gibt noch viele weitere Fragestellungen …

(einige) Antworten in „Rechnernutzung in der Physik“ im 5. Semester, oder „Moderne Methoden der Datenanalyse“ im Master

Modellanpassung mit Nebenbedingungen: - Parameter innerhalb von Grenzen, a < λi < b - Einschränkungen an Parameter, z.B. durch Funktionen ej( λ1, …, λp) = 0

Wie funktioniert numerische Optimierung / Minimierung ?

Konfidenzintervalle / Grenzen für Parameter (statt „Punktschätzung“)

Hypothesentests und Klassifizierung

Methoden des „Machine Learning“

… (und weitere Themen)