# 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 i in range(100):
    220            try:
    221                funcd[key](np.arange(i), 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 = i
    231        n_parms_ls.append(n_parms)
    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                if isinstance(i_prior, Obs):
    244                    loc_priors.append(i_prior)
    245                elif isinstance(i_prior, str):
    246                    loc_val, loc_dval = _extract_val_and_dval(i_prior)
    247                    loc_priors.append(cov_Obs(loc_val, loc_dval ** 2, '#prior' + str(i_n) + f"_{np.random.randint(2147483647):010d}"))
    248                else:
    249                    raise TypeError("Prior entries need to be 'Obs' or 'str'.")
    250
    251            prior_mask = np.arange(len(priors))
    252            output.priors = loc_priors
    253
    254        elif isinstance(priors, dict):
    255            loc_priors = []
    256            prior_mask = []
    257            output.priors = {}
    258            for pos, prior in priors.items():
    259                if isinstance(pos, int):
    260                    prior_mask.append(pos)
    261                else:
    262                    raise TypeError("Prior position needs to be an integer.")
    263                if isinstance(prior, Obs):
    264                    loc_priors.append(prior)
    265                elif isinstance(prior, str):
    266                    loc_val, loc_dval = _extract_val_and_dval(prior)
    267                    loc_priors.append(cov_Obs(loc_val, loc_dval ** 2, '#prior' + str(pos) + f"_{np.random.randint(2147483647):010d}"))
    268                else:
    269                    raise TypeError("Prior entries need to be 'Obs' or 'str'.")
    270                output.priors[pos] = loc_priors[-1]
    271            if max(prior_mask) >= n_parms:
    272                raise ValueError("Prior position out of range.")
    273        else:
    274            raise TypeError("Unkown type for `priors`.")
    275
    276        p_f = [o.value for o in loc_priors]
    277        dp_f = [o.dvalue for o in loc_priors]
    278        if np.any(np.asarray(dp_f) <= 0.0):
    279            raise Exception("No prior errors available, run the gamma method first.")
    280    else:
    281        p_f = dp_f = np.array([])
    282        prior_mask = []
    283        loc_priors = []
    284
    285    if 'initial_guess' in kwargs:
    286        x0 = kwargs.get('initial_guess')
    287        if len(x0) != n_parms:
    288            raise ValueError('Initial guess does not have the correct length: %d vs. %d' % (len(x0), n_parms))
    289    else:
    290        x0 = [0.1] * n_parms
    291
    292    if priors is None:
    293        def general_chisqfunc_uncorr(p, ivars, pr):
    294            model = anp.concatenate([anp.array(funcd[key](p, xd[key])).reshape(-1) for key in key_ls])
    295            return (ivars - model) / dy_f
    296    else:
    297        def general_chisqfunc_uncorr(p, ivars, pr):
    298            model = anp.concatenate([anp.array(funcd[key](p, xd[key])).reshape(-1) for key in key_ls])
    299            return anp.concatenate(((ivars - model) / dy_f, (p[prior_mask] - pr) / dp_f))
    300
    301    def chisqfunc_uncorr(p):
    302        return anp.sum(general_chisqfunc_uncorr(p, y_f, p_f) ** 2)
    303
    304    if kwargs.get('correlated_fit') is True:
    305        corr = covariance(y_all, correlation=True, **kwargs)
    306        covdiag = np.diag(1 / np.asarray(dy_f))
    307        condn = np.linalg.cond(corr)
    308        if condn > 0.1 / np.finfo(float).eps:
    309            raise Exception(f"Cannot invert correlation matrix as its condition number exceeds machine precision ({condn:1.2e})")
    310        if condn > 1e13:
    311            warnings.warn("Correlation matrix may be ill-conditioned, condition number: {%1.2e}" % (condn), RuntimeWarning)
    312        chol = np.linalg.cholesky(corr)
    313        chol_inv = scipy.linalg.solve_triangular(chol, covdiag, lower=True)
    314
    315        def general_chisqfunc(p, ivars, pr):
    316            model = anp.concatenate([anp.array(funcd[key](p, xd[key])).reshape(-1) for key in key_ls])
    317            return anp.concatenate((anp.dot(chol_inv, (ivars - model)), (p[prior_mask] - pr) / dp_f))
    318
    319        def chisqfunc(p):
    320            return anp.sum(general_chisqfunc(p, y_f, p_f) ** 2)
    321    else:
    322        general_chisqfunc = general_chisqfunc_uncorr
    323        chisqfunc = chisqfunc_uncorr
    324
    325    output.method = kwargs.get('method', 'Levenberg-Marquardt')
    326    if not silent:
    327        print('Method:', output.method)
    328
    329    if output.method != 'Levenberg-Marquardt':
    330        if output.method == 'migrad':
    331            tolerance = 1e-4  # default value of 1e-1 set by iminuit can be problematic
    332            if 'tol' in kwargs:
    333                tolerance = kwargs.get('tol')
    334            fit_result = iminuit.minimize(chisqfunc_uncorr, x0, tol=tolerance)  # Stopping criterion 0.002 * tol * errordef
    335            if kwargs.get('correlated_fit') is True:
    336                fit_result = iminuit.minimize(chisqfunc, fit_result.x, tol=tolerance)
    337            output.iterations = fit_result.nfev
    338        else:
    339            tolerance = 1e-12
    340            if 'tol' in kwargs:
    341                tolerance = kwargs.get('tol')
    342            fit_result = scipy.optimize.minimize(chisqfunc_uncorr, x0, method=kwargs.get('method'), tol=tolerance)
    343            if kwargs.get('correlated_fit') is True:
    344                fit_result = scipy.optimize.minimize(chisqfunc, fit_result.x, method=kwargs.get('method'), tol=tolerance)
    345            output.iterations = fit_result.nit
    346
    347        chisquare = fit_result.fun
    348
    349    else:
    350        if 'tol' in kwargs:
    351            print('tol cannot be set for Levenberg-Marquardt')
    352
    353        def chisqfunc_residuals_uncorr(p):
    354            return general_chisqfunc_uncorr(p, y_f, p_f)
    355
    356        fit_result = scipy.optimize.least_squares(chisqfunc_residuals_uncorr, x0, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15)
    357        if kwargs.get('correlated_fit') is True:
    358
    359            def chisqfunc_residuals(p):
    360                return general_chisqfunc(p, y_f, p_f)
    361
    362            fit_result = scipy.optimize.least_squares(chisqfunc_residuals, fit_result.x, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15)
    363
    364        chisquare = np.sum(fit_result.fun ** 2)
    365        assert np.isclose(chisquare, chisqfunc(fit_result.x), atol=1e-14)
    366
    367        output.iterations = fit_result.nfev
    368
    369    if not fit_result.success:
    370        raise Exception('The minimization procedure did not converge.')
    371
    372    if x_all.shape[-1] - n_parms > 0:
    373        output.chisquare = chisquare
    374        output.dof = x_all.shape[-1] - n_parms + len(loc_priors)
    375        output.chisquare_by_dof = output.chisquare / output.dof
    376        output.p_value = 1 - scipy.stats.chi2.cdf(output.chisquare, output.dof)
    377    else:
    378        output.chisquare_by_dof = float('nan')
    379
    380    output.message = fit_result.message
    381    if not silent:
    382        print(fit_result.message)
    383        print('chisquare/d.o.f.:', output.chisquare_by_dof)
    384        print('fit parameters', fit_result.x)
    385
    386    def prepare_hat_matrix():
    387        hat_vector = []
    388        for key in key_ls:
    389            if (len(xd[key]) != 0):
    390                hat_vector.append(jacobian(funcd[key])(fit_result.x, xd[key]))
    391        hat_vector = [item for sublist in hat_vector for item in sublist]
    392        return hat_vector
    393
    394    if kwargs.get('expected_chisquare') is True:
    395        if kwargs.get('correlated_fit') is not True:
    396            W = np.diag(1 / np.asarray(dy_f))
    397            cov = covariance(y_all)
    398            hat_vector = prepare_hat_matrix()
    399            A = W @ hat_vector
    400            P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T
    401            expected_chisquare = np.trace((np.identity(x_all.shape[-1]) - P_phi) @ W @ cov @ W)
    402            output.chisquare_by_expected_chisquare = output.chisquare / expected_chisquare
    403            if not silent:
    404                print('chisquare/expected_chisquare:', output.chisquare_by_expected_chisquare)
    405
    406    fitp = fit_result.x
    407
    408    try:
    409        hess = hessian(chisqfunc)(fitp)
    410    except TypeError:
    411        raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None
    412
    413    len_y = len(y_f)
    414
    415    def chisqfunc_compact(d):
    416        return anp.sum(general_chisqfunc(d[:n_parms], d[n_parms: n_parms + len_y], d[n_parms + len_y:]) ** 2)
    417
    418    jac_jac_y = hessian(chisqfunc_compact)(np.concatenate((fitp, y_f, p_f)))
    419
    420    # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv
    421    try:
    422        deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms, n_parms:])
    423    except np.linalg.LinAlgError:
    424        raise Exception("Cannot invert hessian matrix.")
    425
    426    result = []
    427    for i in range(n_parms):
    428        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])))
    429
    430    output.fit_parameters = result
    431
    432    # Hotelling t-squared p-value for correlated fits.
    433    if kwargs.get('correlated_fit') is True:
    434        n_cov = np.min(np.vectorize(lambda x_all: x_all.N)(y_all))
    435        output.t2_p_value = 1 - scipy.stats.f.cdf((n_cov - output.dof) / (output.dof * (n_cov - 1)) * output.chisquare,
    436                                                  output.dof, n_cov - output.dof)
    437
    438    if kwargs.get('resplot') is True:
    439        for key in key_ls:
    440            residual_plot(xd[key], yd[key], funcd[key], result, title=key)
    441
    442    if kwargs.get('qqplot') is True:
    443        for key in key_ls:
    444            qqplot(xd[key], yd[key], funcd[key], result, title=key)
    445
    446    return output
    447
    448
    449def total_least_squares(x, y, func, silent=False, **kwargs):
    450    r'''Performs a non-linear fit to y = func(x) and returns a list of Obs corresponding to the fit parameters.
    451
    452    Parameters
    453    ----------
    454    x : list
    455        list of Obs, or a tuple of lists of Obs
    456    y : list
    457        list of Obs. The dvalues of the Obs are used as x- and yerror for the fit.
    458    func : object
    459        func has to be of the form
    460
    461        ```python
    462        import autograd.numpy as anp
    463
    464        def func(a, x):
    465            return a[0] + a[1] * x + a[2] * anp.sinh(x)
    466        ```
    467
    468        For multiple x values func can be of the form
    469
    470        ```python
    471        def func(a, x):
    472            (x1, x2) = x
    473            return a[0] * x1 ** 2 + a[1] * x2
    474        ```
    475
    476        It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation
    477        will not work.
    478    silent : bool, optional
    479        If true all output to the console is omitted (default False).
    480    initial_guess : list
    481        can provide an initial guess for the input parameters. Relevant for non-linear
    482        fits with many parameters.
    483    expected_chisquare : bool
    484        If true prints the expected chisquare which is
    485        corrected by effects caused by correlated input data.
    486        This can take a while as the full correlation matrix
    487        has to be calculated (default False).
    488    num_grad : bool
    489        Use numerical differentation instead of automatic differentiation to perform the error propagation (default False).
    490
    491    Notes
    492    -----
    493    Based on the orthogonal distance regression module of scipy.
    494
    495    Returns
    496    -------
    497    output : Fit_result
    498        Parameters and information on the fitted result.
    499    '''
    500
    501    output = Fit_result()
    502
    503    output.fit_function = func
    504
    505    x = np.array(x)
    506
    507    x_shape = x.shape
    508
    509    if kwargs.get('num_grad') is True:
    510        jacobian = num_jacobian
    511        hessian = num_hessian
    512    else:
    513        jacobian = auto_jacobian
    514        hessian = auto_hessian
    515
    516    if not callable(func):
    517        raise TypeError('func has to be a function.')
    518
    519    for i in range(42):
    520        try:
    521            func(np.arange(i), x.T[0])
    522        except TypeError:
    523            continue
    524        except IndexError:
    525            continue
    526        else:
    527            break
    528    else:
    529        raise RuntimeError("Fit function is not valid.")
    530
    531    n_parms = i
    532    if not silent:
    533        print('Fit with', n_parms, 'parameter' + 's' * (n_parms > 1))
    534
    535    x_f = np.vectorize(lambda o: o.value)(x)
    536    dx_f = np.vectorize(lambda o: o.dvalue)(x)
    537    y_f = np.array([o.value for o in y])
    538    dy_f = np.array([o.dvalue for o in y])
    539
    540    if np.any(np.asarray(dx_f) <= 0.0):
    541        raise Exception('No x errors available, run the gamma method first.')
    542
    543    if np.any(np.asarray(dy_f) <= 0.0):
    544        raise Exception('No y errors available, run the gamma method first.')
    545
    546    if 'initial_guess' in kwargs:
    547        x0 = kwargs.get('initial_guess')
    548        if len(x0) != n_parms:
    549            raise Exception('Initial guess does not have the correct length: %d vs. %d' % (len(x0), n_parms))
    550    else:
    551        x0 = [1] * n_parms
    552
    553    data = RealData(x_f, y_f, sx=dx_f, sy=dy_f)
    554    model = Model(func)
    555    odr = ODR(data, model, x0, partol=np.finfo(np.float64).eps)
    556    odr.set_job(fit_type=0, deriv=1)
    557    out = odr.run()
    558
    559    output.residual_variance = out.res_var
    560
    561    output.method = 'ODR'
    562
    563    output.message = out.stopreason
    564
    565    output.xplus = out.xplus
    566
    567    if not silent:
    568        print('Method: ODR')
    569        print(*out.stopreason)
    570        print('Residual variance:', output.residual_variance)
    571
    572    if out.info > 3:
    573        raise Exception('The minimization procedure did not converge.')
    574
    575    m = x_f.size
    576
    577    def odr_chisquare(p):
    578        model = func(p[:n_parms], p[n_parms:].reshape(x_shape))
    579        chisq = anp.sum(((y_f - model) / dy_f) ** 2) + anp.sum(((x_f - p[n_parms:].reshape(x_shape)) / dx_f) ** 2)
    580        return chisq
    581
    582    if kwargs.get('expected_chisquare') is True:
    583        W = np.diag(1 / np.asarray(np.concatenate((dy_f.ravel(), dx_f.ravel()))))
    584
    585        if kwargs.get('covariance') is not None:
    586            cov = kwargs.get('covariance')
    587        else:
    588            cov = covariance(np.concatenate((y, x.ravel())))
    589
    590        number_of_x_parameters = int(m / x_f.shape[-1])
    591
    592        old_jac = jacobian(func)(out.beta, out.xplus)
    593        fused_row1 = np.concatenate((old_jac, np.concatenate((number_of_x_parameters * [np.zeros(old_jac.shape)]), axis=0)))
    594        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])))
    595        new_jac = np.concatenate((fused_row1, fused_row2), axis=1)
    596
    597        A = W @ new_jac
    598        P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T
    599        expected_chisquare = np.trace((np.identity(P_phi.shape[0]) - P_phi) @ W @ cov @ W)
    600        if expected_chisquare <= 0.0:
    601            warnings.warn("Negative expected_chisquare.", RuntimeWarning)
    602            expected_chisquare = np.abs(expected_chisquare)
    603        output.chisquare_by_expected_chisquare = odr_chisquare(np.concatenate((out.beta, out.xplus.ravel()))) / expected_chisquare
    604        if not silent:
    605            print('chisquare/expected_chisquare:',
    606                  output.chisquare_by_expected_chisquare)
    607
    608    fitp = out.beta
    609    try:
    610        hess = hessian(odr_chisquare)(np.concatenate((fitp, out.xplus.ravel())))
    611    except TypeError:
    612        raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None
    613
    614    def odr_chisquare_compact_x(d):
    615        model = func(d[:n_parms], d[n_parms:n_parms + m].reshape(x_shape))
    616        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)
    617        return chisq
    618
    619    jac_jac_x = hessian(odr_chisquare_compact_x)(np.concatenate((fitp, out.xplus.ravel(), x_f.ravel())))
    620
    621    # Compute hess^{-1} @ jac_jac_x[:n_parms + m, n_parms + m:] using LAPACK dgesv
    622    try:
    623        deriv_x = -scipy.linalg.solve(hess, jac_jac_x[:n_parms + m, n_parms + m:])
    624    except np.linalg.LinAlgError:
    625        raise Exception("Cannot invert hessian matrix.")
    626
    627    def odr_chisquare_compact_y(d):
    628        model = func(d[:n_parms], d[n_parms:n_parms + m].reshape(x_shape))
    629        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)
    630        return chisq
    631
    632    jac_jac_y = hessian(odr_chisquare_compact_y)(np.concatenate((fitp, out.xplus.ravel(), y_f)))
    633
    634    # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv
    635    try:
    636        deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms + m, n_parms + m:])
    637    except np.linalg.LinAlgError:
    638        raise Exception("Cannot invert hessian matrix.")
    639
    640    result = []
    641    for i in range(n_parms):
    642        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])))
    643
    644    output.fit_parameters = result
    645
    646    output.odr_chisquare = odr_chisquare(np.concatenate((out.beta, out.xplus.ravel())))
    647    output.dof = x.shape[-1] - n_parms
    648    output.p_value = 1 - scipy.stats.chi2.cdf(output.odr_chisquare, output.dof)
    649
    650    return output
    651
    652
    653def fit_lin(x, y, **kwargs):
    654    """Performs a linear fit to y = n + m * x and returns two Obs n, m.
    655
    656    Parameters
    657    ----------
    658    x : list
    659        Can either be a list of floats in which case no xerror is assumed, or
    660        a list of Obs, where the dvalues of the Obs are used as xerror for the fit.
    661    y : list
    662        List of Obs, the dvalues of the Obs are used as yerror for the fit.
    663
    664    Returns
    665    -------
    666    fit_parameters : list[Obs]
    667        LIist of fitted observables.
    668    """
    669
    670    def f(a, x):
    671        y = a[0] + a[1] * x
    672        return y
    673
    674    if all(isinstance(n, Obs) for n in x):
    675        out = total_least_squares(x, y, f, **kwargs)
    676        return out.fit_parameters
    677    elif all(isinstance(n, float) or isinstance(n, int) for n in x) or isinstance(x, np.ndarray):
    678        out = least_squares(x, y, f, **kwargs)
    679        return out.fit_parameters
    680    else:
    681        raise TypeError('Unsupported types for x')
    682
    683
    684def qqplot(x, o_y, func, p, title=""):
    685    """Generates a quantile-quantile plot of the fit result which can be used to
    686       check if the residuals of the fit are gaussian distributed.
    687
    688    Returns
    689    -------
    690    None
    691    """
    692
    693    residuals = []
    694    for i_x, i_y in zip(x, o_y):
    695        residuals.append((i_y - func(p, i_x)) / i_y.dvalue)
    696    residuals = sorted(residuals)
    697    my_y = [o.value for o in residuals]
    698    probplot = scipy.stats.probplot(my_y)
    699    my_x = probplot[0][0]
    700    plt.figure(figsize=(8, 8 / 1.618))
    701    plt.errorbar(my_x, my_y, fmt='o')
    702    fit_start = my_x[0]
    703    fit_stop = my_x[-1]
    704    samples = np.arange(fit_start, fit_stop, 0.01)
    705    plt.plot(samples, samples, 'k--', zorder=11, label='Standard normal distribution')
    706    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='-')
    707
    708    plt.xlabel('Theoretical quantiles')
    709    plt.ylabel('Ordered Values')
    710    plt.legend(title=title)
    711    plt.draw()
    712
    713
    714def residual_plot(x, y, func, fit_res, title=""):
    715    """Generates a plot which compares the fit to the data and displays the corresponding residuals
    716
    717    For uncorrelated data the residuals are expected to be distributed ~N(0,1).
    718
    719    Returns
    720    -------
    721    None
    722    """
    723    sorted_x = sorted(x)
    724    xstart = sorted_x[0] - 0.5 * (sorted_x[1] - sorted_x[0])
    725    xstop = sorted_x[-1] + 0.5 * (sorted_x[-1] - sorted_x[-2])
    726    x_samples = np.arange(xstart, xstop + 0.01, 0.01)
    727
    728    plt.figure(figsize=(8, 8 / 1.618))
    729    gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
    730    ax0 = plt.subplot(gs[0])
    731    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')
    732    ax0.plot(x_samples, func([o.value for o in fit_res], x_samples), label='Fit', zorder=10, ls='-', ms=0)
    733    ax0.set_xticklabels([])
    734    ax0.set_xlim([xstart, xstop])
    735    ax0.set_xticklabels([])
    736    ax0.legend(title=title)
    737
    738    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])
    739    ax1 = plt.subplot(gs[1])
    740    ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
    741    ax1.tick_params(direction='out')
    742    ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
    743    ax1.axhline(y=0.0, ls='--', color='k', marker=" ")
    744    ax1.fill_between(x_samples, -1.0, 1.0, alpha=0.1, facecolor='k')
    745    ax1.set_xlim([xstart, xstop])
    746    ax1.set_ylabel('Residuals')
    747    plt.subplots_adjust(wspace=None, hspace=None)
    748    plt.draw()
    749
    750
    751def error_band(x, func, beta):
    752    """Calculate the error band for an array of sample values x, for given fit function func with optimized parameters beta.
    753
    754    Returns
    755    -------
    756    err : np.array(Obs)
    757        Error band for an array of sample values x
    758    """
    759    cov = covariance(beta)
    760    if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps):
    761        warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning)
    762
    763    deriv = []
    764    for i, item in enumerate(x):
    765        deriv.append(np.array(egrad(func)([o.value for o in beta], item)))
    766
    767    err = []
    768    for i, item in enumerate(x):
    769        err.append(np.sqrt(deriv[i] @ cov @ deriv[i]))
    770    err = np.array(err)
    771
    772    return err
    773
    774
    775def ks_test(objects=None):
    776    """Performs a Kolmogorov–Smirnov test for the p-values of all fit object.
    777
    778    Parameters
    779    ----------
    780    objects : list
    781        List of fit results to include in the analysis (optional).
    782
    783    Returns
    784    -------
    785    None
    786    """
    787
    788    if objects is None:
    789        obs_list = []
    790        for obj in gc.get_objects():
    791            if isinstance(obj, Fit_result):
    792                obs_list.append(obj)
    793    else:
    794        obs_list = objects
    795
    796    p_values = [o.p_value for o in obs_list]
    797
    798    bins = len(p_values)
    799    x = np.arange(0, 1.001, 0.001)
    800    plt.plot(x, x, 'k', zorder=1)
    801    plt.xlim(0, 1)
    802    plt.ylim(0, 1)
    803    plt.xlabel('p-value')
    804    plt.ylabel('Cumulative probability')
    805    plt.title(str(bins) + ' p-values')
    806
    807    n = np.arange(1, bins + 1) / np.float64(bins)
    808    Xs = np.sort(p_values)
    809    plt.step(Xs, n)
    810    diffs = n - Xs
    811    loc_max_diff = np.argmax(np.abs(diffs))
    812    loc = Xs[loc_max_diff]
    813    plt.annotate('', xy=(loc, loc), xytext=(loc, loc + diffs[loc_max_diff]), arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
    814    plt.draw()
    815
    816    print(scipy.stats.kstest(p_values, 'uniform'))
    817
    818
    819def _extract_val_and_dval(string):
    820    split_string = string.split('(')
    821    if '.' in split_string[0] and '.' not in split_string[1][:-1]:
    822        factor = 10 ** -len(split_string[0].partition('.')[2])
    823    else:
    824        factor = 1
    825    return float(split_string[0]), float(split_string[1][:-1]) * factor

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 i in range(100):
    221            try:
    222                funcd[key](np.arange(i), 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 = i
    232        n_parms_ls.append(n_parms)
    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                if isinstance(i_prior, Obs):
    245                    loc_priors.append(i_prior)
    246                elif isinstance(i_prior, str):
    247                    loc_val, loc_dval = _extract_val_and_dval(i_prior)
    248                    loc_priors.append(cov_Obs(loc_val, loc_dval ** 2, '#prior' + str(i_n) + f"_{np.random.randint(2147483647):010d}"))
    249                else:
    250                    raise TypeError("Prior entries need to be 'Obs' or 'str'.")
    251
    252            prior_mask = np.arange(len(priors))
    253            output.priors = loc_priors
    254
    255        elif isinstance(priors, dict):
    256            loc_priors = []
    257            prior_mask = []
    258            output.priors = {}
    259            for pos, prior in priors.items():
    260                if isinstance(pos, int):
    261                    prior_mask.append(pos)
    262                else:
    263                    raise TypeError("Prior position needs to be an integer.")
    264                if isinstance(prior, Obs):
    265                    loc_priors.append(prior)
    266                elif isinstance(prior, str):
    267                    loc_val, loc_dval = _extract_val_and_dval(prior)
    268                    loc_priors.append(cov_Obs(loc_val, loc_dval ** 2, '#prior' + str(pos) + f"_{np.random.randint(2147483647):010d}"))
    269                else:
    270                    raise TypeError("Prior entries need to be 'Obs' or 'str'.")
    271                output.priors[pos] = loc_priors[-1]
    272            if max(prior_mask) >= n_parms:
    273                raise ValueError("Prior position out of range.")
    274        else:
    275            raise TypeError("Unkown type for `priors`.")
    276
    277        p_f = [o.value for o in loc_priors]
    278        dp_f = [o.dvalue for o in loc_priors]
    279        if np.any(np.asarray(dp_f) <= 0.0):
    280            raise Exception("No prior errors available, run the gamma method first.")
    281    else:
    282        p_f = dp_f = np.array([])
    283        prior_mask = []
    284        loc_priors = []
    285
    286    if 'initial_guess' in kwargs:
    287        x0 = kwargs.get('initial_guess')
    288        if len(x0) != n_parms:
    289            raise ValueError('Initial guess does not have the correct length: %d vs. %d' % (len(x0), n_parms))
    290    else:
    291        x0 = [0.1] * n_parms
    292
    293    if priors is None:
    294        def general_chisqfunc_uncorr(p, ivars, pr):
    295            model = anp.concatenate([anp.array(funcd[key](p, xd[key])).reshape(-1) for key in key_ls])
    296            return (ivars - model) / dy_f
    297    else:
    298        def general_chisqfunc_uncorr(p, ivars, pr):
    299            model = anp.concatenate([anp.array(funcd[key](p, xd[key])).reshape(-1) for key in key_ls])
    300            return anp.concatenate(((ivars - model) / dy_f, (p[prior_mask] - pr) / dp_f))
    301
    302    def chisqfunc_uncorr(p):
    303        return anp.sum(general_chisqfunc_uncorr(p, y_f, p_f) ** 2)
    304
    305    if kwargs.get('correlated_fit') is True:
    306        corr = covariance(y_all, correlation=True, **kwargs)
    307        covdiag = np.diag(1 / np.asarray(dy_f))
    308        condn = np.linalg.cond(corr)
    309        if condn > 0.1 / np.finfo(float).eps:
    310            raise Exception(f"Cannot invert correlation matrix as its condition number exceeds machine precision ({condn:1.2e})")
    311        if condn > 1e13:
    312            warnings.warn("Correlation matrix may be ill-conditioned, condition number: {%1.2e}" % (condn), RuntimeWarning)
    313        chol = np.linalg.cholesky(corr)
    314        chol_inv = scipy.linalg.solve_triangular(chol, covdiag, lower=True)
    315
    316        def general_chisqfunc(p, ivars, pr):
    317            model = anp.concatenate([anp.array(funcd[key](p, xd[key])).reshape(-1) for key in key_ls])
    318            return anp.concatenate((anp.dot(chol_inv, (ivars - model)), (p[prior_mask] - pr) / dp_f))
    319
    320        def chisqfunc(p):
    321            return anp.sum(general_chisqfunc(p, y_f, p_f) ** 2)
    322    else:
    323        general_chisqfunc = general_chisqfunc_uncorr
    324        chisqfunc = chisqfunc_uncorr
    325
    326    output.method = kwargs.get('method', 'Levenberg-Marquardt')
    327    if not silent:
    328        print('Method:', output.method)
    329
    330    if output.method != 'Levenberg-Marquardt':
    331        if output.method == 'migrad':
    332            tolerance = 1e-4  # default value of 1e-1 set by iminuit can be problematic
    333            if 'tol' in kwargs:
    334                tolerance = kwargs.get('tol')
    335            fit_result = iminuit.minimize(chisqfunc_uncorr, x0, tol=tolerance)  # Stopping criterion 0.002 * tol * errordef
    336            if kwargs.get('correlated_fit') is True:
    337                fit_result = iminuit.minimize(chisqfunc, fit_result.x, tol=tolerance)
    338            output.iterations = fit_result.nfev
    339        else:
    340            tolerance = 1e-12
    341            if 'tol' in kwargs:
    342                tolerance = kwargs.get('tol')
    343            fit_result = scipy.optimize.minimize(chisqfunc_uncorr, x0, method=kwargs.get('method'), tol=tolerance)
    344            if kwargs.get('correlated_fit') is True:
    345                fit_result = scipy.optimize.minimize(chisqfunc, fit_result.x, method=kwargs.get('method'), tol=tolerance)
    346            output.iterations = fit_result.nit
    347
    348        chisquare = fit_result.fun
    349
    350    else:
    351        if 'tol' in kwargs:
    352            print('tol cannot be set for Levenberg-Marquardt')
    353
    354        def chisqfunc_residuals_uncorr(p):
    355            return general_chisqfunc_uncorr(p, y_f, p_f)
    356
    357        fit_result = scipy.optimize.least_squares(chisqfunc_residuals_uncorr, x0, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15)
    358        if kwargs.get('correlated_fit') is True:
    359
    360            def chisqfunc_residuals(p):
    361                return general_chisqfunc(p, y_f, p_f)
    362
    363            fit_result = scipy.optimize.least_squares(chisqfunc_residuals, fit_result.x, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15)
    364
    365        chisquare = np.sum(fit_result.fun ** 2)
    366        assert np.isclose(chisquare, chisqfunc(fit_result.x), atol=1e-14)
    367
    368        output.iterations = fit_result.nfev
    369
    370    if not fit_result.success:
    371        raise Exception('The minimization procedure did not converge.')
    372
    373    if x_all.shape[-1] - n_parms > 0:
    374        output.chisquare = chisquare
    375        output.dof = x_all.shape[-1] - n_parms + len(loc_priors)
    376        output.chisquare_by_dof = output.chisquare / output.dof
    377        output.p_value = 1 - scipy.stats.chi2.cdf(output.chisquare, output.dof)
    378    else:
    379        output.chisquare_by_dof = float('nan')
    380
    381    output.message = fit_result.message
    382    if not silent:
    383        print(fit_result.message)
    384        print('chisquare/d.o.f.:', output.chisquare_by_dof)
    385        print('fit parameters', fit_result.x)
    386
    387    def prepare_hat_matrix():
    388        hat_vector = []
    389        for key in key_ls:
    390            if (len(xd[key]) != 0):
    391                hat_vector.append(jacobian(funcd[key])(fit_result.x, xd[key]))
    392        hat_vector = [item for sublist in hat_vector for item in sublist]
    393        return hat_vector
    394
    395    if kwargs.get('expected_chisquare') is True:
    396        if kwargs.get('correlated_fit') is not True:
    397            W = np.diag(1 / np.asarray(dy_f))
    398            cov = covariance(y_all)
    399            hat_vector = prepare_hat_matrix()
    400            A = W @ hat_vector
    401            P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T
    402            expected_chisquare = np.trace((np.identity(x_all.shape[-1]) - P_phi) @ W @ cov @ W)
    403            output.chisquare_by_expected_chisquare = output.chisquare / expected_chisquare
    404            if not silent:
    405                print('chisquare/expected_chisquare:', output.chisquare_by_expected_chisquare)
    406
    407    fitp = fit_result.x
    408
    409    try:
    410        hess = hessian(chisqfunc)(fitp)
    411    except TypeError:
    412        raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None
    413
    414    len_y = len(y_f)
    415
    416    def chisqfunc_compact(d):
    417        return anp.sum(general_chisqfunc(d[:n_parms], d[n_parms: n_parms + len_y], d[n_parms + len_y:]) ** 2)
    418
    419    jac_jac_y = hessian(chisqfunc_compact)(np.concatenate((fitp, y_f, p_f)))
    420
    421    # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv
    422    try:
    423        deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms, n_parms:])
    424    except np.linalg.LinAlgError:
    425        raise Exception("Cannot invert hessian matrix.")
    426
    427    result = []
    428    for i in range(n_parms):
    429        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])))
    430
    431    output.fit_parameters = result
    432
    433    # Hotelling t-squared p-value for correlated fits.
    434    if kwargs.get('correlated_fit') is True:
    435        n_cov = np.min(np.vectorize(lambda x_all: x_all.N)(y_all))
    436        output.t2_p_value = 1 - scipy.stats.f.cdf((n_cov - output.dof) / (output.dof * (n_cov - 1)) * output.chisquare,
    437                                                  output.dof, n_cov - output.dof)
    438
    439    if kwargs.get('resplot') is True:
    440        for key in key_ls:
    441            residual_plot(xd[key], yd[key], funcd[key], result, title=key)
    442
    443    if kwargs.get('qqplot') is True:
    444        for key in key_ls:
    445            qqplot(xd[key], yd[key], funcd[key], result, title=key)
    446
    447    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

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

    654def fit_lin(x, y, **kwargs):
    655    """Performs a linear fit to y = n + m * x and returns two Obs n, m.
    656
    657    Parameters
    658    ----------
    659    x : list
    660        Can either be a list of floats in which case no xerror is assumed, or
    661        a list of Obs, where the dvalues of the Obs are used as xerror for the fit.
    662    y : list
    663        List of Obs, the dvalues of the Obs are used as yerror for the fit.
    664
    665    Returns
    666    -------
    667    fit_parameters : list[Obs]
    668        LIist of fitted observables.
    669    """
    670
    671    def f(a, x):
    672        y = a[0] + a[1] * x
    673        return y
    674
    675    if all(isinstance(n, Obs) for n in x):
    676        out = total_least_squares(x, y, f, **kwargs)
    677        return out.fit_parameters
    678    elif all(isinstance(n, float) or isinstance(n, int) for n in x) or isinstance(x, np.ndarray):
    679        out = least_squares(x, y, f, **kwargs)
    680        return out.fit_parameters
    681    else:
    682        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

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

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

    752def error_band(x, func, beta):
    753    """Calculate the error band for an array of sample values x, for given fit function func with optimized parameters beta.
    754
    755    Returns
    756    -------
    757    err : np.array(Obs)
    758        Error band for an array of sample values x
    759    """
    760    cov = covariance(beta)
    761    if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps):
    762        warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning)
    763
    764    deriv = []
    765    for i, item in enumerate(x):
    766        deriv.append(np.array(egrad(func)([o.value for o in beta], item)))
    767
    768    err = []
    769    for i, item in enumerate(x):
    770        err.append(np.sqrt(deriv[i] @ cov @ deriv[i]))
    771    err = np.array(err)
    772
    773    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

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