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

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

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

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

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

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

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

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