pyerrors.fits

  1import gc
  2from collections.abc import Sequence
  3import warnings
  4import numpy as np
  5import autograd.numpy as anp
  6import scipy.optimize
  7import scipy.stats
  8import matplotlib.pyplot as plt
  9from matplotlib import gridspec
 10from scipy.odr import ODR, Model, RealData
 11from scipy.stats import chi2
 12import iminuit
 13from autograd import jacobian
 14from autograd import elementwise_grad as egrad
 15from .obs import Obs, derived_observable, covariance, cov_Obs
 16
 17
 18class Fit_result(Sequence):
 19    """Represents fit results.
 20
 21    Attributes
 22    ----------
 23    fit_parameters : list
 24        results for the individual fit parameters,
 25        also accessible via indices.
 26    """
 27
 28    def __init__(self):
 29        self.fit_parameters = None
 30
 31    def __getitem__(self, idx):
 32        return self.fit_parameters[idx]
 33
 34    def __len__(self):
 35        return len(self.fit_parameters)
 36
 37    def gamma_method(self, **kwargs):
 38        """Apply the gamma method to all fit parameters"""
 39        [o.gamma_method(**kwargs) for o in self.fit_parameters]
 40
 41    def __str__(self):
 42        my_str = 'Goodness of fit:\n'
 43        if hasattr(self, 'chisquare_by_dof'):
 44            my_str += '\u03C7\u00b2/d.o.f. = ' + f'{self.chisquare_by_dof:2.6f}' + '\n'
 45        elif hasattr(self, 'residual_variance'):
 46            my_str += 'residual variance = ' + f'{self.residual_variance:2.6f}' + '\n'
 47        if hasattr(self, 'chisquare_by_expected_chisquare'):
 48            my_str += '\u03C7\u00b2/\u03C7\u00b2exp  = ' + f'{self.chisquare_by_expected_chisquare:2.6f}' + '\n'
 49        if hasattr(self, 'p_value'):
 50            my_str += 'p-value   = ' + f'{self.p_value:2.4f}' + '\n'
 51        my_str += 'Fit parameters:\n'
 52        for i_par, par in enumerate(self.fit_parameters):
 53            my_str += str(i_par) + '\t' + ' ' * int(par >= 0) + str(par).rjust(int(par < 0.0)) + '\n'
 54        return my_str
 55
 56    def __repr__(self):
 57        m = max(map(len, list(self.__dict__.keys()))) + 1
 58        return '\n'.join([key.rjust(m) + ': ' + repr(value) for key, value in sorted(self.__dict__.items())])
 59
 60
 61def least_squares(x, y, func, priors=None, silent=False, **kwargs):
 62    r'''Performs a non-linear fit to y = func(x).
 63
 64    Parameters
 65    ----------
 66    x : list
 67        list of floats.
 68    y : list
 69        list of Obs.
 70    func : object
 71        fit function, has to be of the form
 72
 73        ```python
 74        import autograd.numpy as anp
 75
 76        def func(a, x):
 77            return a[0] + a[1] * x + a[2] * anp.sinh(x)
 78        ```
 79
 80        For multiple x values func can be of the form
 81
 82        ```python
 83        def func(a, x):
 84            (x1, x2) = x
 85            return a[0] * x1 ** 2 + a[1] * x2
 86        ```
 87
 88        It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation
 89        will not work.
 90    priors : list, optional
 91        priors has to be a list with an entry for every parameter in the fit. The entries can either be
 92        Obs (e.g. results from a previous fit) or strings containing a value and an error formatted like
 93        0.548(23), 500(40) or 0.5(0.4)
 94    silent : bool, optional
 95        If true all output to the console is omitted (default False).
 96    initial_guess : list
 97        can provide an initial guess for the input parameters. Relevant for
 98        non-linear fits with many parameters. In case of correlated fits the guess is used to perform
 99        an uncorrelated fit which then serves as guess for the correlated fit.
100    method : str, optional
101        can be used to choose an alternative method for the minimization of chisquare.
102        The possible methods are the ones which can be used for scipy.optimize.minimize and
103        migrad of iminuit. If no method is specified, Levenberg-Marquard is used.
104        Reliable alternatives are migrad, Powell and Nelder-Mead.
105    correlated_fit : bool
106        If True, use the full inverse covariance matrix in the definition of the chisquare cost function.
107        For details about how the covariance matrix is estimated see `pyerrors.obs.covariance`.
108        In practice the correlation matrix is Cholesky decomposed and inverted (instead of the covariance matrix).
109        This procedure should be numerically more stable as the correlation matrix is typically better conditioned (Jacobi preconditioning).
110        At the moment this option only works for `prior==None` and when no `method` is given.
111    expected_chisquare : bool
112        If True estimates the expected chisquare which is
113        corrected by effects caused by correlated input data (default False).
114    resplot : bool
115        If True, a plot which displays fit, data and residuals is generated (default False).
116    qqplot : bool
117        If True, a quantile-quantile plot of the fit result is generated (default False).
118    '''
119    if priors is not None:
120        return _prior_fit(x, y, func, priors, silent=silent, **kwargs)
121    else:
122        return _standard_fit(x, y, func, silent=silent, **kwargs)
123
124
125def total_least_squares(x, y, func, silent=False, **kwargs):
126    r'''Performs a non-linear fit to y = func(x) and returns a list of Obs corresponding to the fit parameters.
127
128    Parameters
129    ----------
130    x : list
131        list of Obs, or a tuple of lists of Obs
132    y : list
133        list of Obs. The dvalues of the Obs are used as x- and yerror for the fit.
134    func : object
135        func has to be of the form
136
137        ```python
138        import autograd.numpy as anp
139
140        def func(a, x):
141            return a[0] + a[1] * x + a[2] * anp.sinh(x)
142        ```
143
144        For multiple x values func can be of the form
145
146        ```python
147        def func(a, x):
148            (x1, x2) = x
149            return a[0] * x1 ** 2 + a[1] * x2
150        ```
151
152        It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation
153        will not work.
154    silent : bool, optional
155        If true all output to the console is omitted (default False).
156    initial_guess : list
157        can provide an initial guess for the input parameters. Relevant for non-linear
158        fits with many parameters.
159    expected_chisquare : bool
160        If true prints the expected chisquare which is
161        corrected by effects caused by correlated input data.
162        This can take a while as the full correlation matrix
163        has to be calculated (default False).
164
165    Notes
166    -----
167    Based on the orthogonal distance regression module of scipy
168    '''
169
170    output = Fit_result()
171
172    output.fit_function = func
173
174    x = np.array(x)
175
176    x_shape = x.shape
177
178    if not callable(func):
179        raise TypeError('func has to be a function.')
180
181    for i in range(25):
182        try:
183            func(np.arange(i), x.T[0])
184        except Exception:
185            pass
186        else:
187            break
188
189    n_parms = i
190    if not silent:
191        print('Fit with', n_parms, 'parameter' + 's' * (n_parms > 1))
192
193    x_f = np.vectorize(lambda o: o.value)(x)
194    dx_f = np.vectorize(lambda o: o.dvalue)(x)
195    y_f = np.array([o.value for o in y])
196    dy_f = np.array([o.dvalue for o in y])
197
198    if np.any(np.asarray(dx_f) <= 0.0):
199        raise Exception('No x errors available, run the gamma method first.')
200
201    if np.any(np.asarray(dy_f) <= 0.0):
202        raise Exception('No y errors available, run the gamma method first.')
203
204    if 'initial_guess' in kwargs:
205        x0 = kwargs.get('initial_guess')
206        if len(x0) != n_parms:
207            raise Exception('Initial guess does not have the correct length: %d vs. %d' % (len(x0), n_parms))
208    else:
209        x0 = [1] * n_parms
210
211    data = RealData(x_f, y_f, sx=dx_f, sy=dy_f)
212    model = Model(func)
213    odr = ODR(data, model, x0, partol=np.finfo(np.float64).eps)
214    odr.set_job(fit_type=0, deriv=1)
215    out = odr.run()
216
217    output.residual_variance = out.res_var
218
219    output.method = 'ODR'
220
221    output.message = out.stopreason
222
223    output.xplus = out.xplus
224
225    if not silent:
226        print('Method: ODR')
227        print(*out.stopreason)
228        print('Residual variance:', output.residual_variance)
229
230    if out.info > 3:
231        raise Exception('The minimization procedure did not converge.')
232
233    m = x_f.size
234
235    def odr_chisquare(p):
236        model = func(p[:n_parms], p[n_parms:].reshape(x_shape))
237        chisq = anp.sum(((y_f - model) / dy_f) ** 2) + anp.sum(((x_f - p[n_parms:].reshape(x_shape)) / dx_f) ** 2)
238        return chisq
239
240    if kwargs.get('expected_chisquare') is True:
241        W = np.diag(1 / np.asarray(np.concatenate((dy_f.ravel(), dx_f.ravel()))))
242
243        if kwargs.get('covariance') is not None:
244            cov = kwargs.get('covariance')
245        else:
246            cov = covariance(np.concatenate((y, x.ravel())))
247
248        number_of_x_parameters = int(m / x_f.shape[-1])
249
250        old_jac = jacobian(func)(out.beta, out.xplus)
251        fused_row1 = np.concatenate((old_jac, np.concatenate((number_of_x_parameters * [np.zeros(old_jac.shape)]), axis=0)))
252        fused_row2 = np.concatenate((jacobian(lambda x, y: func(y, x))(out.xplus, out.beta).reshape(x_f.shape[-1], x_f.shape[-1] * number_of_x_parameters), np.identity(number_of_x_parameters * old_jac.shape[0])))
253        new_jac = np.concatenate((fused_row1, fused_row2), axis=1)
254
255        A = W @ new_jac
256        P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T
257        expected_chisquare = np.trace((np.identity(P_phi.shape[0]) - P_phi) @ W @ cov @ W)
258        if expected_chisquare <= 0.0:
259            warnings.warn("Negative expected_chisquare.", RuntimeWarning)
260            expected_chisquare = np.abs(expected_chisquare)
261        output.chisquare_by_expected_chisquare = odr_chisquare(np.concatenate((out.beta, out.xplus.ravel()))) / expected_chisquare
262        if not silent:
263            print('chisquare/expected_chisquare:',
264                  output.chisquare_by_expected_chisquare)
265
266    fitp = out.beta
267    try:
268        hess = jacobian(jacobian(odr_chisquare))(np.concatenate((fitp, out.xplus.ravel())))
269    except TypeError:
270        raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None
271
272    def odr_chisquare_compact_x(d):
273        model = func(d[:n_parms], d[n_parms:n_parms + m].reshape(x_shape))
274        chisq = anp.sum(((y_f - model) / dy_f) ** 2) + anp.sum(((d[n_parms + m:].reshape(x_shape) - d[n_parms:n_parms + m].reshape(x_shape)) / dx_f) ** 2)
275        return chisq
276
277    jac_jac_x = jacobian(jacobian(odr_chisquare_compact_x))(np.concatenate((fitp, out.xplus.ravel(), x_f.ravel())))
278
279    # Compute hess^{-1} @ jac_jac_x[:n_parms + m, n_parms + m:] using LAPACK dgesv
280    try:
281        deriv_x = -scipy.linalg.solve(hess, jac_jac_x[:n_parms + m, n_parms + m:])
282    except np.linalg.LinAlgError:
283        raise Exception("Cannot invert hessian matrix.")
284
285    def odr_chisquare_compact_y(d):
286        model = func(d[:n_parms], d[n_parms:n_parms + m].reshape(x_shape))
287        chisq = anp.sum(((d[n_parms + m:] - model) / dy_f) ** 2) + anp.sum(((x_f - d[n_parms:n_parms + m].reshape(x_shape)) / dx_f) ** 2)
288        return chisq
289
290    jac_jac_y = jacobian(jacobian(odr_chisquare_compact_y))(np.concatenate((fitp, out.xplus.ravel(), y_f)))
291
292    # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv
293    try:
294        deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms + m, n_parms + m:])
295    except np.linalg.LinAlgError:
296        raise Exception("Cannot invert hessian matrix.")
297
298    result = []
299    for i in range(n_parms):
300        result.append(derived_observable(lambda my_var, **kwargs: (my_var[0] + np.finfo(np.float64).eps) / (x.ravel()[0].value + np.finfo(np.float64).eps) * out.beta[i], list(x.ravel()) + list(y), man_grad=list(deriv_x[i]) + list(deriv_y[i])))
301
302    output.fit_parameters = result
303
304    output.odr_chisquare = odr_chisquare(np.concatenate((out.beta, out.xplus.ravel())))
305    output.dof = x.shape[-1] - n_parms
306    output.p_value = 1 - chi2.cdf(output.odr_chisquare, output.dof)
307
308    return output
309
310
311def _prior_fit(x, y, func, priors, silent=False, **kwargs):
312    output = Fit_result()
313
314    output.fit_function = func
315
316    x = np.asarray(x)
317
318    if not callable(func):
319        raise TypeError('func has to be a function.')
320
321    for i in range(100):
322        try:
323            func(np.arange(i), 0)
324        except Exception:
325            pass
326        else:
327            break
328
329    n_parms = i
330
331    if n_parms != len(priors):
332        raise Exception('Priors does not have the correct length.')
333
334    def extract_val_and_dval(string):
335        split_string = string.split('(')
336        if '.' in split_string[0] and '.' not in split_string[1][:-1]:
337            factor = 10 ** -len(split_string[0].partition('.')[2])
338        else:
339            factor = 1
340        return float(split_string[0]), float(split_string[1][:-1]) * factor
341
342    loc_priors = []
343    for i_n, i_prior in enumerate(priors):
344        if isinstance(i_prior, Obs):
345            loc_priors.append(i_prior)
346        else:
347            loc_val, loc_dval = extract_val_and_dval(i_prior)
348            loc_priors.append(cov_Obs(loc_val, loc_dval ** 2, '#prior' + str(i_n) + f"_{np.random.randint(2147483647):010d}"))
349
350    output.priors = loc_priors
351
352    if not silent:
353        print('Fit with', n_parms, 'parameter' + 's' * (n_parms > 1))
354
355    y_f = [o.value for o in y]
356    dy_f = [o.dvalue for o in y]
357
358    if np.any(np.asarray(dy_f) <= 0.0):
359        raise Exception('No y errors available, run the gamma method first.')
360
361    p_f = [o.value for o in loc_priors]
362    dp_f = [o.dvalue for o in loc_priors]
363
364    if np.any(np.asarray(dp_f) <= 0.0):
365        raise Exception('No prior errors available, run the gamma method first.')
366
367    if 'initial_guess' in kwargs:
368        x0 = kwargs.get('initial_guess')
369        if len(x0) != n_parms:
370            raise Exception('Initial guess does not have the correct length.')
371    else:
372        x0 = p_f
373
374    def chisqfunc(p):
375        model = func(p, x)
376        chisq = anp.sum(((y_f - model) / dy_f) ** 2) + anp.sum(((p_f - p) / dp_f) ** 2)
377        return chisq
378
379    if not silent:
380        print('Method: migrad')
381
382    m = iminuit.Minuit(chisqfunc, x0)
383    m.errordef = 1
384    m.print_level = 0
385    if 'tol' in kwargs:
386        m.tol = kwargs.get('tol')
387    else:
388        m.tol = 1e-4
389    m.migrad()
390    params = np.asarray(m.values)
391
392    output.chisquare_by_dof = m.fval / len(x)
393
394    output.method = 'migrad'
395
396    if not silent:
397        print('chisquare/d.o.f.:', output.chisquare_by_dof)
398
399    if not m.fmin.is_valid:
400        raise Exception('The minimization procedure did not converge.')
401
402    hess_inv = np.linalg.pinv(jacobian(jacobian(chisqfunc))(params))
403
404    def chisqfunc_compact(d):
405        model = func(d[:n_parms], x)
406        chisq = anp.sum(((d[n_parms: n_parms + len(x)] - model) / dy_f) ** 2) + anp.sum(((d[n_parms + len(x):] - d[:n_parms]) / dp_f) ** 2)
407        return chisq
408
409    jac_jac = jacobian(jacobian(chisqfunc_compact))(np.concatenate((params, y_f, p_f)))
410
411    deriv = -hess_inv @ jac_jac[:n_parms, n_parms:]
412
413    result = []
414    for i in range(n_parms):
415        result.append(derived_observable(lambda x, **kwargs: (x[0] + np.finfo(np.float64).eps) / (y[0].value + np.finfo(np.float64).eps) * params[i], list(y) + list(loc_priors), man_grad=list(deriv[i])))
416
417    output.fit_parameters = result
418    output.chisquare = chisqfunc(np.asarray(params))
419
420    if kwargs.get('resplot') is True:
421        residual_plot(x, y, func, result)
422
423    if kwargs.get('qqplot') is True:
424        qqplot(x, y, func, result)
425
426    return output
427
428
429def _standard_fit(x, y, func, silent=False, **kwargs):
430
431    output = Fit_result()
432
433    output.fit_function = func
434
435    x = np.asarray(x)
436
437    if x.shape[-1] != len(y):
438        raise Exception('x and y input have to have the same length')
439
440    if len(x.shape) > 2:
441        raise Exception('Unknown format for x values')
442
443    if not callable(func):
444        raise TypeError('func has to be a function.')
445
446    for i in range(25):
447        try:
448            func(np.arange(i), x.T[0])
449        except Exception:
450            pass
451        else:
452            break
453
454    n_parms = i
455
456    if not silent:
457        print('Fit with', n_parms, 'parameter' + 's' * (n_parms > 1))
458
459    y_f = [o.value for o in y]
460    dy_f = [o.dvalue for o in y]
461
462    if np.any(np.asarray(dy_f) <= 0.0):
463        raise Exception('No y errors available, run the gamma method first.')
464
465    if 'initial_guess' in kwargs:
466        x0 = kwargs.get('initial_guess')
467        if len(x0) != n_parms:
468            raise Exception('Initial guess does not have the correct length: %d vs. %d' % (len(x0), n_parms))
469    else:
470        x0 = [0.1] * n_parms
471
472    if kwargs.get('correlated_fit') is True:
473        corr = covariance(y, correlation=True, **kwargs)
474        covdiag = np.diag(1 / np.asarray(dy_f))
475        condn = np.linalg.cond(corr)
476        if condn > 0.1 / np.finfo(float).eps:
477            raise Exception(f"Cannot invert correlation matrix as its condition number exceeds machine precision ({condn:1.2e})")
478        if condn > 1 / np.sqrt(np.finfo(float).eps):
479            warnings.warn("Correlation matrix may be ill-conditioned, condition number: {%1.2e}" % (condn), RuntimeWarning)
480        chol = np.linalg.cholesky(corr)
481        chol_inv = scipy.linalg.solve_triangular(chol, covdiag, lower=True)
482
483        def chisqfunc_corr(p):
484            model = func(p, x)
485            chisq = anp.sum(anp.dot(chol_inv, (y_f - model)) ** 2)
486            return chisq
487
488    def chisqfunc(p):
489        model = func(p, x)
490        chisq = anp.sum(((y_f - model) / dy_f) ** 2)
491        return chisq
492
493    output.method = kwargs.get('method', 'Levenberg-Marquardt')
494    if not silent:
495        print('Method:', output.method)
496
497    if output.method != 'Levenberg-Marquardt':
498        if output.method == 'migrad':
499            fit_result = iminuit.minimize(chisqfunc, x0, tol=1e-4)  # Stopping criterion 0.002 * tol * errordef
500            if kwargs.get('correlated_fit') is True:
501                fit_result = iminuit.minimize(chisqfunc_corr, fit_result.x, tol=1e-4)  # Stopping criterion 0.002 * tol * errordef
502            output.iterations = fit_result.nfev
503        else:
504            fit_result = scipy.optimize.minimize(chisqfunc, x0, method=kwargs.get('method'), tol=1e-12)
505            if kwargs.get('correlated_fit') is True:
506                fit_result = scipy.optimize.minimize(chisqfunc_corr, fit_result.x, method=kwargs.get('method'), tol=1e-12)
507            output.iterations = fit_result.nit
508
509        chisquare = fit_result.fun
510
511    else:
512        if kwargs.get('correlated_fit') is True:
513            def chisqfunc_residuals_corr(p):
514                model = func(p, x)
515                chisq = anp.dot(chol_inv, (y_f - model))
516                return chisq
517
518        def chisqfunc_residuals(p):
519            model = func(p, x)
520            chisq = ((y_f - model) / dy_f)
521            return chisq
522
523        fit_result = scipy.optimize.least_squares(chisqfunc_residuals, x0, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15)
524        if kwargs.get('correlated_fit') is True:
525            fit_result = scipy.optimize.least_squares(chisqfunc_residuals_corr, fit_result.x, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15)
526
527        chisquare = np.sum(fit_result.fun ** 2)
528        if kwargs.get('correlated_fit') is True:
529            assert np.isclose(chisquare, chisqfunc_corr(fit_result.x), atol=1e-14)
530        else:
531            assert np.isclose(chisquare, chisqfunc(fit_result.x), atol=1e-14)
532
533        output.iterations = fit_result.nfev
534
535    if not fit_result.success:
536        raise Exception('The minimization procedure did not converge.')
537
538    if x.shape[-1] - n_parms > 0:
539        output.chisquare_by_dof = chisquare / (x.shape[-1] - n_parms)
540    else:
541        output.chisquare_by_dof = float('nan')
542
543    output.message = fit_result.message
544    if not silent:
545        print(fit_result.message)
546        print('chisquare/d.o.f.:', output.chisquare_by_dof)
547
548    if kwargs.get('expected_chisquare') is True:
549        if kwargs.get('correlated_fit') is not True:
550            W = np.diag(1 / np.asarray(dy_f))
551            cov = covariance(y)
552            A = W @ jacobian(func)(fit_result.x, x)
553            P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T
554            expected_chisquare = np.trace((np.identity(x.shape[-1]) - P_phi) @ W @ cov @ W)
555            output.chisquare_by_expected_chisquare = chisquare / expected_chisquare
556            if not silent:
557                print('chisquare/expected_chisquare:',
558                      output.chisquare_by_expected_chisquare)
559
560    fitp = fit_result.x
561    try:
562        if kwargs.get('correlated_fit') is True:
563            hess = jacobian(jacobian(chisqfunc_corr))(fitp)
564        else:
565            hess = jacobian(jacobian(chisqfunc))(fitp)
566    except TypeError:
567        raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None
568
569    if kwargs.get('correlated_fit') is True:
570        def chisqfunc_compact(d):
571            model = func(d[:n_parms], x)
572            chisq = anp.sum(anp.dot(chol_inv, (d[n_parms:] - model)) ** 2)
573            return chisq
574
575    else:
576        def chisqfunc_compact(d):
577            model = func(d[:n_parms], x)
578            chisq = anp.sum(((d[n_parms:] - model) / dy_f) ** 2)
579            return chisq
580
581    jac_jac = jacobian(jacobian(chisqfunc_compact))(np.concatenate((fitp, y_f)))
582
583    # Compute hess^{-1} @ jac_jac[:n_parms, n_parms:] using LAPACK dgesv
584    try:
585        deriv = -scipy.linalg.solve(hess, jac_jac[:n_parms, n_parms:])
586    except np.linalg.LinAlgError:
587        raise Exception("Cannot invert hessian matrix.")
588
589    result = []
590    for i in range(n_parms):
591        result.append(derived_observable(lambda x, **kwargs: (x[0] + np.finfo(np.float64).eps) / (y[0].value + np.finfo(np.float64).eps) * fit_result.x[i], list(y), man_grad=list(deriv[i])))
592
593    output.fit_parameters = result
594
595    output.chisquare = chisquare
596    output.dof = x.shape[-1] - n_parms
597    output.p_value = 1 - chi2.cdf(output.chisquare, output.dof)
598
599    if kwargs.get('resplot') is True:
600        residual_plot(x, y, func, result)
601
602    if kwargs.get('qqplot') is True:
603        qqplot(x, y, func, result)
604
605    return output
606
607
608def fit_lin(x, y, **kwargs):
609    """Performs a linear fit to y = n + m * x and returns two Obs n, m.
610
611    Parameters
612    ----------
613    x : list
614        Can either be a list of floats in which case no xerror is assumed, or
615        a list of Obs, where the dvalues of the Obs are used as xerror for the fit.
616    y : list
617        List of Obs, the dvalues of the Obs are used as yerror for the fit.
618    """
619
620    def f(a, x):
621        y = a[0] + a[1] * x
622        return y
623
624    if all(isinstance(n, Obs) for n in x):
625        out = total_least_squares(x, y, f, **kwargs)
626        return out.fit_parameters
627    elif all(isinstance(n, float) or isinstance(n, int) for n in x) or isinstance(x, np.ndarray):
628        out = least_squares(x, y, f, **kwargs)
629        return out.fit_parameters
630    else:
631        raise Exception('Unsupported types for x')
632
633
634def qqplot(x, o_y, func, p):
635    """Generates a quantile-quantile plot of the fit result which can be used to
636       check if the residuals of the fit are gaussian distributed.
637    """
638
639    residuals = []
640    for i_x, i_y in zip(x, o_y):
641        residuals.append((i_y - func(p, i_x)) / i_y.dvalue)
642    residuals = sorted(residuals)
643    my_y = [o.value for o in residuals]
644    probplot = scipy.stats.probplot(my_y)
645    my_x = probplot[0][0]
646    plt.figure(figsize=(8, 8 / 1.618))
647    plt.errorbar(my_x, my_y, fmt='o')
648    fit_start = my_x[0]
649    fit_stop = my_x[-1]
650    samples = np.arange(fit_start, fit_stop, 0.01)
651    plt.plot(samples, samples, 'k--', zorder=11, label='Standard normal distribution')
652    plt.plot(samples, probplot[1][0] * samples + probplot[1][1], zorder=10, label='Least squares fit, r=' + str(np.around(probplot[1][2], 3)), marker='', ls='-')
653
654    plt.xlabel('Theoretical quantiles')
655    plt.ylabel('Ordered Values')
656    plt.legend()
657    plt.draw()
658
659
660def residual_plot(x, y, func, fit_res):
661    """ Generates a plot which compares the fit to the data and displays the corresponding residuals"""
662    sorted_x = sorted(x)
663    xstart = sorted_x[0] - 0.5 * (sorted_x[1] - sorted_x[0])
664    xstop = sorted_x[-1] + 0.5 * (sorted_x[-1] - sorted_x[-2])
665    x_samples = np.arange(xstart, xstop + 0.01, 0.01)
666
667    plt.figure(figsize=(8, 8 / 1.618))
668    gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
669    ax0 = plt.subplot(gs[0])
670    ax0.errorbar(x, [o.value for o in y], yerr=[o.dvalue for o in y], ls='none', fmt='o', capsize=3, markersize=5, label='Data')
671    ax0.plot(x_samples, func([o.value for o in fit_res], x_samples), label='Fit', zorder=10, ls='-', ms=0)
672    ax0.set_xticklabels([])
673    ax0.set_xlim([xstart, xstop])
674    ax0.set_xticklabels([])
675    ax0.legend()
676
677    residuals = (np.asarray([o.value for o in y]) - func([o.value for o in fit_res], x)) / np.asarray([o.dvalue for o in y])
678    ax1 = plt.subplot(gs[1])
679    ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
680    ax1.tick_params(direction='out')
681    ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
682    ax1.axhline(y=0.0, ls='--', color='k', marker=" ")
683    ax1.fill_between(x_samples, -1.0, 1.0, alpha=0.1, facecolor='k')
684    ax1.set_xlim([xstart, xstop])
685    ax1.set_ylabel('Residuals')
686    plt.subplots_adjust(wspace=None, hspace=None)
687    plt.draw()
688
689
690def error_band(x, func, beta):
691    """Returns the error band for an array of sample values x, for given fit function func with optimized parameters beta."""
692    cov = covariance(beta)
693    if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps):
694        warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning)
695
696    deriv = []
697    for i, item in enumerate(x):
698        deriv.append(np.array(egrad(func)([o.value for o in beta], item)))
699
700    err = []
701    for i, item in enumerate(x):
702        err.append(np.sqrt(deriv[i] @ cov @ deriv[i]))
703    err = np.array(err)
704
705    return err
706
707
708def ks_test(objects=None):
709    """Performs a Kolmogorov–Smirnov test for the p-values of all fit object.
710
711    Parameters
712    ----------
713    objects : list
714        List of fit results to include in the analysis (optional).
715    """
716
717    if objects is None:
718        obs_list = []
719        for obj in gc.get_objects():
720            if isinstance(obj, Fit_result):
721                obs_list.append(obj)
722    else:
723        obs_list = objects
724
725    p_values = [o.p_value for o in obs_list]
726
727    bins = len(p_values)
728    x = np.arange(0, 1.001, 0.001)
729    plt.plot(x, x, 'k', zorder=1)
730    plt.xlim(0, 1)
731    plt.ylim(0, 1)
732    plt.xlabel('p-value')
733    plt.ylabel('Cumulative probability')
734    plt.title(str(bins) + ' p-values')
735
736    n = np.arange(1, bins + 1) / np.float64(bins)
737    Xs = np.sort(p_values)
738    plt.step(Xs, n)
739    diffs = n - Xs
740    loc_max_diff = np.argmax(np.abs(diffs))
741    loc = Xs[loc_max_diff]
742    plt.annotate('', xy=(loc, loc), xytext=(loc, loc + diffs[loc_max_diff]), arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
743    plt.draw()
744
745    print(scipy.stats.kstest(p_values, 'uniform'))
class Fit_result(collections.abc.Sequence):
19class Fit_result(Sequence):
20    """Represents fit results.
21
22    Attributes
23    ----------
24    fit_parameters : list
25        results for the individual fit parameters,
26        also accessible via indices.
27    """
28
29    def __init__(self):
30        self.fit_parameters = None
31
32    def __getitem__(self, idx):
33        return self.fit_parameters[idx]
34
35    def __len__(self):
36        return len(self.fit_parameters)
37
38    def gamma_method(self, **kwargs):
39        """Apply the gamma method to all fit parameters"""
40        [o.gamma_method(**kwargs) for o in self.fit_parameters]
41
42    def __str__(self):
43        my_str = 'Goodness of fit:\n'
44        if hasattr(self, 'chisquare_by_dof'):
45            my_str += '\u03C7\u00b2/d.o.f. = ' + f'{self.chisquare_by_dof:2.6f}' + '\n'
46        elif hasattr(self, 'residual_variance'):
47            my_str += 'residual variance = ' + f'{self.residual_variance:2.6f}' + '\n'
48        if hasattr(self, 'chisquare_by_expected_chisquare'):
49            my_str += '\u03C7\u00b2/\u03C7\u00b2exp  = ' + f'{self.chisquare_by_expected_chisquare:2.6f}' + '\n'
50        if hasattr(self, 'p_value'):
51            my_str += 'p-value   = ' + f'{self.p_value:2.4f}' + '\n'
52        my_str += 'Fit parameters:\n'
53        for i_par, par in enumerate(self.fit_parameters):
54            my_str += str(i_par) + '\t' + ' ' * int(par >= 0) + str(par).rjust(int(par < 0.0)) + '\n'
55        return my_str
56
57    def __repr__(self):
58        m = max(map(len, list(self.__dict__.keys()))) + 1
59        return '\n'.join([key.rjust(m) + ': ' + repr(value) for key, value in sorted(self.__dict__.items())])

Represents fit results.

Attributes
  • fit_parameters (list): results for the individual fit parameters, also accessible via indices.
Fit_result()
29    def __init__(self):
30        self.fit_parameters = None
def gamma_method(self, **kwargs)
38    def gamma_method(self, **kwargs):
39        """Apply the gamma method to all fit parameters"""
40        [o.gamma_method(**kwargs) for o in self.fit_parameters]

Apply the gamma method to all fit parameters

Inherited Members
collections.abc.Sequence
index
count
def least_squares(x, y, func, priors=None, silent=False, **kwargs)
 62def least_squares(x, y, func, priors=None, silent=False, **kwargs):
 63    r'''Performs a non-linear fit to y = func(x).
 64
 65    Parameters
 66    ----------
 67    x : list
 68        list of floats.
 69    y : list
 70        list of Obs.
 71    func : object
 72        fit function, has to be of the form
 73
 74        ```python
 75        import autograd.numpy as anp
 76
 77        def func(a, x):
 78            return a[0] + a[1] * x + a[2] * anp.sinh(x)
 79        ```
 80
 81        For multiple x values func can be of the form
 82
 83        ```python
 84        def func(a, x):
 85            (x1, x2) = x
 86            return a[0] * x1 ** 2 + a[1] * x2
 87        ```
 88
 89        It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation
 90        will not work.
 91    priors : list, optional
 92        priors has to be a list with an entry for every parameter in the fit. The entries can either be
 93        Obs (e.g. results from a previous fit) or strings containing a value and an error formatted like
 94        0.548(23), 500(40) or 0.5(0.4)
 95    silent : bool, optional
 96        If true all output to the console is omitted (default False).
 97    initial_guess : list
 98        can provide an initial guess for the input parameters. Relevant for
 99        non-linear fits with many parameters. In case of correlated fits the guess is used to perform
100        an uncorrelated fit which then serves as guess for the correlated fit.
101    method : str, optional
102        can be used to choose an alternative method for the minimization of chisquare.
103        The possible methods are the ones which can be used for scipy.optimize.minimize and
104        migrad of iminuit. If no method is specified, Levenberg-Marquard is used.
105        Reliable alternatives are migrad, Powell and Nelder-Mead.
106    correlated_fit : bool
107        If True, use the full inverse covariance matrix in the definition of the chisquare cost function.
108        For details about how the covariance matrix is estimated see `pyerrors.obs.covariance`.
109        In practice the correlation matrix is Cholesky decomposed and inverted (instead of the covariance matrix).
110        This procedure should be numerically more stable as the correlation matrix is typically better conditioned (Jacobi preconditioning).
111        At the moment this option only works for `prior==None` and when no `method` is given.
112    expected_chisquare : bool
113        If True estimates the expected chisquare which is
114        corrected by effects caused by correlated input data (default False).
115    resplot : bool
116        If True, a plot which displays fit, data and residuals is generated (default False).
117    qqplot : bool
118        If True, a quantile-quantile plot of the fit result is generated (default False).
119    '''
120    if priors is not None:
121        return _prior_fit(x, y, func, priors, silent=silent, **kwargs)
122    else:
123        return _standard_fit(x, y, func, silent=silent, **kwargs)

Performs a non-linear fit to y = func(x).

Parameters
  • x (list): list of floats.
  • y (list): list of Obs.
  • func (object): fit function, has to be of the form

    import autograd.numpy as anp
    
    def func(a, x):
       return a[0] + a[1] * x + a[2] * anp.sinh(x)
    

    For multiple x values func can be of the form

    def func(a, x):
       (x1, x2) = x
       return a[0] * x1 ** 2 + a[1] * x2
    

    It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation will not work.

  • priors (list, optional): priors has to be a list with an entry for every parameter in the fit. The entries can either be Obs (e.g. results from a previous fit) or strings containing a value and an error formatted like 0.548(23), 500(40) or 0.5(0.4)
  • silent (bool, optional): If true all output to the console is omitted (default False).
  • initial_guess (list): can provide an initial guess for the input parameters. Relevant for non-linear fits with many parameters. In case of correlated fits the guess is used to perform an uncorrelated fit which then serves as guess for the correlated fit.
  • method (str, optional): can be used to choose an alternative method for the minimization of chisquare. The possible methods are the ones which can be used for scipy.optimize.minimize and migrad of iminuit. If no method is specified, Levenberg-Marquard is used. Reliable alternatives are migrad, Powell and Nelder-Mead.
  • correlated_fit (bool): If True, use the full inverse covariance matrix in the definition of the chisquare cost function. For details about how the covariance matrix is estimated see pyerrors.obs.covariance. In practice the correlation matrix is Cholesky decomposed and inverted (instead of the covariance matrix). This procedure should be numerically more stable as the correlation matrix is typically better conditioned (Jacobi preconditioning). At the moment this option only works for prior==None and when no method is given.
  • expected_chisquare (bool): If True estimates the expected chisquare which is corrected by effects caused by correlated input data (default False).
  • resplot (bool): If True, a plot which displays fit, data and residuals is generated (default False).
  • qqplot (bool): If True, a quantile-quantile plot of the fit result is generated (default False).
def total_least_squares(x, y, func, silent=False, **kwargs)
126def total_least_squares(x, y, func, silent=False, **kwargs):
127    r'''Performs a non-linear fit to y = func(x) and returns a list of Obs corresponding to the fit parameters.
128
129    Parameters
130    ----------
131    x : list
132        list of Obs, or a tuple of lists of Obs
133    y : list
134        list of Obs. The dvalues of the Obs are used as x- and yerror for the fit.
135    func : object
136        func has to be of the form
137
138        ```python
139        import autograd.numpy as anp
140
141        def func(a, x):
142            return a[0] + a[1] * x + a[2] * anp.sinh(x)
143        ```
144
145        For multiple x values func can be of the form
146
147        ```python
148        def func(a, x):
149            (x1, x2) = x
150            return a[0] * x1 ** 2 + a[1] * x2
151        ```
152
153        It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation
154        will not work.
155    silent : bool, optional
156        If true all output to the console is omitted (default False).
157    initial_guess : list
158        can provide an initial guess for the input parameters. Relevant for non-linear
159        fits with many parameters.
160    expected_chisquare : bool
161        If true prints the expected chisquare which is
162        corrected by effects caused by correlated input data.
163        This can take a while as the full correlation matrix
164        has to be calculated (default False).
165
166    Notes
167    -----
168    Based on the orthogonal distance regression module of scipy
169    '''
170
171    output = Fit_result()
172
173    output.fit_function = func
174
175    x = np.array(x)
176
177    x_shape = x.shape
178
179    if not callable(func):
180        raise TypeError('func has to be a function.')
181
182    for i in range(25):
183        try:
184            func(np.arange(i), x.T[0])
185        except Exception:
186            pass
187        else:
188            break
189
190    n_parms = i
191    if not silent:
192        print('Fit with', n_parms, 'parameter' + 's' * (n_parms > 1))
193
194    x_f = np.vectorize(lambda o: o.value)(x)
195    dx_f = np.vectorize(lambda o: o.dvalue)(x)
196    y_f = np.array([o.value for o in y])
197    dy_f = np.array([o.dvalue for o in y])
198
199    if np.any(np.asarray(dx_f) <= 0.0):
200        raise Exception('No x errors available, run the gamma method first.')
201
202    if np.any(np.asarray(dy_f) <= 0.0):
203        raise Exception('No y errors available, run the gamma method first.')
204
205    if 'initial_guess' in kwargs:
206        x0 = kwargs.get('initial_guess')
207        if len(x0) != n_parms:
208            raise Exception('Initial guess does not have the correct length: %d vs. %d' % (len(x0), n_parms))
209    else:
210        x0 = [1] * n_parms
211
212    data = RealData(x_f, y_f, sx=dx_f, sy=dy_f)
213    model = Model(func)
214    odr = ODR(data, model, x0, partol=np.finfo(np.float64).eps)
215    odr.set_job(fit_type=0, deriv=1)
216    out = odr.run()
217
218    output.residual_variance = out.res_var
219
220    output.method = 'ODR'
221
222    output.message = out.stopreason
223
224    output.xplus = out.xplus
225
226    if not silent:
227        print('Method: ODR')
228        print(*out.stopreason)
229        print('Residual variance:', output.residual_variance)
230
231    if out.info > 3:
232        raise Exception('The minimization procedure did not converge.')
233
234    m = x_f.size
235
236    def odr_chisquare(p):
237        model = func(p[:n_parms], p[n_parms:].reshape(x_shape))
238        chisq = anp.sum(((y_f - model) / dy_f) ** 2) + anp.sum(((x_f - p[n_parms:].reshape(x_shape)) / dx_f) ** 2)
239        return chisq
240
241    if kwargs.get('expected_chisquare') is True:
242        W = np.diag(1 / np.asarray(np.concatenate((dy_f.ravel(), dx_f.ravel()))))
243
244        if kwargs.get('covariance') is not None:
245            cov = kwargs.get('covariance')
246        else:
247            cov = covariance(np.concatenate((y, x.ravel())))
248
249        number_of_x_parameters = int(m / x_f.shape[-1])
250
251        old_jac = jacobian(func)(out.beta, out.xplus)
252        fused_row1 = np.concatenate((old_jac, np.concatenate((number_of_x_parameters * [np.zeros(old_jac.shape)]), axis=0)))
253        fused_row2 = np.concatenate((jacobian(lambda x, y: func(y, x))(out.xplus, out.beta).reshape(x_f.shape[-1], x_f.shape[-1] * number_of_x_parameters), np.identity(number_of_x_parameters * old_jac.shape[0])))
254        new_jac = np.concatenate((fused_row1, fused_row2), axis=1)
255
256        A = W @ new_jac
257        P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T
258        expected_chisquare = np.trace((np.identity(P_phi.shape[0]) - P_phi) @ W @ cov @ W)
259        if expected_chisquare <= 0.0:
260            warnings.warn("Negative expected_chisquare.", RuntimeWarning)
261            expected_chisquare = np.abs(expected_chisquare)
262        output.chisquare_by_expected_chisquare = odr_chisquare(np.concatenate((out.beta, out.xplus.ravel()))) / expected_chisquare
263        if not silent:
264            print('chisquare/expected_chisquare:',
265                  output.chisquare_by_expected_chisquare)
266
267    fitp = out.beta
268    try:
269        hess = jacobian(jacobian(odr_chisquare))(np.concatenate((fitp, out.xplus.ravel())))
270    except TypeError:
271        raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None
272
273    def odr_chisquare_compact_x(d):
274        model = func(d[:n_parms], d[n_parms:n_parms + m].reshape(x_shape))
275        chisq = anp.sum(((y_f - model) / dy_f) ** 2) + anp.sum(((d[n_parms + m:].reshape(x_shape) - d[n_parms:n_parms + m].reshape(x_shape)) / dx_f) ** 2)
276        return chisq
277
278    jac_jac_x = jacobian(jacobian(odr_chisquare_compact_x))(np.concatenate((fitp, out.xplus.ravel(), x_f.ravel())))
279
280    # Compute hess^{-1} @ jac_jac_x[:n_parms + m, n_parms + m:] using LAPACK dgesv
281    try:
282        deriv_x = -scipy.linalg.solve(hess, jac_jac_x[:n_parms + m, n_parms + m:])
283    except np.linalg.LinAlgError:
284        raise Exception("Cannot invert hessian matrix.")
285
286    def odr_chisquare_compact_y(d):
287        model = func(d[:n_parms], d[n_parms:n_parms + m].reshape(x_shape))
288        chisq = anp.sum(((d[n_parms + m:] - model) / dy_f) ** 2) + anp.sum(((x_f - d[n_parms:n_parms + m].reshape(x_shape)) / dx_f) ** 2)
289        return chisq
290
291    jac_jac_y = jacobian(jacobian(odr_chisquare_compact_y))(np.concatenate((fitp, out.xplus.ravel(), y_f)))
292
293    # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv
294    try:
295        deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms + m, n_parms + m:])
296    except np.linalg.LinAlgError:
297        raise Exception("Cannot invert hessian matrix.")
298
299    result = []
300    for i in range(n_parms):
301        result.append(derived_observable(lambda my_var, **kwargs: (my_var[0] + np.finfo(np.float64).eps) / (x.ravel()[0].value + np.finfo(np.float64).eps) * out.beta[i], list(x.ravel()) + list(y), man_grad=list(deriv_x[i]) + list(deriv_y[i])))
302
303    output.fit_parameters = result
304
305    output.odr_chisquare = odr_chisquare(np.concatenate((out.beta, out.xplus.ravel())))
306    output.dof = x.shape[-1] - n_parms
307    output.p_value = 1 - chi2.cdf(output.odr_chisquare, output.dof)
308
309    return output

Performs a non-linear fit to y = func(x) and returns a list of Obs corresponding to the fit parameters.

Parameters
  • x (list): list of Obs, or a tuple of lists of Obs
  • y (list): list of Obs. The dvalues of the Obs are used as x- and yerror for the fit.
  • func (object): func has to be of the form

    import autograd.numpy as anp
    
    def func(a, x):
       return a[0] + a[1] * x + a[2] * anp.sinh(x)
    

    For multiple x values func can be of the form

    def func(a, x):
       (x1, x2) = x
       return a[0] * x1 ** 2 + a[1] * x2
    

    It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation will not work.

  • silent (bool, optional): If true all output to the console is omitted (default False).
  • initial_guess (list): can provide an initial guess for the input parameters. Relevant for non-linear fits with many parameters.
  • expected_chisquare (bool): If true prints the expected chisquare which is corrected by effects caused by correlated input data. This can take a while as the full correlation matrix has to be calculated (default False).
Notes

Based on the orthogonal distance regression module of scipy

def fit_lin(x, y, **kwargs)
609def fit_lin(x, y, **kwargs):
610    """Performs a linear fit to y = n + m * x and returns two Obs n, m.
611
612    Parameters
613    ----------
614    x : list
615        Can either be a list of floats in which case no xerror is assumed, or
616        a list of Obs, where the dvalues of the Obs are used as xerror for the fit.
617    y : list
618        List of Obs, the dvalues of the Obs are used as yerror for the fit.
619    """
620
621    def f(a, x):
622        y = a[0] + a[1] * x
623        return y
624
625    if all(isinstance(n, Obs) for n in x):
626        out = total_least_squares(x, y, f, **kwargs)
627        return out.fit_parameters
628    elif all(isinstance(n, float) or isinstance(n, int) for n in x) or isinstance(x, np.ndarray):
629        out = least_squares(x, y, f, **kwargs)
630        return out.fit_parameters
631    else:
632        raise Exception('Unsupported types for x')

Performs a linear fit to y = n + m * x and returns two Obs n, m.

Parameters
  • x (list): Can either be a list of floats in which case no xerror is assumed, or a list of Obs, where the dvalues of the Obs are used as xerror for the fit.
  • y (list): List of Obs, the dvalues of the Obs are used as yerror for the fit.
def qqplot(x, o_y, func, p)
635def qqplot(x, o_y, func, p):
636    """Generates a quantile-quantile plot of the fit result which can be used to
637       check if the residuals of the fit are gaussian distributed.
638    """
639
640    residuals = []
641    for i_x, i_y in zip(x, o_y):
642        residuals.append((i_y - func(p, i_x)) / i_y.dvalue)
643    residuals = sorted(residuals)
644    my_y = [o.value for o in residuals]
645    probplot = scipy.stats.probplot(my_y)
646    my_x = probplot[0][0]
647    plt.figure(figsize=(8, 8 / 1.618))
648    plt.errorbar(my_x, my_y, fmt='o')
649    fit_start = my_x[0]
650    fit_stop = my_x[-1]
651    samples = np.arange(fit_start, fit_stop, 0.01)
652    plt.plot(samples, samples, 'k--', zorder=11, label='Standard normal distribution')
653    plt.plot(samples, probplot[1][0] * samples + probplot[1][1], zorder=10, label='Least squares fit, r=' + str(np.around(probplot[1][2], 3)), marker='', ls='-')
654
655    plt.xlabel('Theoretical quantiles')
656    plt.ylabel('Ordered Values')
657    plt.legend()
658    plt.draw()

Generates a quantile-quantile plot of the fit result which can be used to check if the residuals of the fit are gaussian distributed.

def residual_plot(x, y, func, fit_res)
661def residual_plot(x, y, func, fit_res):
662    """ Generates a plot which compares the fit to the data and displays the corresponding residuals"""
663    sorted_x = sorted(x)
664    xstart = sorted_x[0] - 0.5 * (sorted_x[1] - sorted_x[0])
665    xstop = sorted_x[-1] + 0.5 * (sorted_x[-1] - sorted_x[-2])
666    x_samples = np.arange(xstart, xstop + 0.01, 0.01)
667
668    plt.figure(figsize=(8, 8 / 1.618))
669    gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
670    ax0 = plt.subplot(gs[0])
671    ax0.errorbar(x, [o.value for o in y], yerr=[o.dvalue for o in y], ls='none', fmt='o', capsize=3, markersize=5, label='Data')
672    ax0.plot(x_samples, func([o.value for o in fit_res], x_samples), label='Fit', zorder=10, ls='-', ms=0)
673    ax0.set_xticklabels([])
674    ax0.set_xlim([xstart, xstop])
675    ax0.set_xticklabels([])
676    ax0.legend()
677
678    residuals = (np.asarray([o.value for o in y]) - func([o.value for o in fit_res], x)) / np.asarray([o.dvalue for o in y])
679    ax1 = plt.subplot(gs[1])
680    ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
681    ax1.tick_params(direction='out')
682    ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
683    ax1.axhline(y=0.0, ls='--', color='k', marker=" ")
684    ax1.fill_between(x_samples, -1.0, 1.0, alpha=0.1, facecolor='k')
685    ax1.set_xlim([xstart, xstop])
686    ax1.set_ylabel('Residuals')
687    plt.subplots_adjust(wspace=None, hspace=None)
688    plt.draw()

Generates a plot which compares the fit to the data and displays the corresponding residuals

def error_band(x, func, beta)
691def error_band(x, func, beta):
692    """Returns the error band for an array of sample values x, for given fit function func with optimized parameters beta."""
693    cov = covariance(beta)
694    if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps):
695        warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning)
696
697    deriv = []
698    for i, item in enumerate(x):
699        deriv.append(np.array(egrad(func)([o.value for o in beta], item)))
700
701    err = []
702    for i, item in enumerate(x):
703        err.append(np.sqrt(deriv[i] @ cov @ deriv[i]))
704    err = np.array(err)
705
706    return err

Returns the error band for an array of sample values x, for given fit function func with optimized parameters beta.

def ks_test(objects=None)
709def ks_test(objects=None):
710    """Performs a Kolmogorov–Smirnov test for the p-values of all fit object.
711
712    Parameters
713    ----------
714    objects : list
715        List of fit results to include in the analysis (optional).
716    """
717
718    if objects is None:
719        obs_list = []
720        for obj in gc.get_objects():
721            if isinstance(obj, Fit_result):
722                obs_list.append(obj)
723    else:
724        obs_list = objects
725
726    p_values = [o.p_value for o in obs_list]
727
728    bins = len(p_values)
729    x = np.arange(0, 1.001, 0.001)
730    plt.plot(x, x, 'k', zorder=1)
731    plt.xlim(0, 1)
732    plt.ylim(0, 1)
733    plt.xlabel('p-value')
734    plt.ylabel('Cumulative probability')
735    plt.title(str(bins) + ' p-values')
736
737    n = np.arange(1, bins + 1) / np.float64(bins)
738    Xs = np.sort(p_values)
739    plt.step(Xs, n)
740    diffs = n - Xs
741    loc_max_diff = np.argmax(np.abs(diffs))
742    loc = Xs[loc_max_diff]
743    plt.annotate('', xy=(loc, loc), xytext=(loc, loc + diffs[loc_max_diff]), arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
744    plt.draw()
745
746    print(scipy.stats.kstest(p_values, 'uniform'))

Performs a Kolmogorov–Smirnov test for the p-values of all fit object.

Parameters
  • objects (list): List of fit results to include in the analysis (optional).