Fitting models to data

The module nbi_stat provides the function nbi_stat.fit_plot() to visualise data and a function fitted to the data in a straight forward manner. We will not describe that function here, but we will use it throughout the examples below.

\(\chi^2\) calculation

The function nbi_stat.chi2nu() calculates the \(\chi^2\) and number degrees of freedom given a sample, a fitted model function, and the best-fit parameters.

>>> import nbi_stat
>>> np.random.seed(1234567)
>>> x = np.linspace(0,10,51)
>>> y = x*(1+np.random.normal(0,.1,size=x.shape))
>>> e = y-x
>>> f = x
>>> plt.errorbar(x,y,e,fmt='.',label=r'$y=x+\delta$')
>>> plt.plot(x,f,'-',label=r'$y=x$')
>>> chi2,nu = nbi_stat.chi2nu(x,y,lambda x,*p: x,[],e)
>>> plt.text(2,8,r'$\chi^2/\nu={:.1f}/{}={:.1f}$'
...          .format(chi2,nu,chi2/nu))

(Source code, png, hires.png, pdf)

_images/fitting-1.png

Linear fitting

To fit a linear model - that is a model of the form

\[f(x;p) = \sum_{i=1}^n p_i f_i(x)\]

where

\[\forall j : \frac{\partial f_i}{\partial p_j} = 0\]

the module nbi_stat provides the function nbi_stat.linfit(). This function accepts the independent variable \(x\), the dependent variable \(y\), and the list of \(f_i\). Optionally, one can specify the uncertainties on \(y\).

Let’s take an example. Suppose we want to fit the model

\[f(x;a,b) = ax^b\]

which isn’t linear in the parameter \(b\). However, we can rewrite this as

\[\log\left(f(x;a,b)\right) = \log(a)+b\log(x)\]

which is linear in \(\log(a)\) and \(b\). Thus, to fit this model to data, we need to take the logarithm of both \(x\) and \(y\): and fit a straight line to it.

>>> import nbi_stat
>>> np.random.seed(123456)
>>> l        = 10
>>> b        = 0.3
>>> a        = (b+1)/(l**(b+1)-1)
>>> x        =  np.linspace(1,l,20)
>>> r        = ((a+(b+1)*np.random.random(size=8000))/a)**(1/(b+1))
>>> h,c,w,e  = nbi_stat.histogram(r,np.linspace(1,l,20+1),normalize=True)
>>>
>>> f        = [lambda logx: np.ones_like(logx), lambda logx: logx]
>>>
>>> p, covar = nbi_stat.linfit(f,np.log(c),np.log(h),e)
>>> g        = lambda x,loga,b: np.exp(loga)*x**b
>>>
>>> plt.plot(x,a*x**b,'--',label=r'${:.2f}x^{{{:.1f}}}$'.format(a,b))
>>> nbi_stat.plot_fit(c,h,e,g,p,covar,
...                  parameters=[r"\log(a)",r"b"],
...                  data={'label':'Data','fmt':'o'},
...                  fit={'label':'Fit'},
...                  table={'loc':'lower right'},
...                  legend={'loc':'upper left'})

(Source code, png, hires.png, pdf)

_images/fitting-2.png

Non-linear fitting

The module nbi_stat provides the drop-in function nbi_stat.curve_fit() for scipy.optimize.curve_fit, which allows for uncertainties in the independent variable (\(x\)). Here, we will not will first show a regular fit.

>>> import nbi_stat
>>> y = np.array([ 7,  2,  6, 12,  15,  18,  31,  29,  27,  27, 41, 35,
...               37, 37, 63, 71, 102,  95, 115, 202, 190, 113, 86, 68,
...               74, 79, 75, 79,  68,  62,  69,  81,  79,  85, 87, 68,
...               70, 89, 77, 70,  71,  62,  85,  62,  73,  70, 59, 61,
...               77, 61, 62, 73,  67,  71,  75,  66,  73,  71, 71, 49])
>>> x = np.linspace(0,3,len(y),endpoint=False)
>>> ey = np.sqrt(y)
>>>
>>> def peakBackground(x,a1,a2,a3,a0,gamma,e0):
>>>     return a1+a2*x+a3*x**2+a0*(gamma/(2*np.pi))/((x-e0)**2+(gamma/2)**2)
>>>
>>> p0 = np.array([0,0,0,1000,.1,1])
>>> p, pcov = nbi_stat.curve_fit(peakBackground,x,y,p0,ey)
>>>
>>> nbi_stat.plot_fit(x,y,ey,peakBackground,p,pcov,nsig=1,
...                  parameters=["a_1","a_2","a_3","A_0",
...                              {'label':r"\Gamma",'unit'='GeV'},
...                              {'label':r"E_0",   'unit'='GeV'}],
...                  data={"label":"Data","fmt":"."},
...                  fit={"label":"Fit"},
...                  legend={"loc":"upper left"});

(Source code)

To fit a function to data with uncertainties on the independent variable, we can do

>>> import nbi_stat
>>> data = np.array ([[22000 , 440, -4.017 ,  0.5],
...                   [22930 , 470, -2.742 ,  0.25],
...                   [23880 , 500, -1.1478 , 0.08],
...                   [25130 , 530,  1.491 ,  0.09],
...                   [26390 , 540,  6.873 ,  1.90]])
>>>
>>> def f(omega,a,b): return a*omega + b/omega
>>>
>>> omega,domega = data[:,0], data[:,1]
>>> cott,dcott   = data[:,2], data[:,3]
>>> p, cov = nbi_stat.curve_fit(f,omega,cott,(1,1),dcott,domega,
...                             absolute_sigma=True)
>>>
>>> nbi_stat.plot_fit(omega, cott, dcott, f, p, cov, domega,
...                  parameters=['a','b'],
...                  fit ={ 'label':'Fit'},
...                  data ={'fmt':'none','label':'Data'})

(Source code, png, hires.png, pdf)

_images/fitting-4.png

Maximum likelihood fits

The function nbi_stat.mlefit() performs a Maximum Likelihood Estimator fit of a PDF (or log-PDF) to data.

The function can do

  • Regular MLE

  • Extended MLE (E-MLE)

  • Binned MLE (B-MLE)

  • Extended and Binned MLE (EB-MLE)

More on these modes will follow. The function nbi_stat.llh() (nbi_stat.binned_llh()) calculate the, possibly extended, (binned) likelihood, and can be used stand-alone - for example in likelihood ratios.

See also the description of mle_fit() for more on the different estimates.

>>> import nbi_stat
>>> b = np.linspace(0,3,16)
>>> x = (b[1:]+b[:-1])/2
>>> y = np.random.exponential(size=100)
>>> y = y[y<3]
>>> h = np.histogram(y,b,density=True)[0]
>>> e = np.sqrt(h/np.diff(b)/len(y))
>>>
>>> f = lambda t,tau: \
...     -(np.log(tau)+np.log(1-np.exp(-3/tau)))-t/tau if tau > 0 else -np.inf
>>> g = lambda t,tau: np.exp(-t/tau)/(1-np.exp(-3/tau))
>>>
>>> p, cov = nbi_stat.mlefit(f,y,[1],logpdf=True)
>>> nbi_stat.plot_fit(x,h,e,g,p,cov,
...                  data={'fmt':'none','label':'Data'},
...                  fit={'label':'Fit'},
...                  legend={'loc':'lower center'})

(Source code)

Unified fit interface

The function nbi_stat.fit() provides us with a unified interface for fitting models to data.

>>> import nbi_stat
>>> xlin = np.arange(0,105,5)
>>> ylin = [-.849, -.738, -.537, -.354, -.196, -.019, 0.262,
...         0.413, 0.734, 0.882, 1.258, 1.305, 1.541, 1.768,
...         1.935, 2.147, 2.456, 2.676, 2.994, 3.200, 3.318]
>>> ylsq = [  7,   2,   6,  12,  15,  18,  31,  29,  27,  27,  41,  35,
...          37,  37,  63,  71, 102,  95, 115, 202, 190, 113,  86,  68,
...          74,  79,  75,  79,  68,  62,  69,  81,  79,  85,  87,  68,
...          70,  89,  77,  70,  71,  62,  85,  62,  73,  70,  59,  61,
...          77,  61,  62,  73,  67,  71,  75,  66,  73,  71,  71,  49]
>>> xlsq = np.linspace(0,3,len(ylsq),endpoint=False)
>>> bmle = np.linspace(0,3,16)
>>> xmle = (bmle[1:]+bmle[:-1])/2
>>> ymle = np.random.exponential(size=100)
>>> ymle = ymle[ymle<3]
>>>
>>> hmle = np.histogram(ymle,bmle,density=True)[0]
>>> emle = np.sqrt(hmle/np.diff(bmle)/len(ymle))
>>>
>>> flin = [lambda x: np.ones_like(x), lambda x: x, lambda x: x**2]
>>> plin = lambda x,*p : p[0]+p[1]*x+p[2]*x**2
>>> flsq = lambda x,a1,a2,a3,a0,gamma,e0: \
...     a1+a2*x+a3*x**2+a0*(gamma/(2*np.pi))/((x-e0)**2+(gamma/2)**2)
>>> fmle = lambda t,tau: \
...     -(np.log(tau)+np.log(1-np.exp(-3/tau)))-t/tau if tau > 0 else -np.inf
>>> pmle = lambda t,tau: np.exp(-t/tau)/(1-np.exp(-3/tau))
>>>
>>> tests = [{'t': 'LIN',
...           'f': flin,
...           'args': (xlin,                                    # X
...                    np.array(ylin),                          # Y
...                    np.ones_like(ylin)*0.05),                # delta
...           'F':    plin},
...          {'t': 'LSQ',
...           'f': flsq,
...           'args': (xlsq,                                      # X
...                    np.array(ylsq),                            # Y
...                    np.array([0,0,0,1000,.1,1]),               # p0
...                    np.sqrt(ylsq)),                            # delta
...           'plot': (xlsq,ylsq,np.sqrt(ylsq)) },
...          {'t': 'MLE',
...           'f': fmle,
...           'args': (ymle,                                     # X
...                     [1]),                                     # p0
...           'kwargs': {'logpdf':True},
...           'F':    pmle,
...           'plot': (xmle,hmle,emle)
...          }]
>>>
>>> fig, ax = plt.subplots(ncols=len(tests),figsize=(10,5))
>>> for t, a in zip(tests,ax):
...     p, cov = nbi_stat.fit(t['f'],*t['args'],**t.get('kwargs',{}))
...
...     nbi_stat.plot_fit(*t.get('plot',t['args']),
...                      t.get('F', t['f']), p, cov, axes=a,
...                      data={'fmt':'none','label':'Data'},
...                      fit={'label':'Fit'},
...                      legend={'loc':'lower center'})
...     a.set_title(t['t'])

(Source code, png, hires.png, pdf)

_images/fitting-6.png

Confidence contours

Using the function nbi_stat.nsigma_contour2() one can calculate the \(n\sigma\) confidence contours of any two parameters using the parameter values, uncertainties, and correlation coefficient.

The function nbi_stat.nsigma_contour() calculates _all_ \(n\sigma\) confidence contours of a fit (or similar) from the vector of parameter values and their covariance matrix. The function returns a triangular list of the confidence contours.

>>> import nbi_stat
>>> y = np.array([ 7,  2,  6, 12,  15,  18,  31,  29,  27,  27, 41, 35,
...               37, 37, 63, 71, 102,  95, 115, 202, 190, 113, 86, 68,
...               74, 79, 75, 79,  68,  62,  69,  81,  79,  85, 87, 68,
...               70, 89, 77, 70,  71,  62,  85,  62,  73,  70, 59, 61,
...               77, 61, 62, 73,  67,  71,  75,  66,  73,  71, 71, 49])
>>> x = np.linspace(0,3,len(y),endpoint=False)
>>> ey = np.sqrt(y)
>>>
>>> def peakBackground(x,a1,a2,a3,a0,gamma,e0):
>>>     return a1+a2*x+a3*x**2+a0*(gamma/(2*np.pi))/((x-e0)**2+(gamma/2)**2)
>>>
>>> p0 = np.array([0,0,0,1000,.1,1])
>>> p, pcov = nbi_stat.curve_fit(peakBackground,x,y,p0,ey)
>>> cont   = nbi_stat.nsigma_contour(p,pcov)
>>> n      = len(cont)
>>>
>>> fig,ax = plt.subplots(ncols=n,nrows=n,
...                       sharey='row',sharex='col',
...                       gridspec_kw=dict(wspace=0,hspace=0))
>>> names  = ["a_1","a_2","a_3","A_0",r"\Gamma","E_0"]
>>> for i,j in zip(*np.triu_indices(n,1)):
...     ax[i,j].remove()
>>>
>>> for i, (l, ar, ny) in enumerate(zip(cont, ax, names[1:])):
...     for j, (c, a, nx),  in enumerate(zip(l, ar, names[:-1])):
...         a.plot(c[:,0],c[:,1],'-')
...         if j == 0:
...             a.set_ylabel(fr'${ny}$')
...         if i == n - 1:
...             a.set_xlabel(fr'${nx}$')

(Source code, png, hires.png, pdf)

_images/fitting-7.png