Vorlesung:

Rechnernutzung in der Physik Anpassung von Parametern

Günter Quast Fakultät für Physik Institut für Experimentelle Kernphysik

KIT – Die Forschungsuniversität in der Helmholtz-Gemeinschaft

WS 2016/17

www.kit.edu

Parameteranpassung

χ2 und Likelihood Parameterschätung (= Funktionsanpassung) mit der Methode der kleinsten Fehlerquadrate (χ2-Methode) Maximum Likelihood – Methode Anpassungen mit kafe Fertige Programmpakete

Parameterschätzung

(neudeutsch „fitten“)

Anpassung von Modellen = parameterabhängige an statistische Daten = Messwerte

Funktionen

⃗y

Messdaten ( xi , yi ) mit angepasster Funktion f(x;p)

p ist der Vektor der Parameter der Funktion

Benötigen ein „Abstandsmaß“ für Übereinstimmung von (xi, yi) und f(xi) Mathematisch: yi und fi = f(xi) sind Elemente eines Vektorraums: Vektoren y und f 2

Abstand = Betrag des Differenzvektors, d = (y also: Finde Parameterwerte, die d2 minimieren

- f) ∙ (y – f)

(Skalarprodukt)

(für geeignet definiertes Skalarprodukt)

Frage: wie werden Messfehler berücksichtigt?

Beziehung

Messung ↔ Funktionswert

Wahrscheinlichkeitsverteilung um den wahren Wert wahrer Wert

Messung mit Fehlerbalken

übliche Darstellung

= Messpunkt

Fehlerbalken entspricht ±1σ dieser Gaußkurve

eigentliche Bedeutung

Messung mit Fehlerbalken bedeutet: Beobachtung eines Messergebnisses, das der Summe aus einem „wahren“ Wert und einer Zufallszahl aus einer Verteilungsdichte entspricht

Parameterschätzung

2

Methode der kleinsten Quadrate (χ - Methode)

5

Parameterschätzung mit der Methode der kleinsten Quadrate („χ2-Methode“) Entwickelt zu Beginn des 19. Jahrhunderts von Legendre, Gauß und Laplace. Ziel: finden der besten Funktion, die Fehler-behaftete Datenpunkte annähert. zunächst für Geraden (Ausgleichsgerade oder lineare Regression), aber mit numerischen Minimierungsalgorithmen (z.B. MINUIT) auf beliebige Funktionen und Anzahl Parameter anwendbar !

Minimiere bzgl. der Parameter {p} S: „Summe der Residuen-Quadrate“

= f (x; {p})

σi2 sind die Varianzen der N Messunge Falls die Fehler korreliert sind, ersetze 1/σi2 → cov -1 (Inverse der Kovarianzmatrix)

Methode ist auch aus „Likelihood-Prinzip“ herleitbar

→ später

Minimieren von S N Messungen k Parameter

Wie minimiert man

●

analytisch:

als notwendige Bedingung für ein Minimum

lösbar in Spezialfällen, z. B. für lineare Probleme: ●

i. a. numerisch „numerische Optimierung: Algorithmen zur Suche nach dem (eine Minimum einer Skalaren Funktion im k-dimensionalen Paramete In der Praxis werden heute auch für lineare Probleme die (effizienten) numerischen Minimierungsmethoden verwendet. (außer in Spezialfällen, z. B. bei zeitkritischen oder immer wieder vorkommenden Problemstellungen)

Beispiel mit einem Parameter Mittelwert von 10 Messungen yi mit Unsicherheiten σ entspricht der Anpassung einer konstanten Funktion f(x;c)=c

analytisch:

identisch zum „Mittelwert“

„numerisch“:

berechnen und grafisch darstellen

Script PlotAverage-withChi2.py

Mittelung bei unterschiedlichen Unsicherheiten

N Messungen yi mit individuellen Unsicherheiten σi

„1/σ - quadrat gewichtete Summe „

Wichtiges Ergebnis zur Mittelwertbildung bei Messungen mit unterschiedlichen Unsicherheiten

Zusammenhang von S mit der χ2 - Verteilung Quadrat einer standard-normal verteilten Größe; folgt χ2Verteilung Smin, die „gewichtete Summe der Residuen-Quadrate“ am Minimum, folgt bei Gauß-verteilten Fehlern σi einer

χ2-Verteilung mit nf = N-k Freiheitsgraden. Erwartungswert: =nf oder = 1 Die χ2-Wahrscheinlichkeit

dient zur Quantifizierung der Qualität einer Anpassung Aussage, mit welcher Wahrscheinlichkeit ein größerer Wert von χ2 am Minimum als der tatsächlich beobachtete zu erwarten wäre.

Güte der Anpassung χ2-Wahrscheinlichkeit

in python: scipy.stat.chi2.cdf(S, nf)

E[Smin] = nf → E[Smin / nf] = 1 Die Wahrscheinlichkeit, einen χ2 - Wert größer nf zu erhalten, ist ca. 40% (ziemlich unabhängig von n) Mit wachsender Zahl der Freiheitsgrade sollte Smin / nf immer näher bei 1 liegen

in python: chi2prb= 1. - scipy.stat.chi2.cdf(S, nf)

Das (allg.) Lineare Problem lineares Problem: anzupassende Funktion ist linear in den Parametern

Beispiele:

Allgemein lassen sich also die N Messwerte yi durch eine Linearkombination von Funktionswerten von k Funktionen an Stützstellen xi darstellen:

Dies lässt sich sehr kompakt in Matrix-Scheibweise darstellen:

mit Kovarianzmatrix V der N Messungen ergibt sich für S(aj ):

dann: Minimieren von S(aj ) bzgl. der Parameter aj , d.h.

Lösung des allg. linearen Problems Notwendigen Bedingung für ein Minimum von S(aj ):

Für lineare Probleme analytisch lösbar. Anm.: häufig führt aber eine Erweiterung des anzupassenden Modells dazu, dass Parameter in einer der nicht-linearen Funktionen auftreten, dann muss auf andere (=numerische) Methoden zurück gegriffen werden. [ s. nächste Folien ]

Lösung

ist lineare Funktion der Messwerte:

Kovarianzmatrix der Parameter durch Fehlerfortpflanzung:

Analytische Lösung des

allg. linearen Problems

Kleine (Matrix-) Rechnung: Produktregel, ∇a auf 1. und 2. Term mit a anwenden

nach dem Vektor der Schätzwerte auflösen

Die Lösung für die Parameter ist also eine Linearkombination der Messwerte, mit der „Lösungsmatrix“ B:

Die Kovarianz-Matrix ist damit durch das Fehlerfortpflanzungsgesetz gegeben: Einsetzten der expliziten Form von B und ausnutzen, dass V und ATA symmetrische Matrizen sind

Spezialfall:

Lineare Regression (Geradenanpassung)

Durch Einsetzen in die Lösungsformeln von vorhin erhält man mit den Abkürzungen

die Lösung:

Diese Formeln waren für Generationen von Studierenden die Basis einer jeden Anpassung („Regression“).

Anmerkung:

Lineare Regression

Das Kovarianzelement V12 verschwindet für Sx=0, d. h. wenn der Erwartungswert der Abszissenwerte 0 ist

Dies lässt sich durch Änderung der Parametrisierung erreichen: setze Man erhält dann diese einfacheren, unkorrelierten Lösungen für die Parameter

Bestimmung der Parameterfehler χ2

χ2 2

Je schärfer das Minimum von χ (a), desto kleiner die Parameterfehler:

a

a

scharfes Mimimum: große Krümmung

→

flaches Mimimum: kleine Krümmung

Parameterunsicherheiten sind umgekehrt proportional zu Krümmung(en) von χ2(a) am Minimum

bzw. bei mehreren Parametern .

Erinnerung:

lineares Problem

Für das lineare Problem mit ergab sich die Kovarianzmatrix der Parameter durch Fehlerfortpflanzung

Bilden der zweiten Ableitungen von S, liefert das gleiche Ergebnis:

Beispiel: Fehler des Mittelwerts Erinnerung: Mittelwertbestimmung

von 10 Messungen yi mit Unsicherheiten σ

Unsicherheit nach Formel:

Identisch zu früherem Ergebnis (das wir mittels Fehlerfortpflanzung erhalten hatten)

Beispiel: Qualität der Anpassung Erinnerung Mittelwertbestimmung:

Numerische Lösung Script PlotAverage-withChi2.py

Qualität der Anpassung: χ2 = 6.23 bei 10 Messwerten und 1 Parameter, d. h. nf = 9 Freiheitsgrade

→ χ2 - Wahrscheinlichkeit = 0.72 bestimmt mit der kumulativen Verteilung der χ2-Verteilung chi2prb= 1. - scipy.stat.chi2.cdf(6.2, 9)

Numerische Verfahren

Numerische Anpassung I. a. (heute) Anpassungen mittels numerischer Optimierung durch: Gittermethode: ●

z.B. Minuit: SCAN-Algorithmus

k gleichverteilte Testwerte pro Dimension ➔

Erfordert kd Berechnungen bei d Dimensionen Ungeeignet für große d

Monte-Carlo-Methode: ●

Funktionsberechnung an zufällig verteilten Testpunkten ➔

Auch bei großen d geeignet gut für Bestimmung von Startwerten

Einfache Parametervariation: ●

Eindimensionale Minimierung in einem Parameter, dann Minimierung in nächstem Parameter → Iteration ➔

I.a. schnelle Konvergenz, wenn Minimierung in Richtung der Eigenvektoren der Matrix der 2. Ableitung erfolgt

weitere numerische (deterministische) Verfahren: ●

●

Simplex-Verfahren (ohne Ableitungen) Methoden mit angepasster Schrittweite, numerische Auswertung der Ableitungen

z.B. Minuit: SIMPLEX-Algorithmus MIGRAD-Algorithmus

Beispiel

Simplex-Verfahren

Illustration der Simplex-Methode am Beispiel der Rosenbrock'schen „Bananenfunktion“ Simplex: n-dimensionaler Polyeder aus n+1 Punkten im Rn Bei jedem Schritt wird der schlechteste Punkt xr mit größtem Funktionswert F(xr) durch einen neuen, besseren, ersetzt, oder der Simplex ggf. verkleinert oder vergrößert.

[** Simplex-Verfahren ] n+1 Punkte x1 , ..., xn + 1 im Rn n-dimensionaler Polyeder oder „Simplex“

●

Sortierung, so dass F(x1 ) < … < F(xn + 1 )

●

Schwerpunkt der n besten Punkte: c = ∑ i = 1 … n xi / n

●

Spiegelung des schlechtesten Punktes an c: xr = c +  (c – xn + 1 ) Falls F(x1 ) < F(xr ) < F(xn ) → xr ersetzt xn + 1 Falls F(xr ) < F(x1 ) gute Richtung → Streckung: xs =c+(xr – c),  > 1 Falls F(xs ) < F(xr ): xs ersetzt xn+1 , ansonsten xr ersetzt xn + 1 Falls F(xr ) > F(xn ): Simplex zu groß → Abflachung: xs = c –  (c – xn + 1 ), 0 <  < 1 Falls F(xs ) < F(xn + 1 ): xs ersetzt xn + 1

●

Ansonsten Kontraktion um x1 : xj = x1 +  (xj – x1 ), 0 <  < 1

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 quadratischen Funktion an Datenpunkte mit Unsicherheiten nur in Ordinatenrichtung. siehe Script fitexample_kafe.py

68%-Konfidenzband der Modellanpassung

Behandlung von Fehlern in x- und y-Richtung 1. Ableitung nutzen, um x-Fehler in y-Fehler umzurechnen und Quadratisch zu den y-Fehlen addieren Iteratives Verfahren: 1. Anpassung ohne x-Fehler 2. f '(xi) bilden und neue Fehler bestimmen: 3. Schritt analog 2 zur Kontrolle wiederholen; χ2 am Minimum darf sich nicht stark ändern !

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

Implementeirung: Methode Fit() der Root-Klasse TGraphErrors

Parameteranpassung mit ROOT Klasse TGraphErrors void TGraphFit() { //Draw a graph with error bars and fit a function to it. //set global options gStyle->SetOptFit(111); //superimpose fit results // make nice Canvas TCanvas *c1 = new TCanvas("c1","Daten",200,10,700,500); c1->SetGrid(); //define some data points const Int_t n = 10; Float_t x[n] = {-0.22, 0.1, 0.25, 0.35, 0.5, 0.61, 0.7, 0.85, 0.89, 1.1}; Float_t y[n] = {0.7, 2.9, 5.6, 7.4, 9., 9.6, 8.7, 6.3, 4.5, 1.1}; Float_t ey[n] = {.8,.7,.6,.5,.4,.4,.5,.6,.7,.8}; Float_t ex[n] = {.05,.1,.07,.07,.04,.05,.06,.07,.08,.05}; // copy data to TGraphErros object TGraphErrors *gr = new TGraphErrors(n,x,y,ex,ey); gr->SetTitle("TGraphErrors mit Fit"); gr->Draw("AP"); // now perform a fit(with errors in x and y!) gr->Fit("pol3"); c1->Update();

Script TGraphFit.C

Anpassung mittels Root Macro: > root TGraphFit.C

Behandlung von Fehlern in x- und y-Richtung (2) Allgemein bei korrelierten Fehlern: Kovarianzmatrizen Cx und Cy im zweiten Schritt x-Fehler zur Kovarianzmatrix addieren:

Geschrieben mit dem Residuenvektor und dem Vektor der 1. Ableitungen, ergibt sich

Implementiert in kafe

kafe: Beispiele zu kafe gibt es eine Anzahl von gut dokumentierten Beispielen, s. http://www.ekp.kit.edu/~quast/kafe/html

oder

https://github.com/dsavoiu/kafe

Hinweis:

kafe update

Es gibt eine neue Version von kafe mit einigen (meist kosmetischen) Änderungen; !!! setzt matplotlib ≥1.5.0 voraus Installation:

pip install -- upgrade kafe

bzw. sudo pip install --upgrade kafe

(Linux systemweit)

Eine beliebige kafe-Version im Nuztzer-Bereich installiert man mittels pip install --user kafe==1.x.x

(x.x: 0.1 oder 1.0)

(kann auf den CIP-Pool Rechnern notwendig sein)

Wesentliche Neuerung: Fit() enthält eine Option: fit_name='name' name wird für Log-Files und Grafiken verwendet, erleichtert die Buchführung bei vielen Fits im gleichen Skript

Beispielcode:

numerisch Anpassung mit kafe siehe Script fitexample_kafe.py

# example fit with kafe import kafe from kafe.function_tools import FitFunction, LaTeX, ASCII

(

# fit function definition (with decorators for nice output) @ASCII(expression='a * x^2 + b * x + c') @LaTeX(name='f', parameter_names=('a','b','c'), expression=r'a\,x^2+b\,x+c') @FitFunction def poly2(x, a=1.0, b=0.0, c=0.0): return a * x**2 + b * x + c # --------- begin of workflow ---------------------# set data xm = [.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 = [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 = [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] # create a kafe Dataset kdata = kafe.Dataset(data=(xm, ym), basename='kData', title='example data') kdata.add_error_source('y', 'simple', ye) # add uncertainties kfit=kafe.Fit(kdata, poly2) # create the fit from data & fit function kfit.do_fit() # perform fit kplot=kafe.Plot(kfit) # create plot object kplot.axis_labels = [r'$St\"utzstellen \, x $', r'$data\,\&\,f(x)$'] kplot.plot_all() # make plots kplot.show() # show the plots

kafe nutzt das CERN-Paket MINUIT zur numerischen Optimierung

) nur zur Verschönerung der wichtige Code für kafe

Beispielcode:

numerische Anpassung (scipy )

Tool curve_fit aus scipy.optimize import matplotlib.pyplot as plt import numpy as np zur χ2-Anpassung und numerischen Optimierung 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__': # –----------curvefit_example.py chi2_fit()

Unsicherheiten mit Toy-MC Die Unsicherheiten können ganz allgemein bestimmt werden, indem man die Anpassung für viele Stichproben wiederholt und die Verteilungen der Parameterwerte analysiert.

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.

Ergebnisse der Anpassung eines Polynoms 2. Grades an 10'000 den beobachteten Messwerten entsprechende Stichproben toyMCErros.py

Maximum Likelihood Methode

Kleinste Quadrate & Maximum Likelihood Zwei Möglichkeiten, die am besten zur Messung passende Verteilung zu finden:

kleinster Abstand Messung ↔ Erwartungswert

Minimiere Abstand vom Sollwert

Messung Maximale Wahrscheinlichkeit

Maximiere Höhe der pdf

Likelihood-Methode

Beide Methoden bevorzugen in diesem Beispiel die durchgezogene Verteilung

Likelihood-Methode Mehrere unabhängige Messgrößen xi, i =1, …, n einer Größe folgen einer Verteilungsdichte p(x,a) mit Parametern ai. Als Verteilungsdichte ist p(x;a) positiv und normiert (bzgl. x, aber nicht bzgl. a !)

Beispiel von 7 Messungen und Likelihood bzgl. zweier Verteilungen Skript likelihood-pdf.py

„Likelihood“ ist das Produkt der Wahrscheinlichkeiten p(xi |a) aller Messungen

Anmerkung Likelihood-Methode Mit Hilfe der Likelihood-Methode können neben der Parameterabhängigkeit auch verschiedene Verteilungen verglichen werden:

Wichtig für „Hypothesentests“

Maximum Likelihood-Prinzip Likelihood-Funktion: Produkt der Werte der Wahrscheinlichkeitsdichte, Pi , für n unabhängige Messungen xi :

hängt nur noch von den Parametern a ab ! Maximum-Likelihood-Prinzip: Der beste Schätzwert für den Parametervektor, der die Likelihood-Funktion maximiert

, ist derjenige,

Anmerkungen zur Likelihood-Methode Maximum Likelihood (ML) Eingeführt von Sir Ronald A. Fisher 1912 zur „bestmöglichen“ Schätzung von Parametern von Verteilungen aus Stichproben engl. „likelihood“ = Wahrscheinlichkeit (ähnlich „probability“); im Deutschen wird der Begriff „Likelihood“ übernommen Basis der mathematischen Statistik (→ Fisher-Information) sehr universell einsetzbare Methode ( ältere Methode der kleinsten Quadrate ist Spezialfall) Likelihood ist eine zentrale Größe in der klassischen (frequentistischen) und in der Bayes'schen Statistik erlaubt die Schätzung von Konfidenzintervallen für Parameterwerte Für gegebene Stichprobenwerte xi aus einer zu Grunde liegenden Verteilungsdichte f ist die Likelihood-Funktion eine Funktion der Modellparameter: (L „lebt“ im Parameterraum)

Maximum Likelihood in der Praxis Technische und theoretische Gründe: Minimiere den negativen Logarithmus der Likelihood-Funktion:

Likelihood-Gleichung definiert Schätzwert(e)

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 – lnL Ist f(xi , a) zusätzlich eine lineare Funktion der Parameter a = (ai ), dann ist – lnL (ai | yi) ( und χ2 (ai , yi) ) eine Parabel – lnL (ai | yi) kann in der Nähe des

Minimums häufig durch eine Parabel approximiert werden !

Parameterunsicherheiten mit Hilfe der Likelihood-Funktion

Maximum-Likelihood:

Prameterunsicherheiten

Parabel aus Krümmung am Minimum

F(a) näherungsweise quadratisch um das Minimum;

± 1σ - Intervall (=68%) aus ΔF = 0.5

±1σ

1. Ableitung näherungsweise linear, =0 am Minimum

2. Ableitung ~ konstant Varianz ≈ 1 / Krümmung 1/σ2 ≈ ∂2F / ∂a2 bei mehreren Parametern ai: (cov-1)ij ≈ ∂2F / ∂ai ∂aj Typischer Verlauf einer negativen log-Likelihood Funktion und ihrer 1. und 2. Ableitungen

Mathematisch exakt: die angegebenen Fehlerabschätzungen sind Untergrenzen

Nur für Parabel-förmigen Verlauf von F(a) sind die beiden Fehlerdefinitionen äquivalent

Prameterunsicherheiten (2) Plausibilitätserklärung (kein Beweis, Stichwort „Cramer-Rao-Frechet Grenze“) nur ein Parameter a, betrachten Taylor-Entwicklung von F(a) um Minimum:

näherungsweise parabelförmig Likelihood = exp(-F(a)) als Verteilungsdichte in a auffassen:

ist Gauß-Verteilung (mit Normierungsfaktor A)

Standardabweichung gegeben durch

weiter gilt mit dieser Beziehung für σ:

Einschub:

Parabel-Eigenschaften

allg. Darstellung einer Parabel

Wenn F(a) eine negative Log-Likelihood Funktion ist, dann ist

Prameterunsicherheiten (3) zur Fehlerbestimmung aus F(a):

Log-Likelihood- Differenz bestimmt Fehler. Vorteil dieser Methode: invariant unter Variablentransformation a → a'(a): F(a) – F(â) = F( a'(a) ) – F( a'(â) )

Das so bestimme Unsicherheitsintervall entspricht 1σ einer Gaußverteilung (68% Konfidenz-Intervall)

Fehlerbestimmung: |Δ (-ln L )| Δχ2 _____________________

1σ 2σ 3σ nσ

| | | |

0.5 2.0 4.5 n2/2

| | | |

Wichtig, wenn –ln L nicht parabelförmig in der Nähe des Minimums: Angabe eines asymmetrischen Fehlerintervalls

1 4 9 n2

Übersicht: Methoden zur Bestimmung der 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:

Nicht-lineare Probleme mit nicht-parabolischer Likelihood am Minimum:

- Scan der (Profil-) Likelihood in der Nähe des Minimums

2. Ableitungen am Minimum:

(für Grenzfall hoher Statistik)

Bei Unklarheit oder sehr kleiner Statistik:

Monte-Carlo-Studie: Anpassung an viele der Genauigkeit und Verteilung der Daten entsprechende Stichproben, Verteilungen der Parameter studieren.

Achtung: nur die letzten beiden Methoden liefern „Konfidenz-Intervalle“ (z. B. 1σ ≙ 68%)

Beispiel:

Likelihood für Poisson-Prozess

Ein typisches Beispiel für die Anwendung der negLogL-Methode sind kleine Stichproben von Poission-verteilten Beobachtungen mit kleinem Erwartungswert:

Stichprobe: Poisson-verteilte Zufallszahlen

{ki }= { 0, 2 , 4 , 4 , 2 ,1, 1, 1, 1 }

allg. Implementierung ist etwas anders, s. Codebeispiel nlLExample.py

Darstellungen von

–ln L (λ | {ki , i ≤ n })

für die jeweils n ersten Elemente der Stichprobe

Beispiel: Likelihood für Poisson-Prozess (2) Im Grenzfall großer Ereigniszahlen ist -lnL von vorne herein parabelförmig (in solch einem Fall wäre auch eine χ2 -Anpassung gerechtfertigt)

2. Stichprobe: Poisson-verteilte Zufallszahlen

{ki}= {41,48,62,49,58,45,38,52,46}

Darstellungen von

–ln L (λ | {ki , i ≤ n })

für die jeweils n ersten Elemente der Stichprobe

Zusammenhang -ln L und χ2 Für Gauß-förmig um f(xi; a) verteilte Messungen yi ist die χ2 Methode äquivalent zur -lnL-Methode :

χ2

Fehlerbestimmung: |Δ (-ln L)| Δχ2 const. bzgl. a 2

Minimieren von -ln L ↔ Minimieren von χ ∆(-ln L) = ½ ∆χ2 ∂2(-ln L) / ∂ai∂aj = ½ ∂χ2 / ∂ai∂aj

____________________________

1σ 2σ 3σ nσ

| | | |

0.5 2.0 4.5 n2/2

| | | |

1 4 9 n2

Bei anderen als Gauß-förmigen Fehlerverteilungen ist χ2 eine eigenständige Methode; - bei unbekannter Fehlerverteilung haben wir keine bessere - χ2 ist optimal für die Anpassung von Linearkombinationen von Fit-Funktionen

Maximum Likelihood vs. Kleinste Quadrate Vergleich: Voraussetzung

Maximum - Likelihood Kleinste Quadrate Mittelwert und Varianz PDF exakt bekannt bekannt

Methode

Höhe der PDF

Abweichung vom Mittelwert der PDF

Effizienz

maximal

maximal bei linearen Problemen

Komplexität Robustheit korrelierte Datenfehler Güte der Anpassung

aufwändig, meist nichtoft linear und exakt lösbar linear nein - PDF muss exakt bekannt sein nein („Ausreißer“) u.U. kompliziert einfach über Kovarianzmatrix nein

ja: χ2-Wahrscheinlichkeit 53

Parameterschätzung

Parameteranpassung mit Root

54

Parameteranpassung mit ROOT ROOT enthält einige Minimierungsalgorithmen, u.a. das aus der FORTRAN-Zeit stammende (und nach C++ umgeschriebene) MINUIT, entstanden am CERN - gut getestet und anerkannt

Standard in der Teilchenphysik

Auswahl von Minimierern über die Klasse TvirtualFitter (s. später) In vielen Fällen reicht die von ROOT per GUI zur Verfügung gestellte Funktionalität: • -2 logL (!) und χ2-Anpassungen an Histogramme • χ2-Anpassungen in Klasse TGraphErrors • Vordefinierte Funktionen: Polynome bis zum 9. Grad, Gauss, Exponential- und Landauverteilung χ2-Wert FCN=82.016 FROM MIGRAD STATUS=CONVERGED 12 CALLS1 3 TOTAL EDM=9.82506e-16 STRATEGY= 1 ERROR MATRIX ACCURATE EXT PARAMETER STEP FIRST NO. NAME VALUE ERROR SIZE DERIVATIVE 1 p0 1.54084e+02 3.51831e+00 1.56525e-02 1.25994e-08

Text-Ausgabe von MINUIT

55

Anpassen mit ROOT GUI Root-Klassen TH1, TH2 und TH3 sowie TGraph und TgraphErrors enthalten eine Methode .Fit() zur Funktionsanpassung (TF1 bzw.TF2 u.TF3) In der graphischen Darstellung kann durch Rechtsklick das FitPanel aktiviert werden, um Anpassungen interaktiv vorzunehmen. root[0]TGraphErrors ∗gr=new TGraphErrors("ExampleData.txt"); root[1]gr−>Draw("AP"); // Achsen und Daten root[2]gr->Fit("pol1"); // Gerade

Ein einfaches Beispiel

56

Parameteranpassung mit ROOT Man kann auch eigene Fit-Funktionen aus dem Repertoire von ROOT definieren: Root > TF1 *myfit = new TF1("myfit","[0]*sin(x) + [1]*exp(-[2]*x)", 0, 2); // set parameters Root > myfit->SetParameters(1,0.05,0.2); // Fitten Root > hist->Fit("myfit");

… oder beliebige Funktionen selbst schreiben: Double_t myfunc (Double_t *x, Double_t *par) { // IHR CODE !!! } Root > TF1 *myf = new TF1("myf",myfunc, , , ); // set parameters Root > myf->SetParameters(, … , hist->Fit("myf");

Im allgemeinsten Fall kann man auch die χ2- oder -2lnL- Funktion selbst vorgeben: Die Methode von TVirtualFitter void SetFCN(void (*fcn)(Int_t &, Double_t *, Double_t &f, Double_t *, Int_t))

erlaubt das Setzen der zu minimierenden Funktion (mit Namen „fcn“) void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)

f ist dabei der Rückgabewert, im Array par übergibt Root die Parmeter 57

Parameteranpassung mit ROOT

Ein Code-Fragment zur Berechnung von χ2 mit Kovarianz-Matrix

// fcn function for fit with covariance matrix double arr_poi; // global array pointer TMatrix *icovp; // global matrix pointer int nval = 3; // global for number of measurements int main() { arr_poi = new double[3]; TMatrix cov(3,3); cov(0,0) = ... , cov(2,2) = ...; // Invert matrix TMatrix icov = cov.Invert(); icovp = &icov; ...; } void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) { ... for (int i; i < nval; i++ ) { for (int j; j < nval; j++ ) { chi += (arr_poi[i] - fit_func(i)) * (arr_poi[j] - fit_func(j)) * (*icovp)(i,j) ; } } ... } 58

Anpassen von Funktionen an Histogramme

b: Binbreite 1. χ2 -Methode: Probleme: * Bins mit Ni=0, d.h Fehler 0 Abhilfe: setze * Ni folgen Binomial- oder Poisson-Verteilung, χ2-Methode nimmt Gauß-Fehler an: Methode nur gut für „große“ Ni 59

Anpassen an Histogramme: „Binned Likelihood-Fit“ 2. Binned Likelihood-Fit

hängt nicht von der Funktion ab

⇒

:

Schätzung für Parameter a

χ2 – und Likilood-Anpassung für Histogramme in ROOT implementiert ( χ2 ist default! ) Achtung: ROOT verwendet -2 ln L 60

Anpassen mit ROOT – Beispiel Histogramm-Fit

Vorsicht b. kleinen Zahlen!

10'000 Histogramme dieser Art, jeweils Gauß(x;μ,σ) anpassen mit der χ2-Methode Test: Pull p =

Beispiel: pull.c

(fitted mean – μ)

/error

muss standard-normalverteilt sein ist es aber in diesem Fall nicht ! Grund: angenommene Fehler in jedem Bin sind √ni - falls ni nach unten fluktuiert, wird auch der Fehler kleiner angenommen, die Folge ist einer verzerrte Parameterschätzung (betrifft hier vor allem den Parameter σ !)

Abhilfe: -Log L – Anpassung: Mean 0.0020 RMS 1.007 61

Anpassung von Parametern: abschließende Bemerkungen Hier konnte nicht alles angesprochen werden; es gibt noch viele weitere Fragestellungen … Anpassen mit Nebenbedingungen: - Parameter innerhalb von Grenzen, a< λi