Getting started¶

[1]:

import pyobs
import numpy

pyobs.set_verbose('mfit') # prints additional information
pyobs.set_verbose('diff') # prints additional information

import matplotlib.pyplot as plt
%matplotlib notebook

We create synthetic data from a model distribution that we try to fit using pyobs.

The model function is a simple polynomial $f = p_0 + p_1 x$

[2]:

# true parameters
P1 = [2.3, 4.5]

The utility function diff provides automatically first and second derivatives, and is based on symbolic differentiation

[3]:

xax1=[1,2,3,4,5,6,7,8]
[f1, df1, _] = pyobs.symbolic.diff('p0 + p1*x','x','p0,p1')

$\displaystyle p_{0} + p_{1} x$

$\displaystyle \left[ 1, \ x\right]$

$\displaystyle \left[ \left[ 0, \ 0\right], \ \left[ 0, \ 0\right]\right]$

Now let’s generate synthetic autocorrelated data without correlations, ie all points are uncorrelated: we do so by passing a cov matrix which is simply a vector. The function acrandn interprets it as a diagonal covariance matrix.

[4]:

# generate autocorrelated data
N=500 # number of configs

# central values
mu=[f1(x,P1[0],P1[1]) for x in xax1]

# diagonal cov matrix
cov=[(mu[i]*0.2)**2 for i in range(len(mu))]

# autocorrelation time
tau=4.0

rng = pyobs.random.generator('fits')
data = rng.markov_chain(mu,cov,tau,N)
data = numpy.reshape(data,numpy.size(data),)

Random generator initialized with seed = 2874984293 [fits]

[5]:

yobs1 = pyobs.observable(description='test-data')
yobs1.create('A',data,shape=(len(mu),))
print(yobs1)

7.2(1.1)        11.7(2.1)       14.6(3.2)       16.3(3.8)       24.6(4.7)       21.9(4.3)       38.6(5.9)       32.2(7.2)

[6]:

plt.figure()
[y,dy] = yobs1.error()
plt.errorbar(xax1,y,dy,fmt='.')

[6]:

<ErrorbarContainer object of 3 artists>

Now we create one instance of the mfit class. We need the kinematic coordinates xax1, the function and its derivatives obtained from diff, and the metric of the $\chi^2$, the matrix $W$, which is taken to be the inverse covariance matrix for a proper statistical interpretation of the $\chi^2$. In this case, since we have uncorrelated data points, we only needs their errors.

[7]:

[y,dy] = yobs1.error()
W=1./dy**2
fit1 = pyobs.mfit(xax1,W,f1,df1)

The minimization of the $\chi^2$ is performed by simply calling the mfit class with the corresponding observable as argument. The output is new observable with the fit parameters.

[8]:

pars = fit1(yobs1)
print('fit1 = ', pars)
print('Ratio with true parameters:',pars/numpy.array(P1))

chisquare = 3.8841897075027094
minimizer iterations = 2
minimizer status: Levenberg-Marquardt: converged 1.0e-06 per-cent tolerance on fun
mfit executed in 0.00343609 secs
fit1 =  3.4(1.3)        3.83(56)

Ratio with true parameters: 1.46(56)    0.85(12)

Using the fit parameters we can evaluate the fit function at arbitrary coordinates and check its compatibility with the original data points, or extrapolate values from it.

[9]:

plt.figure()
[y,dy] = yobs1.error()
plt.errorbar(xax1,y,dy,fmt='.',color='C0')

xax=numpy.arange(0,9,0.2)
yeval = fit1.eval(xax,pars)
[y,dy] = yeval.error()
plt.fill_between(xax,y+dy,y-dy,lw=0.1,alpha=0.5,color='C1')

[yzero] = fit1.eval(numpy.array([1.0]), pars)
print('extrapolated value at 1.0 = ', yzero)

extrapolated value at 1.0 =  7.20(93)

[ ]: