# pyerrors.fits API documentation

# pyerrors.fits

View Source

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

class Fit_result(collections.abc.Sequence): View Source

    21class Fit_result(Sequence):
    22    """Represents fit results.
    23
    24    Attributes
    25    ----------
    26    fit_parameters : list
    27        results for the individual fit parameters,
    28        also accessible via indices.
    29    chisquare_by_dof : float
    30        reduced chisquare.
    31    p_value : float
    32        p-value of the fit
    33    t2_p_value : float
    34        Hotelling t-squared p-value for correlated fits.
    35    """
    36
    37    def __init__(self):
    38        self.fit_parameters = None
    39
    40    def __getitem__(self, idx):
    41        return self.fit_parameters[idx]
    42
    43    def __len__(self):
    44        return len(self.fit_parameters)
    45
    46    def gamma_method(self, **kwargs):
    47        """Apply the gamma method to all fit parameters"""
    48        [o.gamma_method(**kwargs) for o in self.fit_parameters]
    49
    50    gm = gamma_method
    51
    52    def __str__(self):
    53        my_str = 'Goodness of fit:\n'
    54        if hasattr(self, 'chisquare_by_dof'):
    55            my_str += '\u03C7\u00b2/d.o.f. = ' + f'{self.chisquare_by_dof:2.6f}' + '\n'
    56        elif hasattr(self, 'residual_variance'):
    57            my_str += 'residual variance = ' + f'{self.residual_variance:2.6f}' + '\n'
    58        if hasattr(self, 'chisquare_by_expected_chisquare'):
    59            my_str += '\u03C7\u00b2/\u03C7\u00b2exp  = ' + f'{self.chisquare_by_expected_chisquare:2.6f}' + '\n'
    60        if hasattr(self, 'p_value'):
    61            my_str += 'p-value   = ' + f'{self.p_value:2.4f}' + '\n'
    62        if hasattr(self, 't2_p_value'):
    63            my_str += 't\u00B2p-value = ' + f'{self.t2_p_value:2.4f}' + '\n'
    64        my_str += 'Fit parameters:\n'
    65        for i_par, par in enumerate(self.fit_parameters):
    66            my_str += str(i_par) + '\t' + ' ' * int(par >= 0) + str(par).rjust(int(par < 0.0)) + '\n'
    67        return my_str
    68
    69    def __repr__(self):
    70        m = max(map(len, list(self.__dict__.keys()))) + 1
    71        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.
- chisquare_by_dof (float): reduced chisquare.
- p_value (float): p-value of the fit
- t2_p_value (float): Hotelling t-squared p-value for correlated fits.

def gamma_method(self, **kwargs): View Source

    46    def gamma_method(self, **kwargs):
    47        """Apply the gamma method to all fit parameters"""
    48        [o.gamma_method(**kwargs) for o in self.fit_parameters]

Apply the gamma method to all fit parameters

def gm(self, **kwargs): View Source

    46    def gamma_method(self, **kwargs):
    47        """Apply the gamma method to all fit parameters"""
    48        [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): View Source

     74def least_squares(x, y, func, priors=None, silent=False, **kwargs):
     75    r'''Performs a non-linear fit to y = func(x).
     76        ```
     77
     78    Parameters
     79    ----------
     80    For an uncombined fit:
     81
     82    x : list
     83        list of floats.
     84    y : list
     85        list of Obs.
     86    func : object
     87        fit function, has to be of the form
     88
     89        ```python
     90        import autograd.numpy as anp
     91
     92        def func(a, x):
     93            return a[0] + a[1] * x + a[2] * anp.sinh(x)
     94        ```
     95
     96        For multiple x values func can be of the form
     97
     98        ```python
     99        def func(a, x):
    100            (x1, x2) = x
    101            return a[0] * x1 ** 2 + a[1] * x2
    102        ```
    103        It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation
    104        will not work.
    105
    106    OR For a combined fit:
    107
    108    x : dict
    109        dict of lists.
    110    y : dict
    111        dict of lists of Obs.
    112    funcs : dict
    113        dict of objects
    114        fit functions have to be of the form (here a[0] is the common fit parameter)
    115        ```python
    116        import autograd.numpy as anp
    117        funcs = {"a": func_a,
    118                "b": func_b}
    119
    120        def func_a(a, x):
    121            return a[1] * anp.exp(-a[0] * x)
    122
    123        def func_b(a, x):
    124            return a[2] * anp.exp(-a[0] * x)
    125
    126        It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation
    127        will not work.
    128
    129    priors : dict or list, optional
    130        priors can either be a dictionary with integer keys and the corresponding priors as values or
    131        a list with an entry for every parameter in the fit. The entries can either be
    132        Obs (e.g. results from a previous fit) or strings containing a value and an error formatted like
    133        0.548(23), 500(40) or 0.5(0.4)
    134    silent : bool, optional
    135        If true all output to the console is omitted (default False).
    136    initial_guess : list
    137        can provide an initial guess for the input parameters. Relevant for
    138        non-linear fits with many parameters. In case of correlated fits the guess is used to perform
    139        an uncorrelated fit which then serves as guess for the correlated fit.
    140    method : str, optional
    141        can be used to choose an alternative method for the minimization of chisquare.
    142        The possible methods are the ones which can be used for scipy.optimize.minimize and
    143        migrad of iminuit. If no method is specified, Levenberg-Marquard is used.
    144        Reliable alternatives are migrad, Powell and Nelder-Mead.
    145    tol: float, optional
    146        can be used (only for combined fits and methods other than Levenberg-Marquard) to set the tolerance for convergence
    147        to a different value to either speed up convergence at the cost of a larger error on the fitted parameters (and possibly
    148        invalid estimates for parameter uncertainties) or smaller values to get more accurate parameter values
    149        The stopping criterion depends on the method, e.g. migrad: edm_max = 0.002 * tol * errordef (EDM criterion: edm < edm_max)
    150    correlated_fit : bool
    151        If True, use the full inverse covariance matrix in the definition of the chisquare cost function.
    152        For details about how the covariance matrix is estimated see `pyerrors.obs.covariance`.
    153        In practice the correlation matrix is Cholesky decomposed and inverted (instead of the covariance matrix).
    154        This procedure should be numerically more stable as the correlation matrix is typically better conditioned (Jacobi preconditioning).
    155    expected_chisquare : bool
    156        If True estimates the expected chisquare which is
    157        corrected by effects caused by correlated input data (default False).
    158    resplot : bool
    159        If True, a plot which displays fit, data and residuals is generated (default False).
    160    qqplot : bool
    161        If True, a quantile-quantile plot of the fit result is generated (default False).
    162    num_grad : bool
    163        Use numerical differentation instead of automatic differentiation to perform the error propagation (default False).
    164
    165    Returns
    166    -------
    167    output : Fit_result
    168        Parameters and information on the fitted result.
    169    '''
    170    output = Fit_result()
    171
    172    if (type(x) == dict and type(y) == dict and type(func) == dict):
    173        xd = {key: anp.asarray(x[key]) for key in x}
    174        yd = y
    175        funcd = func
    176        output.fit_function = func
    177    elif (type(x) == dict or type(y) == dict or type(func) == dict):
    178        raise TypeError("All arguments have to be dictionaries in order to perform a combined fit.")
    179    else:
    180        x = np.asarray(x)
    181        xd = {"": x}
    182        yd = {"": y}
    183        funcd = {"": func}
    184        output.fit_function = func
    185
    186    if kwargs.get('num_grad') is True:
    187        jacobian = num_jacobian
    188        hessian = num_hessian
    189    else:
    190        jacobian = auto_jacobian
    191        hessian = auto_hessian
    192
    193    key_ls = sorted(list(xd.keys()))
    194
    195    if sorted(list(yd.keys())) != key_ls:
    196        raise ValueError('x and y dictionaries do not contain the same keys.')
    197
    198    if sorted(list(funcd.keys())) != key_ls:
    199        raise ValueError('x and func dictionaries do not contain the same keys.')
    200
    201    x_all = np.concatenate([np.array(xd[key]) for key in key_ls])
    202    y_all = np.concatenate([np.array(yd[key]) for key in key_ls])
    203
    204    y_f = [o.value for o in y_all]
    205    dy_f = [o.dvalue for o in y_all]
    206
    207    if len(x_all.shape) > 2:
    208        raise ValueError("Unknown format for x values")
    209
    210    if np.any(np.asarray(dy_f) <= 0.0):
    211        raise Exception("No y errors available, run the gamma method first.")
    212
    213    # number of fit parameters
    214    n_parms_ls = []
    215    for key in key_ls:
    216        if not callable(funcd[key]):
    217            raise TypeError('func (key=' + key + ') is not a function.')
    218        if np.asarray(xd[key]).shape[-1] != len(yd[key]):
    219            raise ValueError('x and y input (key=' + key + ') do not have the same length')
    220        for n_loc in range(100):
    221            try:
    222                funcd[key](np.arange(n_loc), x_all.T[0])
    223            except TypeError:
    224                continue
    225            except IndexError:
    226                continue
    227            else:
    228                break
    229        else:
    230            raise RuntimeError("Fit function (key=" + key + ") is not valid.")
    231        n_parms_ls.append(n_loc)
    232
    233    n_parms = max(n_parms_ls)
    234    if not silent:
    235        print('Fit with', n_parms, 'parameter' + 's' * (n_parms > 1))
    236
    237    if priors is not None:
    238        if isinstance(priors, (list, np.ndarray)):
    239            if n_parms != len(priors):
    240                raise ValueError("'priors' does not have the correct length.")
    241
    242            loc_priors = []
    243            for i_n, i_prior in enumerate(priors):
    244                loc_priors.append(_construct_prior_obs(i_prior, i_n))
    245
    246            prior_mask = np.arange(len(priors))
    247            output.priors = loc_priors
    248
    249        elif isinstance(priors, dict):
    250            loc_priors = []
    251            prior_mask = []
    252            output.priors = {}
    253            for pos, prior in priors.items():
    254                if isinstance(pos, int):
    255                    prior_mask.append(pos)
    256                else:
    257                    raise TypeError("Prior position needs to be an integer.")
    258                loc_priors.append(_construct_prior_obs(prior, pos))
    259
    260                output.priors[pos] = loc_priors[-1]
    261            if max(prior_mask) >= n_parms:
    262                raise ValueError("Prior position out of range.")
    263        else:
    264            raise TypeError("Unkown type for `priors`.")
    265
    266        p_f = [o.value for o in loc_priors]
    267        dp_f = [o.dvalue for o in loc_priors]
    268        if np.any(np.asarray(dp_f) <= 0.0):
    269            raise Exception("No prior errors available, run the gamma method first.")
    270    else:
    271        p_f = dp_f = np.array([])
    272        prior_mask = []
    273        loc_priors = []
    274
    275    if 'initial_guess' in kwargs:
    276        x0 = kwargs.get('initial_guess')
    277        if len(x0) != n_parms:
    278            raise ValueError('Initial guess does not have the correct length: %d vs. %d' % (len(x0), n_parms))
    279    else:
    280        x0 = [0.1] * n_parms
    281
    282    if priors is None:
    283        def general_chisqfunc_uncorr(p, ivars, pr):
    284            model = anp.concatenate([anp.array(funcd[key](p, xd[key])).reshape(-1) for key in key_ls])
    285            return (ivars - model) / dy_f
    286    else:
    287        def general_chisqfunc_uncorr(p, ivars, pr):
    288            model = anp.concatenate([anp.array(funcd[key](p, xd[key])).reshape(-1) for key in key_ls])
    289            return anp.concatenate(((ivars - model) / dy_f, (p[prior_mask] - pr) / dp_f))
    290
    291    def chisqfunc_uncorr(p):
    292        return anp.sum(general_chisqfunc_uncorr(p, y_f, p_f) ** 2)
    293
    294    if kwargs.get('correlated_fit') is True:
    295        corr = covariance(y_all, correlation=True, **kwargs)
    296        covdiag = np.diag(1 / np.asarray(dy_f))
    297        condn = np.linalg.cond(corr)
    298        if condn > 0.1 / np.finfo(float).eps:
    299            raise Exception(f"Cannot invert correlation matrix as its condition number exceeds machine precision ({condn:1.2e})")
    300        if condn > 1e13:
    301            warnings.warn("Correlation matrix may be ill-conditioned, condition number: {%1.2e}" % (condn), RuntimeWarning)
    302        chol = np.linalg.cholesky(corr)
    303        chol_inv = scipy.linalg.solve_triangular(chol, covdiag, lower=True)
    304
    305        def general_chisqfunc(p, ivars, pr):
    306            model = anp.concatenate([anp.array(funcd[key](p, xd[key])).reshape(-1) for key in key_ls])
    307            return anp.concatenate((anp.dot(chol_inv, (ivars - model)), (p[prior_mask] - pr) / dp_f))
    308
    309        def chisqfunc(p):
    310            return anp.sum(general_chisqfunc(p, y_f, p_f) ** 2)
    311    else:
    312        general_chisqfunc = general_chisqfunc_uncorr
    313        chisqfunc = chisqfunc_uncorr
    314
    315    output.method = kwargs.get('method', 'Levenberg-Marquardt')
    316    if not silent:
    317        print('Method:', output.method)
    318
    319    if output.method != 'Levenberg-Marquardt':
    320        if output.method == 'migrad':
    321            tolerance = 1e-4  # default value of 1e-1 set by iminuit can be problematic
    322            if 'tol' in kwargs:
    323                tolerance = kwargs.get('tol')
    324            fit_result = iminuit.minimize(chisqfunc_uncorr, x0, tol=tolerance)  # Stopping criterion 0.002 * tol * errordef
    325            if kwargs.get('correlated_fit') is True:
    326                fit_result = iminuit.minimize(chisqfunc, fit_result.x, tol=tolerance)
    327            output.iterations = fit_result.nfev
    328        else:
    329            tolerance = 1e-12
    330            if 'tol' in kwargs:
    331                tolerance = kwargs.get('tol')
    332            fit_result = scipy.optimize.minimize(chisqfunc_uncorr, x0, method=kwargs.get('method'), tol=tolerance)
    333            if kwargs.get('correlated_fit') is True:
    334                fit_result = scipy.optimize.minimize(chisqfunc, fit_result.x, method=kwargs.get('method'), tol=tolerance)
    335            output.iterations = fit_result.nit
    336
    337        chisquare = fit_result.fun
    338
    339    else:
    340        if 'tol' in kwargs:
    341            print('tol cannot be set for Levenberg-Marquardt')
    342
    343        def chisqfunc_residuals_uncorr(p):
    344            return general_chisqfunc_uncorr(p, y_f, p_f)
    345
    346        fit_result = scipy.optimize.least_squares(chisqfunc_residuals_uncorr, x0, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15)
    347        if kwargs.get('correlated_fit') is True:
    348
    349            def chisqfunc_residuals(p):
    350                return general_chisqfunc(p, y_f, p_f)
    351
    352            fit_result = scipy.optimize.least_squares(chisqfunc_residuals, fit_result.x, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15)
    353
    354        chisquare = np.sum(fit_result.fun ** 2)
    355        assert np.isclose(chisquare, chisqfunc(fit_result.x), atol=1e-14)
    356
    357        output.iterations = fit_result.nfev
    358
    359    if not fit_result.success:
    360        raise Exception('The minimization procedure did not converge.')
    361
    362    output.chisquare = chisquare
    363    output.dof = x_all.shape[-1] - n_parms + len(loc_priors)
    364    output.p_value = 1 - scipy.stats.chi2.cdf(output.chisquare, output.dof)
    365    if output.dof > 0:
    366        output.chisquare_by_dof = output.chisquare / output.dof
    367    else:
    368        output.chisquare_by_dof = float('nan')
    369
    370    output.message = fit_result.message
    371    if not silent:
    372        print(fit_result.message)
    373        print('chisquare/d.o.f.:', output.chisquare_by_dof)
    374        print('fit parameters', fit_result.x)
    375
    376    def prepare_hat_matrix():
    377        hat_vector = []
    378        for key in key_ls:
    379            if (len(xd[key]) != 0):
    380                hat_vector.append(jacobian(funcd[key])(fit_result.x, xd[key]))
    381        hat_vector = [item for sublist in hat_vector for item in sublist]
    382        return hat_vector
    383
    384    if kwargs.get('expected_chisquare') is True:
    385        if kwargs.get('correlated_fit') is not True:
    386            W = np.diag(1 / np.asarray(dy_f))
    387            cov = covariance(y_all)
    388            hat_vector = prepare_hat_matrix()
    389            A = W @ hat_vector
    390            P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T
    391            expected_chisquare = np.trace((np.identity(x_all.shape[-1]) - P_phi) @ W @ cov @ W)
    392            output.chisquare_by_expected_chisquare = output.chisquare / expected_chisquare
    393            if not silent:
    394                print('chisquare/expected_chisquare:', output.chisquare_by_expected_chisquare)
    395
    396    fitp = fit_result.x
    397
    398    try:
    399        hess = hessian(chisqfunc)(fitp)
    400    except TypeError:
    401        raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None
    402
    403    len_y = len(y_f)
    404
    405    def chisqfunc_compact(d):
    406        return anp.sum(general_chisqfunc(d[:n_parms], d[n_parms: n_parms + len_y], d[n_parms + len_y:]) ** 2)
    407
    408    jac_jac_y = hessian(chisqfunc_compact)(np.concatenate((fitp, y_f, p_f)))
    409
    410    # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv
    411    try:
    412        deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms, n_parms:])
    413    except np.linalg.LinAlgError:
    414        raise Exception("Cannot invert hessian matrix.")
    415
    416    result = []
    417    for i in range(n_parms):
    418        result.append(derived_observable(lambda x_all, **kwargs: (x_all[0] + np.finfo(np.float64).eps) / (y_all[0].value + np.finfo(np.float64).eps) * fitp[i], list(y_all) + loc_priors, man_grad=list(deriv_y[i])))
    419
    420    output.fit_parameters = result
    421
    422    # Hotelling t-squared p-value for correlated fits.
    423    if kwargs.get('correlated_fit') is True:
    424        n_cov = np.min(np.vectorize(lambda x_all: x_all.N)(y_all))
    425        output.t2_p_value = 1 - scipy.stats.f.cdf((n_cov - output.dof) / (output.dof * (n_cov - 1)) * output.chisquare,
    426                                                  output.dof, n_cov - output.dof)
    427
    428    if kwargs.get('resplot') is True:
    429        for key in key_ls:
    430            residual_plot(xd[key], yd[key], funcd[key], result, title=key)
    431
    432    if kwargs.get('qqplot') is True:
    433        for key in key_ls:
    434            qqplot(xd[key], yd[key], funcd[key], result, title=key)
    435
    436    return output

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

###### Parameters

- For an uncombined fit:
- 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.

- OR For a combined fit:
- x (dict): dict of lists.
- y (dict): dict of lists of Obs.

- funcs (dict): dict of objects fit functions have to be of the form (here a[0] is the common fit parameter) ```python import autograd.numpy as anp funcs = {"a": func_a, "b": func_b}

  def func_a(a, x): return a[1] * anp.exp(-a[0] * x)

  def func_b(a, x): return a[2] * anp.exp(-a[0] * x)

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

- priors (dict or list, optional): priors can either be a dictionary with integer keys and the corresponding priors as values or 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.
- tol (float, optional): can be used (only for combined fits and methods other than Levenberg-Marquard) to set the tolerance for convergence to a different value to either speed up convergence at the cost of a larger error on the fitted parameters (and possibly invalid estimates for parameter uncertainties) or smaller values to get more accurate parameter values The stopping criterion depends on the method, e.g. migrad: edm_max = 0.002 * tol * errordef (EDM criterion: edm < edm_max)
- 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).
- 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).
- num_grad (bool): Use numerical differentation instead of automatic differentiation to perform the error propagation (default False).

###### Returns

- output (Fit_result): Parameters and information on the fitted result.

def total_least_squares(x, y, func, silent=False, **kwargs): View Source

    439def total_least_squares(x, y, func, silent=False, **kwargs):
    440    r'''Performs a non-linear fit to y = func(x) and returns a list of Obs corresponding to the fit parameters.
    441
    442    Parameters
    443    ----------
    444    x : list
    445        list of Obs, or a tuple of lists of Obs
    446    y : list
    447        list of Obs. The dvalues of the Obs are used as x- and yerror for the fit.
    448    func : object
    449        func has to be of the form
    450
    451        ```python
    452        import autograd.numpy as anp
    453
    454        def func(a, x):
    455            return a[0] + a[1] * x + a[2] * anp.sinh(x)
    456        ```
    457
    458        For multiple x values func can be of the form
    459
    460        ```python
    461        def func(a, x):
    462            (x1, x2) = x
    463            return a[0] * x1 ** 2 + a[1] * x2
    464        ```
    465
    466        It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation
    467        will not work.
    468    silent : bool, optional
    469        If true all output to the console is omitted (default False).
    470    initial_guess : list
    471        can provide an initial guess for the input parameters. Relevant for non-linear
    472        fits with many parameters.
    473    expected_chisquare : bool
    474        If true prints the expected chisquare which is
    475        corrected by effects caused by correlated input data.
    476        This can take a while as the full correlation matrix
    477        has to be calculated (default False).
    478    num_grad : bool
    479        Use numerical differentation instead of automatic differentiation to perform the error propagation (default False).
    480
    481    Notes
    482    -----
    483    Based on the orthogonal distance regression module of scipy.
    484
    485    Returns
    486    -------
    487    output : Fit_result
    488        Parameters and information on the fitted result.
    489    '''
    490
    491    output = Fit_result()
    492
    493    output.fit_function = func
    494
    495    x = np.array(x)
    496
    497    x_shape = x.shape
    498
    499    if kwargs.get('num_grad') is True:
    500        jacobian = num_jacobian
    501        hessian = num_hessian
    502    else:
    503        jacobian = auto_jacobian
    504        hessian = auto_hessian
    505
    506    if not callable(func):
    507        raise TypeError('func has to be a function.')
    508
    509    for i in range(42):
    510        try:
    511            func(np.arange(i), x.T[0])
    512        except TypeError:
    513            continue
    514        except IndexError:
    515            continue
    516        else:
    517            break
    518    else:
    519        raise RuntimeError("Fit function is not valid.")
    520
    521    n_parms = i
    522    if not silent:
    523        print('Fit with', n_parms, 'parameter' + 's' * (n_parms > 1))
    524
    525    x_f = np.vectorize(lambda o: o.value)(x)
    526    dx_f = np.vectorize(lambda o: o.dvalue)(x)
    527    y_f = np.array([o.value for o in y])
    528    dy_f = np.array([o.dvalue for o in y])
    529
    530    if np.any(np.asarray(dx_f) <= 0.0):
    531        raise Exception('No x errors available, run the gamma method first.')
    532
    533    if np.any(np.asarray(dy_f) <= 0.0):
    534        raise Exception('No y errors available, run the gamma method first.')
    535
    536    if 'initial_guess' in kwargs:
    537        x0 = kwargs.get('initial_guess')
    538        if len(x0) != n_parms:
    539            raise Exception('Initial guess does not have the correct length: %d vs. %d' % (len(x0), n_parms))
    540    else:
    541        x0 = [1] * n_parms
    542
    543    data = RealData(x_f, y_f, sx=dx_f, sy=dy_f)
    544    model = Model(func)
    545    odr = ODR(data, model, x0, partol=np.finfo(np.float64).eps)
    546    odr.set_job(fit_type=0, deriv=1)
    547    out = odr.run()
    548
    549    output.residual_variance = out.res_var
    550
    551    output.method = 'ODR'
    552
    553    output.message = out.stopreason
    554
    555    output.xplus = out.xplus
    556
    557    if not silent:
    558        print('Method: ODR')
    559        print(*out.stopreason)
    560        print('Residual variance:', output.residual_variance)
    561
    562    if out.info > 3:
    563        raise Exception('The minimization procedure did not converge.')
    564
    565    m = x_f.size
    566
    567    def odr_chisquare(p):
    568        model = func(p[:n_parms], p[n_parms:].reshape(x_shape))
    569        chisq = anp.sum(((y_f - model) / dy_f) ** 2) + anp.sum(((x_f - p[n_parms:].reshape(x_shape)) / dx_f) ** 2)
    570        return chisq
    571
    572    if kwargs.get('expected_chisquare') is True:
    573        W = np.diag(1 / np.asarray(np.concatenate((dy_f.ravel(), dx_f.ravel()))))
    574
    575        if kwargs.get('covariance') is not None:
    576            cov = kwargs.get('covariance')
    577        else:
    578            cov = covariance(np.concatenate((y, x.ravel())))
    579
    580        number_of_x_parameters = int(m / x_f.shape[-1])
    581
    582        old_jac = jacobian(func)(out.beta, out.xplus)
    583        fused_row1 = np.concatenate((old_jac, np.concatenate((number_of_x_parameters * [np.zeros(old_jac.shape)]), axis=0)))
    584        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])))
    585        new_jac = np.concatenate((fused_row1, fused_row2), axis=1)
    586
    587        A = W @ new_jac
    588        P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T
    589        expected_chisquare = np.trace((np.identity(P_phi.shape[0]) - P_phi) @ W @ cov @ W)
    590        if expected_chisquare <= 0.0:
    591            warnings.warn("Negative expected_chisquare.", RuntimeWarning)
    592            expected_chisquare = np.abs(expected_chisquare)
    593        output.chisquare_by_expected_chisquare = odr_chisquare(np.concatenate((out.beta, out.xplus.ravel()))) / expected_chisquare
    594        if not silent:
    595            print('chisquare/expected_chisquare:',
    596                  output.chisquare_by_expected_chisquare)
    597
    598    fitp = out.beta
    599    try:
    600        hess = hessian(odr_chisquare)(np.concatenate((fitp, out.xplus.ravel())))
    601    except TypeError:
    602        raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None
    603
    604    def odr_chisquare_compact_x(d):
    605        model = func(d[:n_parms], d[n_parms:n_parms + m].reshape(x_shape))
    606        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)
    607        return chisq
    608
    609    jac_jac_x = hessian(odr_chisquare_compact_x)(np.concatenate((fitp, out.xplus.ravel(), x_f.ravel())))
    610
    611    # Compute hess^{-1} @ jac_jac_x[:n_parms + m, n_parms + m:] using LAPACK dgesv
    612    try:
    613        deriv_x = -scipy.linalg.solve(hess, jac_jac_x[:n_parms + m, n_parms + m:])
    614    except np.linalg.LinAlgError:
    615        raise Exception("Cannot invert hessian matrix.")
    616
    617    def odr_chisquare_compact_y(d):
    618        model = func(d[:n_parms], d[n_parms:n_parms + m].reshape(x_shape))
    619        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)
    620        return chisq
    621
    622    jac_jac_y = hessian(odr_chisquare_compact_y)(np.concatenate((fitp, out.xplus.ravel(), y_f)))
    623
    624    # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv
    625    try:
    626        deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms + m, n_parms + m:])
    627    except np.linalg.LinAlgError:
    628        raise Exception("Cannot invert hessian matrix.")
    629
    630    result = []
    631    for i in range(n_parms):
    632        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])))
    633
    634    output.fit_parameters = result
    635
    636    output.odr_chisquare = odr_chisquare(np.concatenate((out.beta, out.xplus.ravel())))
    637    output.dof = x.shape[-1] - n_parms
    638    output.p_value = 1 - scipy.stats.chi2.cdf(output.odr_chisquare, output.dof)
    639
    640    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).
- num_grad (bool): Use numerical differentation instead of automatic differentiation to perform the error propagation (default False).

###### Notes

Based on the orthogonal distance regression module of scipy.

###### Returns

- output (Fit_result): Parameters and information on the fitted result.

def fit_lin(x, y, **kwargs): View Source

    643def fit_lin(x, y, **kwargs):
    644    """Performs a linear fit to y = n + m * x and returns two Obs n, m.
    645
    646    Parameters
    647    ----------
    648    x : list
    649        Can either be a list of floats in which case no xerror is assumed, or
    650        a list of Obs, where the dvalues of the Obs are used as xerror for the fit.
    651    y : list
    652        List of Obs, the dvalues of the Obs are used as yerror for the fit.
    653
    654    Returns
    655    -------
    656    fit_parameters : list[Obs]
    657        LIist of fitted observables.
    658    """
    659
    660    def f(a, x):
    661        y = a[0] + a[1] * x
    662        return y
    663
    664    if all(isinstance(n, Obs) for n in x):
    665        out = total_least_squares(x, y, f, **kwargs)
    666        return out.fit_parameters
    667    elif all(isinstance(n, float) or isinstance(n, int) for n in x) or isinstance(x, np.ndarray):
    668        out = least_squares(x, y, f, **kwargs)
    669        return out.fit_parameters
    670    else:
    671        raise TypeError('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.

###### Returns

- fit_parameters (list[Obs]): LIist of fitted observables.

def qqplot(x, o_y, func, p, title=''): View Source

    674def qqplot(x, o_y, func, p, title=""):
    675    """Generates a quantile-quantile plot of the fit result which can be used to
    676       check if the residuals of the fit are gaussian distributed.
    677
    678    Returns
    679    -------
    680    None
    681    """
    682
    683    residuals = []
    684    for i_x, i_y in zip(x, o_y):
    685        residuals.append((i_y - func(p, i_x)) / i_y.dvalue)
    686    residuals = sorted(residuals)
    687    my_y = [o.value for o in residuals]
    688    probplot = scipy.stats.probplot(my_y)
    689    my_x = probplot[0][0]
    690    plt.figure(figsize=(8, 8 / 1.618))
    691    plt.errorbar(my_x, my_y, fmt='o')
    692    fit_start = my_x[0]
    693    fit_stop = my_x[-1]
    694    samples = np.arange(fit_start, fit_stop, 0.01)
    695    plt.plot(samples, samples, 'k--', zorder=11, label='Standard normal distribution')
    696    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='-')
    697
    698    plt.xlabel('Theoretical quantiles')
    699    plt.ylabel('Ordered Values')
    700    plt.legend(title=title)
    701    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.

###### Returns

- None

def residual_plot(x, y, func, fit_res, title=''): View Source

    704def residual_plot(x, y, func, fit_res, title=""):
    705    """Generates a plot which compares the fit to the data and displays the corresponding residuals
    706
    707    For uncorrelated data the residuals are expected to be distributed ~N(0,1).
    708
    709    Returns
    710    -------
    711    None
    712    """
    713    sorted_x = sorted(x)
    714    xstart = sorted_x[0] - 0.5 * (sorted_x[1] - sorted_x[0])
    715    xstop = sorted_x[-1] + 0.5 * (sorted_x[-1] - sorted_x[-2])
    716    x_samples = np.arange(xstart, xstop + 0.01, 0.01)
    717
    718    plt.figure(figsize=(8, 8 / 1.618))
    719    gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
    720    ax0 = plt.subplot(gs[0])
    721    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')
    722    ax0.plot(x_samples, func([o.value for o in fit_res], x_samples), label='Fit', zorder=10, ls='-', ms=0)
    723    ax0.set_xticklabels([])
    724    ax0.set_xlim([xstart, xstop])
    725    ax0.set_xticklabels([])
    726    ax0.legend(title=title)
    727
    728    residuals = (np.asarray([o.value for o in y]) - func([o.value for o in fit_res], np.asarray(x))) / np.asarray([o.dvalue for o in y])
    729    ax1 = plt.subplot(gs[1])
    730    ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
    731    ax1.tick_params(direction='out')
    732    ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
    733    ax1.axhline(y=0.0, ls='--', color='k', marker=" ")
    734    ax1.fill_between(x_samples, -1.0, 1.0, alpha=0.1, facecolor='k')
    735    ax1.set_xlim([xstart, xstop])
    736    ax1.set_ylabel('Residuals')
    737    plt.subplots_adjust(wspace=None, hspace=None)
    738    plt.draw()

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

For uncorrelated data the residuals are expected to be distributed ~N(0,1).

###### Returns

- None

def error_band(x, func, beta): View Source

    741def error_band(x, func, beta):
    742    """Calculate the error band for an array of sample values x, for given fit function func with optimized parameters beta.
    743
    744    Returns
    745    -------
    746    err : np.array(Obs)
    747        Error band for an array of sample values x
    748    """
    749    cov = covariance(beta)
    750    if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps):
    751        warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning)
    752
    753    deriv = []
    754    for i, item in enumerate(x):
    755        deriv.append(np.array(egrad(func)([o.value for o in beta], item)))
    756
    757    err = []
    758    for i, item in enumerate(x):
    759        err.append(np.sqrt(deriv[i] @ cov @ deriv[i]))
    760    err = np.array(err)
    761
    762    return err

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

###### Returns

- err (np.array(Obs)): Error band for an array of sample values x

def ks_test(objects=None): View Source

    765def ks_test(objects=None):
    766    """Performs a Kolmogorov–Smirnov test for the p-values of all fit object.
    767
    768    Parameters
    769    ----------
    770    objects : list
    771        List of fit results to include in the analysis (optional).
    772
    773    Returns
    774    -------
    775    None
    776    """
    777
    778    if objects is None:
    779        obs_list = []
    780        for obj in gc.get_objects():
    781            if isinstance(obj, Fit_result):
    782                obs_list.append(obj)
    783    else:
    784        obs_list = objects
    785
    786    p_values = [o.p_value for o in obs_list]
    787
    788    bins = len(p_values)
    789    x = np.arange(0, 1.001, 0.001)
    790    plt.plot(x, x, 'k', zorder=1)
    791    plt.xlim(0, 1)
    792    plt.ylim(0, 1)
    793    plt.xlabel('p-value')
    794    plt.ylabel('Cumulative probability')
    795    plt.title(str(bins) + ' p-values')
    796
    797    n = np.arange(1, bins + 1) / np.float64(bins)
    798    Xs = np.sort(p_values)
    799    plt.step(Xs, n)
    800    diffs = n - Xs
    801    loc_max_diff = np.argmax(np.abs(diffs))
    802    loc = Xs[loc_max_diff]
    803    plt.annotate('', xy=(loc, loc), xytext=(loc, loc + diffs[loc_max_diff]), arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
    804    plt.draw()
    805
    806    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).

###### Returns

- None