Concluding remarks .edu

Jupyter notebook tutorial: Advanced Fitting

G?nter Quast, April 2021

Jupyter Notebook Fundamentals

This file of type .ipynb contains a tutorial as a Jupyter notebook . Jupyter provides a browser interface with a (simple) development environment for Python code and explanatory texts in intuitive Markdown format. The input of formulas in LaTeX format is also supported.

A summary of the most important commands for using Jupyter as a working environment can be found in the notebook JupyterCheatsheet.ipynb (German). Basics for statistical data analysis can be found in the notebooks IntroStatistik.ipynb (German) and Fehlerrechnung.ipynb (German).

In Jupyter, code and text are entered in individual cells. Active cells are indicated by a blue bar in the margin. They can be in two states: in edit mode the input field is white, in command mode it is grayed out. Clicking in the border area selects the command mode, clicking in the text field of a code cell switches to edit mode. The esc key can also be used to leave the edit mode.

Pressing a in command mode creates a new empty cell above the active cell, b creates one below. Entering dd deletes the corresponding cell.

Cells can be either of the type Markdown or Code . Entering m in command mode sets the type Markdown, entering y selects the type Code.

The cell content is processed - i.e. setting text or executing code - by entering shift+return , or alt+return if a new, empty cell should also be created.

The settings mentioned here as well as the insertion, deletion or execution of cells can also be executed via the pull-down menu at the top.

Section 1: Introduction

Fitting data to models is a general task in physics. Typically, parameter-dependent models represent theoretical knowledge, which is confronted with empirical data, i.e. measurements. Measured data are affected by unavoidable experimental uncertainties, and the question then is whether there exists a set of model parameters for which the agreement with the data is acceptable. Speaking in term of statistics, the step of comparing data with uncertainties to a parameter-dependent model is a hypothesis test. Usually, a quantity called goodness-of-fit is used to quantify the level of agreement between the data and the theoretical model, expressed as a p-value, that is the probability to observe a goodness-of-fit which is worse than the one actually observed.

Recap: least-squares method

An

example

is

the

value

of

the

sum

of

least

squares,

S,

at

the

values

of

the

optimal

parameters,

known

as

2

.

The

corresponding p-value is the so-called 2-probability

p2 =

, 2

(x, ndof )dx

Smin

where S ist given by the deviation of the model function f(x, p), depending on parameter values p evaluated at positions xi,

from the measurements yi normalised to the measurement uncertainties i,

N

2

f (xi ; p ) - yi

S = (

) (Equ. 1.1)

i=1

i

S is minimized with respect to the parameter values p, resulting in an optimal parameter set p and the value S_min

0

(xi ,

yi ,

i ,

p

)

0

If p2 is below a pre-defined level, the hypothesis that the model describes the data is discarded. If the model is not rejected, the

parameter set p is accepted as the best-fit set of model parameters. The procedure skeched here is known as the least0

squares method originally proposed by Carl Friedrich Gau?. It is the basis of many fits in fundamental laboratory courses, and

many implementations of this simple method exist in a variety of tools and libraries.

The least-squares method can easily extended to cope with correlated uncertainties by introducing the covariance matrix V of the measurement uncertainties:

T

S (y , V f (x, p )) = (y - f (x, p ))

-1

V (y

- f (x, p ))

(Equ. 1.2)

The covariance matrix of the parameter uncertainties can be determined by evaluating the second derivatives of S at the minimum:

1

-1

(Vp^ ) =

ij

2

2

S(p )

(Equ. 1.3)

pi pj

p^i p^j

Many of the wide-spread tools provide an implementaion, e.g. the easy-to-use function curve-fit in the scipy.optimize package.

1.1 More complex use cases

However, there are more complex use cases. Uncertainties on the abscissa values in addition to uncertainties on the ordinate can be implemented by transforming variations on the abscissa to variations on the ordinate by means of a first-order Taylor expansion. There are not many tools providing such functionality, among them the package ODR (orthogonal distance regression) in scipy , the Method TGraphErrors in the CERN root package. Correlated uncertainties on the abscissa are, however, not treated by these packages.

Another complication consists in relative uncertainties, which are quite common in physics. A prominent example are calibration uncertainties of measurement devices, which lead to correlated, relative uncertainties of all measurements taken with the same device. Another common example are approximations of Poisson uncertainties by a Gaussian distribution by taking the square root of the mean as the uncertainty, \sigma = \sqrt{\mu}. Such relative uncertainties should be related to the true and not to the measured values in order to avoid biases. In doing so, the standard version of the least-squares method is no longer strictly valid, as the normalisation term of the Gaussian distribution becomes dependent on the parameters and can no longer be neglected. In this case, a fit based on the full likelihood of all involved uncertainties must be used, as will be explained below.

To the author's knowledge, presently only the tool kafe2 developed at KIT provides the means to simultaneously handle all of the use cases and error types mentioned above. The package PhyPraKit.phyFit, also provides a light-weight, transparent layer implementing the required precedures on top of the function minimization and uncertainty analyis tool iminuit. Many of the examples described below will therefore use the kafe2 or phyFit packages.

1.2 Prerequisites and other tutorials

The following material assumes familiarity with the basic principles of statistical data analysis and applications of the leastsquares method, as explained in other jupyter tutorial of this series:

1. Introduction to Jupyter notebooks: JupyterCheatsheet.ipynb A. Introduction to Python : PythonIntro.ipynb B. Cheat-sheet for Python : PythonCheatsheet.ipynb C. Basics of matplotlib : matplotlibTutorial.ipynb D. Basics of Statistics: IntroStatistik.ipynb E. Error analysis in laboratory courses: Fehlerrechnung.ipynb F. Introduction to kafe2: kafe2Tutorial.ipynb

1.3 General settings for kafe2

In [1]: from __future__ import division, print_function # python2-compatibility import sys, os

# Lines with % or %% at the beginning are so-called "magic commands", # that specify the cell type or options for displaying graphics. %matplotlib inline

Imports and Presets for kafe2:

In [2]: from kafe2 import config, XYContainer, Fit, Plot

import numpy as np, matplotlib.pyplot as plt

1.4 A simple fit example

To start, a simple, very general example may serve as an introduction to the usage of kafe2. The following, well-documented piece of code is a complete python script performing a fit of a quadratic function to data with errors in y-direction.

Run the following cell by activating it via clicking, then type shift+return :

In [3]:

''' general example for fitting with kafe2 - set-up arrays for data - perform fit (2nd order polynomial) - show and save output

''' # Imports # from kafe2 import XYContainer, Fit, Plot import numpy as np, matplotlib.pyplot as plt

### define the model function def poly2(x, a=1.0, b=0.0, c=0.0):

return a * x**2 + b * x + c

### Workflow # # 1. the data x, y, e = ( [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], [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] )

# 2. convert to kafe2 data structure and add uncertainties

xy_data = XYContainer(x, y)

xy_data.add_error('y', e)

# independent erros y

# -- set meaningful names

xy_data.label = 'Beispieldaten'

xy_data.axis_labels = ['x', 'data & f(x)']

# 3. create the Fit object my_fit = Fit(xy_data, poly2) # set meaningful names for model my_fit.model_label = 'Parabel-Modell'

# 4. perform the fit my_fit.do_fit()

# 5. report fit results my_fit.report()

# 6. create and draw plots my_plot = Plot(my_fit) my_plot.plot()

# 7. show or save plots # ## plt.savefig('kafe_fitexample.pdf') plt.show()

######## # Data # ########

X Data ====== 0.05 0.36 0.68 0.8 1.09 1.46 1.71 1.83 2.44 2.09 3.72 4.36 4.6

Y Data ====== 0.35 0.26 0.52 0.44 0.48 0.55 0.66 0.48 0.75 0.7 0.75 0.8 0.9

Y Data Error ============ 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

Y Data Correlation Matrix ======================================== [1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1.]

######### # Model # #########

Model Function ==============

poly2(x; a, b, c)

X Model ======= 0.05 0.36 0.68 0.8 1.09 1.46 1.71 1.83 2.44 2.09 3.72 4.36 4.6

Y Model ======= 0.3251 0.3829 0.4392 0.4595 0.5065 0.5624 0.5976 0.6138 0.6886 0.6472 0.8057 0.8439 0.8548

############### # Fit Results # ###############

Model Parameters ================

a = -0.017 +/- 0.013 b = 0.193 +/- 0.059 c = 0.315 +/- 0.046

Model Parameter Correlations ============================

a

b

c

======= ======= ======

a 1.0

-0.9559 0.742

b -0.9559 1.0

-0.855

c 0.742 -0.855 1.0

Cost Function =============

Cost function: chi-square (with covariance matrix) chi2 / ndf = 9.646 / 10 = 0.9646 chi2 probability = 0.472

The default output shows the data, the best-fit function and the uncertainty band of the model in graphical form, and text output

with

names,

values

and

uncertainties

of

the

best-fit

values

of

the

parameters

as

well

as

the

value

of

2

per

degree

of

freedom

and the 2-probability of the fit. Customisations of the output, like colors, marker type and size, axis labels are easily possible -

see the kafe2 documentation and the tutorial kafe2Tutorial.ipynb.

1.5 A more complex kafe2 example: comparing two models

One of the strengths of kafe2 is its ability to handle various kinds of errors and multiple fits to the same data set, which are ideally suited to compare different model assumptions. An example, taken from the kafe2 jupyter tutorial, is shown in the code cell below. Compared to the previous example, errors in the x-direction are added, and two models are defined and fit to the data.

In [4]: # Comparison of two models with kafe2 # -> insert code here

# Our first model is a simple linear function: def linear_model(x, a, b):

return a * x + b

# Our second model is a simple exponential function. # The kwargs in the function header specify parameter defaults. def exponential_model(x, A0=1., x0=5.):

return A0 * np.exp(x/x0)

# The data for this exercise: x = [19.8, 3.0, 5.1, 16.1, 8.2, 11.7, 6.2, 10.1] y = [23.2, 3.2, 4.5, 19.9, 7.1, 12.5, 4.5, 7.2] data2 = XYContainer(x_data=x, y_data=y) data2.add_error(axis='x', err_val=0.3) data2.add_error(axis='y', err_val=0.15, relative=True)

data2.label = 'my data points' data2.axis_labels=['my x values', 'my y values']

# Create 2 Fit objects with the same data but with different model functions: linear_fit = Fit(data2, model_function=linear_model) exponential_fit = Fit(data2, model_function=exponential_model)

# Optional: Assign LaTeX strings to parameters and model functions. linear_fit.assign_parameter_latex_names(a='a', b='b') linear_fit.assign_model_function_latex_expression("{a}{x} + {b}") linear_fit.model_label = 'my linear model' exponential_fit.assign_parameter_latex_names(A0='A_0', x0='x_0') exponential_fit.assign_model_function_latex_expression("{A0} e^{{{x}/{x0}}}") exponential_fit.model_label = 'my exponential model'

# Perform the fits: linear_fit.do_fit() exponential_fit.do_fit()

# Optional: Print out a report on the result of each fit. linear_fit.report() exponential_fit.report()

# Optional: Create a plot of the fit results using Plot. p = Plot(fit_objects=[linear_fit, exponential_fit], separate_figures=False) p.plot(fit_info=True)

# Show the fit results: plt.show()

######## # Data # ########

X Data ====== 19.8 3.0 5.1 16.1 8.2 11.7 6.2 10.1

X Data Error ============ 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3

X Data Correlation Matrix ========================= [1. 0. 0. 0. 0. 0. 0. 0.] [0. 1. 0. 0. 0. 0. 0. 0.] [0. 0. 1. 0. 0. 0. 0. 0.] [0. 0. 0. 1. 0. 0. 0. 0.] [0. 0. 0. 0. 1. 0. 0. 0.] [0. 0. 0. 0. 0. 1. 0. 0.] [0. 0. 0. 0. 0. 0. 1. 0.] [0. 0. 0. 0. 0. 0. 0. 1.]

Y Data ====== 23.2 3.2 4.5 19.9 7.1 12.5 4.5 7.2

Y Data Error ============ 3.48 0.48 0.675 2.985 1.065 1.875 0.675 1.08

Y Data Correlation Matrix ========================= [1. 0. 0. 0. 0. 0. 0. 0.] [0. 1. 0. 0. 0. 0. 0. 0.] [0. 0. 1. 0. 0. 0. 0. 0.] [0. 0. 0. 1. 0. 0. 0. 0.] [0. 0. 0. 0. 1. 0. 0. 0.] [0. 0. 0. 0. 0. 1. 0. 0.] [0. 0. 0. 0. 0. 0. 1. 0.] [0. 0. 0. 0. 0. 0. 0. 1.]

######### # Model # #########

Model Function ==============

linear_model(x; a, b)

X Model ======= 19.8 3.0 5.1 16.1 8.2 11.7 6.2 10.1

Y Model ======= 18.33 2.501 4.481 14.85 7.402 10.7 5.517 9.193

############### # Fit Results # ###############

Model Parameters ================

a = 0.94 +/- 0.12 b = -0.33 +/- 0.75

Model Parameter Correlations ============================

a

b

======= =======

a 1.0

-0.8964

b -0.8964 1.0

Cost Function =============

Cost function: chi-square (with covariance matrix)

chi2 / ndf = 12.45 / 6 = 2.075

chi2 probability = 0.0527

######## # Data # ########

X Data ====== 19.8 3.0 5.1 16.1 8.2 11.7 6.2 10.1

X Data Error ============ 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3

X Data Correlation Matrix ========================= [1. 0. 0. 0. 0. 0. 0. 0.] [0. 1. 0. 0. 0. 0. 0. 0.] [0. 0. 1. 0. 0. 0. 0. 0.] [0. 0. 0. 1. 0. 0. 0. 0.] [0. 0. 0. 0. 1. 0. 0. 0.] [0. 0. 0. 0. 0. 1. 0. 0.] [0. 0. 0. 0. 0. 0. 1. 0.] [0. 0. 0. 0. 0. 0. 0. 1.]

Y Data ====== 23.2 3.2 4.5 19.9 7.1 12.5 4.5 7.2

Y Data Error ============ 3.48 0.48 0.675 2.985 1.065 1.875 0.675 1.08

Y Data Correlation Matrix ========================= [1. 0. 0. 0. 0. 0. 0. 0.] [0. 1. 0. 0. 0. 0. 0. 0.] [0. 0. 1. 0. 0. 0. 0. 0.] [0. 0. 0. 1. 0. 0. 0. 0.] [0. 0. 0. 0. 1. 0. 0. 0.] [0. 0. 0. 0. 0. 1. 0. 0.] [0. 0. 0. 0. 0. 0. 1. 0.] [0. 0. 0. 0. 0. 0. 0. 1.]

######### # Model # #########

Model Function ==============

exponential_model(x; A0, x0)

X Model ======= 19.8 3.0 5.1 16.1 8.2 11.7 6.2 10.1

Y Model ======= 26.4 3.349 4.335 16.75 6.345 9.755 4.962 8.014

############### # Fit Results # ###############

Model Parameters ================

A0 = 2.32 +/- 0.26 x0 = 8.14 +/- 0.64

Model Parameter Correlations ============================

A0

x0

====== ======

A0 1.0 0.8678

x0 0.8678 1.0

Cost Function =============

Cost function: chi-square (with covariance matrix)

chi2 / ndf = 5.508 / 6 = 0.9179

chi2 probability = 0.481

Examining the graphical output, ist is clear that the linear model only marginally describes the data. The inspection of the 2probability shows a value of only 5.3%, slightly above the usually chosen critical value of 5% for rejection of a model. On the other hand, the exponential model fits perfectly.

Proposed Exercise: You might try to add a quadratic model - you will see that this would also perfectly match the data.

1.6 Examining non-liner models

A further, very useful feature of kafe2 consists in producing profile scans of S around the minimum and confidence contours for pairs of parameters. Our second model is non-linear in the sense that one of the parameters, x0, enters in a non-linear way. Such problems typically lead to a non-parabolic behaviour of S as a function of the parameters, and to deviations of the contours corresponding to the one- and two- sigma levels from the elliptic form.

Execute the code in the cell below and study the output.

In [5]: # Checking the nonlinear fits # -> enter code here

from kafe2 import ContoursProfiler

# Create contour plot and profiles for the exponential fit cpf = ContoursProfiler(exponential_fit) cpf.plot_profiles_contours_matrix(show_grid_for='contours') plt.show()

In this case, the profiles of 2 = S(x, y, p ) - S(x, y, p0 are approximated acceptably well up to 2=1, but show

stronger deviations for larger values. In practice, this means that the values of the parameter uncertainties are trustworthy, but should not be extrapolated to larger values - i.e. two-times the given uncertainty does not correspond to the two-sigma level. In case of stronger deviations another feature of kafe2 should be used to determine confidence regions for the values of the fit parameters that are asymmetric around the best-fit values:

exponential_fit.report(asymmetric_parameter_errors=True)

The method used to derive such confidence regions relies on the correspondence between S and the negative logarithm of the likelihood function for Gaussian uncertainties on the data points (see explanation below).

Section 2: Least Squares and Likelihood methods

The least-squares method, often also called "2-method", is directly linked to the very general likelihood principle. Like S in the

section above, the likelihood L is a function of the data, x, y, and a model function f (x; p ) characterised by a set of parameters p, i.e. L (x, y , f(x, p )). If the measurements are independent, the likelihood is constructed from the product of the values of the probability density function D(y , f (x, p )) describing the fluctuations of the data y around the model prediction f (x, p). Very often, owing to the Central Limit Theorem of probability theory, D is the Gaussian distribution. In the most general

case, the multivariate Gaussian (or normal) distribution in d dimensions is the appropriate one to use:

D = N (y , V , f (x, p )) =

1

1

exp(- (y

-

T

) V

-1

(y

-

))

with

mean

values

(2)d det(V )

2

=

f

(x,

p).

(Equ.

2.1)

Instead of the product of likelihood factors, one uses the negative sum of the logarithm of the likelihood factors, and the "negative log-likelihood" or "nlL" for Gaussian distributed data becomes:

$ -2\, \ln{\cal{L_{\rm Gau?}}}\left( \vec y, V, \vec{f}(\vec x, \vec {p})\right) \,=\, \left(\vec y - \vec f(\vec x; \vec p ) \right)^T V^{-1} \left(\vec y - \vec f(\vec x; \vec p ) \right)

\ln(\det(V)) + n\,\ln(2\pi) \,. $ (Equ. 2.2)

Minimisation of this function w.r.t. p results in the best estimate of the parameter values according to the likelihood principle. It

can be shown that this estimation is asymptotically optimal and unbiased, i.e. the best among all possible estimators in the sense that is leads to the smallest variance and vanishing bias of the estimated parameter values in the limit of infinite sample sizes. These desirable properties already apply for small samples for the exponential family of distributions. The covariance matrix of the parameters for an optimal and unbiased likelihood-estimator is given by

2

-1

Vij

=

ln L()

pi dpj

where p^ indicates the best-fit value of the parameter p (Equ. 2.3).

p^ p^

ij

When minimizing nlL, constant terms can be omitted. If the covariance matrix is independent of the parameters, i.e. in the

absence of relative errors or of errors in the x-direction, -2 ln L (y , V , f(x, p )) becomes equal to S(y , V f(x, p), the

squared sum of redsiduals in Equ. 1.1. In other cases, the full likelihood as given in Equ. 2.2 must be used.

The evaluation of the determinant in every step of an iterative numerical optimization is, however, computationally expensive. On modern computers and thanks to very efficient numerical libraries this is not a problem at all any more, and therefore still existing simplified, less optimal methods should be abandoned in favour of the correct likelihood treatment. Usage of the full nlL is the default in the tools kafe2 and in phyFit. These packages effectively fall back to the traditional 2-method in cases where it is equivalent, or on explicit user request.

Note

that

sometimes

-2

ln

L

("n2lL")

is

used

instead

of

nlL

to

keep

the

similarity

to

the

2

method.

The

goodness-of-fit

is

still

correctly

quantified

by

2

(S(y , V

, f (x, p 0)),

. ndof )

This

can

be

justified

from

the

likelihood

principle: the best-possible model would describe the data perfectly, i.e. all residuals are zero in this case, and the smallest possible value of n2lL is ln det(V ). The difference between the observed best-fit value and the best possible value is therfore exactly 2(Smin, ndof ). As we will see below, such differences in n2lL are used to determine confidence intervals, and even in the general case 2 corresponds to such a nlL difference and therefore is a good measure of the goodness-of-fit that can be translated to a p-value.

2.1 Gaussian likelihood in kafe2

As long as the uncertainties are independent of the parameter values, the nlL and 2 methods are equivalent, and therefore fits with kafe2 yield results that are - apart from numerical precision - identical to those obtained with other tools. The example in the code cell below illustrates the effect by re-using the data of the previous example, but this time adding a relative uncertainty. The two fit results refer to a fit where the relative uncertainties are defined w.r.t. the measured values, or in the other case to the model values. To make the difference in the uncertainties visible, the x-data are shifted by 0.02 in the second fit. The example is a slightly modified version of the one contained in the kafe2 jupyter tutorial.

In [6]: # same linear fit as above, but with relative errors w.r.t. the model values.

x_1 = np.array(x)-0.02 data2_1 = XYContainer(x_data=x_1, y_data=y) data2_1.add_error(axis='x', err_val=0.3) ## relative y-errors w.r.t. model are added below, via method of fit class

data2_1.label = 'same data, rel. uncertainties w.r.t. model'

# Create Fit: linear_fit2_1 = Fit(data2_1, model_function=linear_model) # --> relative uncertainties with reference to model specified here: linear_fit2_1.add_error(axis='y', err_val=0.15, relative=True, reference='model')

# Optional: Assign LaTeX strings to parameters and model functions. linear_fit2_1.assign_parameter_latex_names(x='x', a='a', b='b') linear_fit2_1.assign_model_function_latex_expression("{a}{x} + {b}")

# Optional: Assign LaTeX strings to parameters and model functions. linear_fit2_1.assign_parameter_latex_names(a='a', b='b') linear_fit2_1.assign_model_function_latex_expression("{a}{x} + {b}") linear_fit2_1.model_label = 'same model, rel. uncertainties w.r.t. model'

# Perform the fit: linear_fit2_1.do_fit()

# Optional: print report. #linear_fit2_1.report(asymmetric_parameter_errors=True)

# Optional: Create a plot of the fit results using Plot. p2_1 = Plot([linear_fit, linear_fit2_1]) # Assign colors to data ... p2_1.customize('data', 'marker', [(0, 'o'), (1,'o')]) p2_1.customize('data', 'markersize', [(0, 5), (1, 5)]) p2_1.customize('data', 'color', [(0, 'grey'), (1,'red')]) # note: although 2nd label is suppressed p2_1.customize('data', 'ecolor', [(0, 'grey'), (1, 'red')]) # note: although 2nd label is suppressed # ... and model. p2_1.customize('model_line', 'color', [(0, 'mistyrose'),(1, 'orange')]) p2_1.customize('model_error_band', 'label', [(0, r'$\pm 1 \sigma$'),(1, r'$\pm 1 \sigma$')]) p2_1.customize('model_error_band', 'color', [(0, 'mistyrose'),(1, 'orange')])

p2_1.plot(asymmetric_parameter_errors=True)

# Create contour plot and profiles for the linear fit: cpf2_1 = ContoursProfiler(linear_fit2_1) cpf2_1.plot_profiles_contours_matrix(show_grid_for='contours')

# Show the fit results. plt.show()

There is a significant shift compared to the parameter uncertainties when taking relative uncertainties w.r.t. the model values. The reason is obvious, because clearly visibly different uncertainties are assigned to the data in the two cases. The profile-likelihood and the contour show that model-dependent uncertainties introduce some non-parabolic behaviour; even the fit of a linear function is no longer a "linear" problem and usually numerical methods are needed to find the optimum of nlL.

2.2 Confidence regions for parameters

All fitting tools implement Equ. 1.3 to determine the covariance matrix describing the parameter uncertainties. It should be noted that correlations of parameter uncertainties are quite common and hard to avoid - and therefore the covariance or correlation matrix of the fit results must be examined and quoted as part of the results if correlations exceed ~10%. However, in non-linear cases the description of the conficence regions of the parameter values by the covariance matrix alone is not sufficient. An example was shown above for the exponential model - the scan of the likelihood contour deviated from a parabola, and the contour plot revealed a non-elliptical behaviour. In such cases Equ. 1.3 does not provide the a fully satisfying descritption of the confidence regions of the best-fit parameter values. In such cases, more advanced methods must be applied.

The challenge is to estimate valid confidence regions for the fit parameters, taking into account parameter inter-dependencies, i.e. general correlations. Reporting confidence regions instead of point estimates of parameters with symmetric errors has now become standard in science. However, there are not many tools implementing this feature.

One way out is to use simulated data with fluctuations as given by the uncertainty model; this requires running a large number of fits and hence is computationally expensive, but exact. Another approach consists in marginalisation of the multi-dimensional likelihood function by integrating out all parameters except one and analysing its marginal distribution. This method also is computationally expensive.

A more favourable method relies on the so-called profile likelihood, which is determined for each parameter in turn by varying it around the best-fit value and minimizing the likelihood w.r.t. all other parameters. The resulting curves are exactly the profile curves that where already shown in the kafe2 examples above. A change in nlL by + 1 corresponds to a range of ?1 around

2

the mean of a Gaussian distribution, i.e. a confidence region of 68.3%. Since the change in the value of the likelihood is invariant under parameter transformations, this procedure also works for non-parabolic (i.e. non-Gaussian) cases. The same method works, with a significantly larger numerical effort, for contours in two dimensions, which thus become "confidence contours".

An efficient implementation of the profile-liklelihood method exists in the minimizer and uncertainty-examination tool MINUIT developed at CERN, which is used by phyFit and is the default option in kafe2. In phyFit, confidence ranges determined by profiling of the likelihood are the default, in kafe2 a special parameter is needed in the report function:

Fit.report(asymmetric_parameter_errors=True)

Example: non-linear fit

The code cell below show an extremely non-linear example. In such a case, asymmetric uncertainties must be reported as the fit result, and the confidence contours should be shown as well.

In [7]: def exponential(x, A0=1, tau=1): return A0 * np.exp(-x/tau)

# define the data as simple Python lists x = [8.0e-01, 2.34e+00, 1.44e+00, 1.28e+00, 2.84e+00, 3.49e+00, 3.78e+00, 4.56e+00, 4.48e+00, 5.38e+00] xerr = 3.0e-01 y = [4.35e-01, 1.51e-01, 8.1e-02, 1.7e-01, 5.3e-02, 1.98e-02, 2.07e-02, 1.230e-02, 9.7e-03, 2.41e-03] yerr = [9.3e-02, 4.1e-02, 3.8e-02, 3.9e-02, 3.1e-02, 1.04e-02, 0.95e-02, 0.89e-02, 0.79e-02, 0.20e-02]

# create a fit object from the data arrays fit = Fit(data=[x, y], model_function=exponential) fit.add_error(axis='x', err_val=xerr) # add the x-error to the fit fit.add_error(axis='y', err_val=yerr) # add the y-errors to the fit

fit.do_fit() # perform the fit fit.report(asymmetric_parameter_errors=True) # print a report with asymmetric uncertainties

# Optional: create a plot plot = Plot(fit) plot.plot(asymmetric_parameter_errors=True, ratio=True) # add the ratio data/function and asymmetric e

# Optional: create the contours profiler cpf = ContoursProfiler(fit) cpf.plot_profiles_contours_matrix() # plot the contour profile matrix for all parameters

plt.show()

######## # Data # ########

X Data ====== 0.8 2.34 1.44 1.28 2.84 3.49 3.78 4.56 4.48 5.38

X Data Error ============ 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3

X Data Correlation Matrix =============================== [1. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 1. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 1. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 1. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 1. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 1. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 1. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 1. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 1. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 1.]

Y Data ======= 0.435 0.151 0.081 0.17 0.053 0.0198 0.0207 0.0123 0.0097 0.00241

Y Data Error ============ 0.093 0.041 0.038 0.039 0.031 0.0104 0.0095 0.0089 0.0079 0.002

Y Data Correlation Matrix =============================== [1. 0. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 1. 0. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 1. 0. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 1. 0. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 1. 0. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 1. 0. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 1. 0. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 1. 0. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 1. 0.] [0. 0. 0. 0. 0. 0. 0. 0. 0. 1.]

######### # Model # #########

Model Function ==============

exponential(x; A0, tau)

X Model ======= 0.8 2.34 1.44 1.28 2.84 3.49 3.78 4.56 4.48 5.38

Y Model ======== 0.3144 0.06921 0.1676 0.1961 0.04234 0.02235 0.01681 0.00781 0.008449 0.003489

############### # Fit Results # ###############

Model Parameters ================

A0 = 0.69 + 0.28 (up) - 0.19 (down) tau = 1.02 + 0.13 (up) - 0.12 (down)

Model Parameter Correlations ============================

A0

tau

======= =======

A0 1.0

-0.8336

tau -0.8336 1.0

Cost Function =============

Cost function: chi-square (with covariance matrix)

chi2 / ndf = 6.873 / 8 = 0.8591

chi2 probability = 0.550

2.3 Practical example with non-trivial covariance matrix: Characteristics of a diode

In many cases the set of uncertainties needed to describe the behaviour of measruements is not trivial. There may be independent and/or correlated absolute and/or relative uncertainties in the x- and/or y-direction. Typically, more than one of these eight kinds of uncertainties are relevant for a given problem, and getting the covariance matrix right is often not an easy task. kafe2 supports constructing the full covariance matrix with the method

add_error(err_val, axis=?, name=None, correlation=0, relative=False, reference='data') ,

which allows to specify various kinds of uncertainties for all or for sub-sets of the data one after the other. The resulting individual covariance matrices are all added to form the full covariance matrix used in the fit. This covariance matrix is re-calculated from the uncertainty components in each minimization step, if needed.

As an example, consider a typical digital amperemeter or voltmeter. Decive characteristics are specified as 4000 Counts, +/(0.5% + 3 digits), where the first, the relative calibration uncertainty, refers to the true value of the current or voltage and is fully correlated among all measurements, while sedond one, the digitisation uncertainty, is independent for each of the measurements. There often is an additional noise component, which also constitutes a set of independent uncertainties.

The code in the cell below shows how these uncertainty components for a set of voltage and current measurements with a diode can be specified for a fit of the Shockley diode equation. Most of the code is indeed needed to specify the uncertainties, the fit itself and the ouput of the results are very similar to the examples discussed already.

In [8]: '''Fit of Shockley equation to I-U characteristic of a diode '''

# model function to fit def Shockley(U, I_s = 0.5, U0 = 0.03):

"""Parametrisation of a diode characteristic

U0 should be limited such that U/U0 ................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download