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

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    inv_chol_cov_matrix [array,list], optional
    156        array: shape = (no of y values) X (no of y values)
    157        list:   for an uncombined fit: [""]
    158                for a combined fit: list of keys belonging to the corr_matrix saved in the array, must be the same as the keys of the y dict in alphabetical order
    159        If correlated_fit=True is set as well, can provide an inverse covariance matrix (y errors, dy_f included!) of your own choosing for a correlated fit.
    160        The matrix must be a lower triangular matrix constructed from a Cholesky decomposition: The function invert_corr_cov_cholesky(corr, inverrdiag) can be
    161        used to construct it from a correlation matrix (corr) and the errors dy_f of the data points (inverrdiag = np.diag(1 / np.asarray(dy_f))). For the correct
    162        ordering the correlation matrix (corr) can be sorted via the function sort_corr(corr, kl, yd) where kl is the list of keys and yd the y dict.
    163    expected_chisquare : bool
    164        If True estimates the expected chisquare which is
    165        corrected by effects caused by correlated input data (default False).
    166    resplot : bool
    167        If True, a plot which displays fit, data and residuals is generated (default False).
    168    qqplot : bool
    169        If True, a quantile-quantile plot of the fit result is generated (default False).
    170    num_grad : bool
    171        Use numerical differentation instead of automatic differentiation to perform the error propagation (default False).
    172
    173    Returns
    174    -------
    175    output : Fit_result
    176        Parameters and information on the fitted result.
    177    Examples
    178    ------
    179    >>> # Example of a correlated (correlated_fit = True, inv_chol_cov_matrix handed over) combined fit, based on a randomly generated data set
    180    >>> import numpy as np
    181    >>> from scipy.stats import norm
    182    >>> from scipy.linalg import cholesky
    183    >>> import pyerrors as pe
    184    >>> # generating the random data set
    185    >>> num_samples = 400
    186    >>> N = 3
    187    >>> x = np.arange(N)
    188    >>> x1 = norm.rvs(size=(N, num_samples)) # generate random numbers
    189    >>> x2 = norm.rvs(size=(N, num_samples)) # generate random numbers
    190    >>> r = r1 = r2 = np.zeros((N, N))
    191    >>> y = {}
    192    >>> for i in range(N):
    193    >>>    for j in range(N):
    194    >>>        r[i, j] = np.exp(-0.8 * np.fabs(i - j)) # element in correlation matrix
    195    >>> errl = np.sqrt([3.4, 2.5, 3.6]) # set y errors
    196    >>> for i in range(N):
    197    >>>    for j in range(N):
    198    >>>        r[i, j] *= errl[i] * errl[j] # element in covariance matrix
    199    >>> c = cholesky(r, lower=True)
    200    >>> y = {'a': np.dot(c, x1), 'b': np.dot(c, x2)} # generate y data with the covariance matrix defined
    201    >>> # random data set has been generated, now the dictionaries and the inverse covariance matrix to be handed over are built
    202    >>> x_dict = {}
    203    >>> y_dict = {}
    204    >>> chol_inv_dict = {}
    205    >>> data = []
    206    >>> for key in y.keys():
    207    >>>    x_dict[key] = x
    208    >>>    for i in range(N):
    209    >>>        data.append(pe.Obs([[i + 1 + o for o in y[key][i]]], ['ens'])) # generate y Obs from the y data
    210    >>>    [o.gamma_method() for o in data]
    211    >>>    corr = pe.covariance(data, correlation=True)
    212    >>>    inverrdiag = np.diag(1 / np.asarray([o.dvalue for o in data]))
    213    >>>    chol_inv = pe.obs.invert_corr_cov_cholesky(corr, inverrdiag) # gives form of the inverse covariance matrix needed for the combined correlated fit below
    214    >>> y_dict = {'a': data[:3], 'b': data[3:]}
    215    >>> # common fit parameter p[0] in combined fit
    216    >>> def fit1(p, x):
    217    >>>    return p[0] + p[1] * x
    218    >>> def fit2(p, x):
    219    >>>    return p[0] + p[2] * x
    220    >>> fitf_dict = {'a': fit1, 'b':fit2}
    221    >>> fitp_inv_cov_combined_fit = pe.least_squares(x_dict,y_dict, fitf_dict, correlated_fit = True, inv_chol_cov_matrix = [chol_inv,['a','b']])
    222    Fit with 3 parameters
    223    Method: Levenberg-Marquardt
    224    `ftol` termination condition is satisfied.
    225    chisquare/d.o.f.: 0.5388013574561786 # random
    226    fit parameters [1.11897846 0.96361162 0.92325319] # random
    227
    228    '''
    229    output = Fit_result()
    230
    231    if (isinstance(x, dict) and isinstance(y, dict) and isinstance(func, dict)):
    232        xd = {key: anp.asarray(x[key]) for key in x}
    233        yd = y
    234        funcd = func
    235        output.fit_function = func
    236    elif (isinstance(x, dict) or isinstance(y, dict) or isinstance(func, dict)):
    237        raise TypeError("All arguments have to be dictionaries in order to perform a combined fit.")
    238    else:
    239        x = np.asarray(x)
    240        xd = {"": x}
    241        yd = {"": y}
    242        funcd = {"": func}
    243        output.fit_function = func
    244
    245    if kwargs.get('num_grad') is True:
    246        jacobian = num_jacobian
    247        hessian = num_hessian
    248    else:
    249        jacobian = auto_jacobian
    250        hessian = auto_hessian
    251
    252    key_ls = sorted(list(xd.keys()))
    253
    254    if sorted(list(yd.keys())) != key_ls:
    255        raise ValueError('x and y dictionaries do not contain the same keys.')
    256
    257    if sorted(list(funcd.keys())) != key_ls:
    258        raise ValueError('x and func dictionaries do not contain the same keys.')
    259
    260    x_all = np.concatenate([np.array(xd[key]).transpose() for key in key_ls]).transpose()
    261    y_all = np.concatenate([np.array(yd[key]) for key in key_ls])
    262
    263    y_f = [o.value for o in y_all]
    264    dy_f = [o.dvalue for o in y_all]
    265
    266    if len(x_all.shape) > 2:
    267        raise ValueError("Unknown format for x values")
    268
    269    if np.any(np.asarray(dy_f) <= 0.0):
    270        raise Exception("No y errors available, run the gamma method first.")
    271
    272    # number of fit parameters
    273    n_parms_ls = []
    274    for key in key_ls:
    275        if not callable(funcd[key]):
    276            raise TypeError('func (key=' + key + ') is not a function.')
    277        if np.asarray(xd[key]).shape[-1] != len(yd[key]):
    278            raise ValueError('x and y input (key=' + key + ') do not have the same length')
    279        for n_loc in range(100):
    280            try:
    281                funcd[key](np.arange(n_loc), x_all.T[0])
    282            except TypeError:
    283                continue
    284            except IndexError:
    285                continue
    286            else:
    287                break
    288        else:
    289            raise RuntimeError("Fit function (key=" + key + ") is not valid.")
    290        n_parms_ls.append(n_loc)
    291
    292    n_parms = max(n_parms_ls)
    293
    294    if len(key_ls) > 1:
    295        for key in key_ls:
    296            if np.asarray(yd[key]).shape != funcd[key](np.arange(n_parms), xd[key]).shape:
    297                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.")
    298
    299    if not silent:
    300        print('Fit with', n_parms, 'parameter' + 's' * (n_parms > 1))
    301
    302    if priors is not None:
    303        if isinstance(priors, (list, np.ndarray)):
    304            if n_parms != len(priors):
    305                raise ValueError("'priors' does not have the correct length.")
    306
    307            loc_priors = []
    308            for i_n, i_prior in enumerate(priors):
    309                loc_priors.append(_construct_prior_obs(i_prior, i_n))
    310
    311            prior_mask = np.arange(len(priors))
    312            output.priors = loc_priors
    313
    314        elif isinstance(priors, dict):
    315            loc_priors = []
    316            prior_mask = []
    317            output.priors = {}
    318            for pos, prior in priors.items():
    319                if isinstance(pos, int):
    320                    prior_mask.append(pos)
    321                else:
    322                    raise TypeError("Prior position needs to be an integer.")
    323                loc_priors.append(_construct_prior_obs(prior, pos))
    324
    325                output.priors[pos] = loc_priors[-1]
    326            if max(prior_mask) >= n_parms:
    327                raise ValueError("Prior position out of range.")
    328        else:
    329            raise TypeError("Unkown type for `priors`.")
    330
    331        p_f = [o.value for o in loc_priors]
    332        dp_f = [o.dvalue for o in loc_priors]
    333        if np.any(np.asarray(dp_f) <= 0.0):
    334            raise Exception("No prior errors available, run the gamma method first.")
    335    else:
    336        p_f = dp_f = np.array([])
    337        prior_mask = []
    338        loc_priors = []
    339
    340    if 'initial_guess' in kwargs:
    341        x0 = kwargs.get('initial_guess')
    342        if len(x0) != n_parms:
    343            raise ValueError('Initial guess does not have the correct length: %d vs. %d' % (len(x0), n_parms))
    344    else:
    345        x0 = [0.1] * n_parms
    346
    347    if priors is None:
    348        def general_chisqfunc_uncorr(p, ivars, pr):
    349            model = anp.concatenate([anp.array(funcd[key](p, xd[key])).reshape(-1) for key in key_ls])
    350            return (ivars - model) / dy_f
    351    else:
    352        def general_chisqfunc_uncorr(p, ivars, pr):
    353            model = anp.concatenate([anp.array(funcd[key](p, xd[key])).reshape(-1) for key in key_ls])
    354            return anp.concatenate(((ivars - model) / dy_f, (p[prior_mask] - pr) / dp_f))
    355
    356    def chisqfunc_uncorr(p):
    357        return anp.sum(general_chisqfunc_uncorr(p, y_f, p_f) ** 2)
    358
    359    if kwargs.get('correlated_fit') is True:
    360        if 'inv_chol_cov_matrix' in kwargs:
    361            chol_inv = kwargs.get('inv_chol_cov_matrix')
    362            if (chol_inv[0].shape[0] != len(dy_f)):
    363                raise TypeError('The number of columns of the inverse covariance matrix handed over needs to be equal to the number of y errors.')
    364            if (chol_inv[0].shape[0] != chol_inv[0].shape[1]):
    365                raise TypeError('The inverse covariance matrix handed over needs to have the same number of rows as columns.')
    366            if (chol_inv[1] != key_ls):
    367                raise ValueError('The keys of inverse covariance matrix are not the same or do not appear in the same order as the x and y values.')
    368            chol_inv = chol_inv[0]
    369        else:
    370            corr = covariance(y_all, correlation=True, **kwargs)
    371            inverrdiag = np.diag(1 / np.asarray(dy_f))
    372            chol_inv = invert_corr_cov_cholesky(corr, inverrdiag)
    373
    374        def general_chisqfunc(p, ivars, pr):
    375            model = anp.concatenate([anp.array(funcd[key](p, xd[key])).reshape(-1) for key in key_ls])
    376            return anp.concatenate((anp.dot(chol_inv, (ivars - model)), (p[prior_mask] - pr) / dp_f))
    377
    378        def chisqfunc(p):
    379            return anp.sum(general_chisqfunc(p, y_f, p_f) ** 2)
    380    else:
    381        general_chisqfunc = general_chisqfunc_uncorr
    382        chisqfunc = chisqfunc_uncorr
    383
    384    output.method = kwargs.get('method', 'Levenberg-Marquardt')
    385    if not silent:
    386        print('Method:', output.method)
    387
    388    if output.method != 'Levenberg-Marquardt':
    389        if output.method == 'migrad':
    390            tolerance = 1e-4  # default value of 1e-1 set by iminuit can be problematic
    391            if 'tol' in kwargs:
    392                tolerance = kwargs.get('tol')
    393            fit_result = iminuit.minimize(chisqfunc_uncorr, x0, tol=tolerance)  # Stopping criterion 0.002 * tol * errordef
    394            if kwargs.get('correlated_fit') is True:
    395                fit_result = iminuit.minimize(chisqfunc, fit_result.x, tol=tolerance)
    396            output.iterations = fit_result.nfev
    397        else:
    398            tolerance = 1e-12
    399            if 'tol' in kwargs:
    400                tolerance = kwargs.get('tol')
    401            fit_result = scipy.optimize.minimize(chisqfunc_uncorr, x0, method=kwargs.get('method'), tol=tolerance)
    402            if kwargs.get('correlated_fit') is True:
    403                fit_result = scipy.optimize.minimize(chisqfunc, fit_result.x, method=kwargs.get('method'), tol=tolerance)
    404            output.iterations = fit_result.nit
    405
    406        chisquare = fit_result.fun
    407
    408    else:
    409        if 'tol' in kwargs:
    410            print('tol cannot be set for Levenberg-Marquardt')
    411
    412        def chisqfunc_residuals_uncorr(p):
    413            return general_chisqfunc_uncorr(p, y_f, p_f)
    414
    415        fit_result = scipy.optimize.least_squares(chisqfunc_residuals_uncorr, x0, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15)
    416        if kwargs.get('correlated_fit') is True:
    417            def chisqfunc_residuals(p):
    418                return general_chisqfunc(p, y_f, p_f)
    419
    420            fit_result = scipy.optimize.least_squares(chisqfunc_residuals, fit_result.x, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15)
    421
    422        chisquare = np.sum(fit_result.fun ** 2)
    423        assert np.isclose(chisquare, chisqfunc(fit_result.x), atol=1e-14)
    424
    425        output.iterations = fit_result.nfev
    426
    427    if not fit_result.success:
    428        raise Exception('The minimization procedure did not converge.')
    429
    430    output.chisquare = chisquare
    431    output.dof = y_all.shape[-1] - n_parms + len(loc_priors)
    432    output.p_value = 1 - scipy.stats.chi2.cdf(output.chisquare, output.dof)
    433    if output.dof > 0:
    434        output.chisquare_by_dof = output.chisquare / output.dof
    435    else:
    436        output.chisquare_by_dof = float('nan')
    437
    438    output.message = fit_result.message
    439    if not silent:
    440        print(fit_result.message)
    441        print('chisquare/d.o.f.:', output.chisquare_by_dof)
    442        print('fit parameters', fit_result.x)
    443
    444    def prepare_hat_matrix():
    445        hat_vector = []
    446        for key in key_ls:
    447            if (len(xd[key]) != 0):
    448                hat_vector.append(jacobian(funcd[key])(fit_result.x, xd[key]))
    449        hat_vector = [item for sublist in hat_vector for item in sublist]
    450        return hat_vector
    451
    452    if kwargs.get('expected_chisquare') is True:
    453        if kwargs.get('correlated_fit') is not True:
    454            W = np.diag(1 / np.asarray(dy_f))
    455            cov = covariance(y_all)
    456            hat_vector = prepare_hat_matrix()
    457            A = W @ hat_vector
    458            P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T
    459            expected_chisquare = np.trace((np.identity(y_all.shape[-1]) - P_phi) @ W @ cov @ W)
    460            output.chisquare_by_expected_chisquare = output.chisquare / expected_chisquare
    461            if not silent:
    462                print('chisquare/expected_chisquare:', output.chisquare_by_expected_chisquare)
    463
    464    fitp = fit_result.x
    465
    466    try:
    467        hess = hessian(chisqfunc)(fitp)
    468    except TypeError:
    469        raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None
    470
    471    len_y = len(y_f)
    472
    473    def chisqfunc_compact(d):
    474        return anp.sum(general_chisqfunc(d[:n_parms], d[n_parms: n_parms + len_y], d[n_parms + len_y:]) ** 2)
    475
    476    jac_jac_y = hessian(chisqfunc_compact)(np.concatenate((fitp, y_f, p_f)))
    477
    478    # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv
    479    try:
    480        deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms, n_parms:])
    481    except np.linalg.LinAlgError:
    482        raise Exception("Cannot invert hessian matrix.")
    483
    484    result = []
    485    for i in range(n_parms):
    486        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])))
    487
    488    output.fit_parameters = result
    489
    490    # Hotelling t-squared p-value for correlated fits.
    491    if kwargs.get('correlated_fit') is True:
    492        n_cov = np.min(np.vectorize(lambda x_all: x_all.N)(y_all))
    493        output.t2_p_value = 1 - scipy.stats.f.cdf((n_cov - output.dof) / (output.dof * (n_cov - 1)) * output.chisquare,
    494                                                  output.dof, n_cov - output.dof)
    495
    496    if kwargs.get('resplot') is True:
    497        for key in key_ls:
    498            residual_plot(xd[key], yd[key], funcd[key], result, title=key)
    499
    500    if kwargs.get('qqplot') is True:
    501        for key in key_ls:
    502            qqplot(xd[key], yd[key], funcd[key], result, title=key)
    503
    504    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).
- inv_chol_cov_matrix [array,list], optional: array: shape = (no of y values) X (no of y values) list: for an uncombined fit: [""] for a combined fit: list of keys belonging to the corr_matrix saved in the array, must be the same as the keys of the y dict in alphabetical order If correlated_fit=True is set as well, can provide an inverse covariance matrix (y errors, dy_f included!) of your own choosing for a correlated fit. The matrix must be a lower triangular matrix constructed from a Cholesky decomposition: The function invert_corr_cov_cholesky(corr, inverrdiag) can be used to construct it from a correlation matrix (corr) and the errors dy_f of the data points (inverrdiag = np.diag(1 / np.asarray(dy_f))). For the correct ordering the correlation matrix (corr) can be sorted via the function sort_corr(corr, kl, yd) where kl is the list of keys and yd the y dict.
- 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.

###### Examples

    >>> # Example of a correlated (correlated_fit = True, inv_chol_cov_matrix handed over) combined fit, based on a randomly generated data set
    >>> import numpy as np
    >>> from scipy.stats import norm
    >>> from scipy.linalg import cholesky
    >>> import pyerrors as pe
    >>> # generating the random data set
    >>> num_samples = 400
    >>> N = 3
    >>> x = np.arange(N)
    >>> x1 = norm.rvs(size=(N, num_samples)) # generate random numbers
    >>> x2 = norm.rvs(size=(N, num_samples)) # generate random numbers
    >>> r = r1 = r2 = np.zeros((N, N))
    >>> y = {}
    >>> for i in range(N):
    >>>    for j in range(N):
    >>>        r[i, j] = np.exp(-0.8 * np.fabs(i - j)) # element in correlation matrix
    >>> errl = np.sqrt([3.4, 2.5, 3.6]) # set y errors
    >>> for i in range(N):
    >>>    for j in range(N):
    >>>        r[i, j] *= errl[i] * errl[j] # element in covariance matrix
    >>> c = cholesky(r, lower=True)
    >>> y = {'a': np.dot(c, x1), 'b': np.dot(c, x2)} # generate y data with the covariance matrix defined
    >>> # random data set has been generated, now the dictionaries and the inverse covariance matrix to be handed over are built
    >>> x_dict = {}
    >>> y_dict = {}
    >>> chol_inv_dict = {}
    >>> data = []
    >>> for key in y.keys():
    >>>    x_dict[key] = x
    >>>    for i in range(N):
    >>>        data.append(pe.Obs([[i + 1 + o for o in y[key][i]]], ['ens'])) # generate y Obs from the y data
    >>>    [o.gamma_method() for o in data]
    >>>    corr = pe.covariance(data, correlation=True)
    >>>    inverrdiag = np.diag(1 / np.asarray([o.dvalue for o in data]))
    >>>    chol_inv = pe.obs.invert_corr_cov_cholesky(corr, inverrdiag) # gives form of the inverse covariance matrix needed for the combined correlated fit below
    >>> y_dict = {'a': data[:3], 'b': data[3:]}
    >>> # common fit parameter p[0] in combined fit
    >>> def fit1(p, x):
    >>>    return p[0] + p[1] * x
    >>> def fit2(p, x):
    >>>    return p[0] + p[2] * x
    >>> fitf_dict = {'a': fit1, 'b':fit2}
    >>> fitp_inv_cov_combined_fit = pe.least_squares(x_dict,y_dict, fitf_dict, correlated_fit = True, inv_chol_cov_matrix = [chol_inv,['a','b']])
    Fit with 3 parameters
    Method: Levenberg-Marquardt
    `ftol` termination condition is satisfied.
    chisquare/d.o.f.: 0.5388013574561786 # random
    fit parameters [1.11897846 0.96361162 0.92325319] # random

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

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

    711def fit_lin(x, y, **kwargs):
    712    """Performs a linear fit to y = n + m * x and returns two Obs n, m.
    713
    714    Parameters
    715    ----------
    716    x : list
    717        Can either be a list of floats in which case no xerror is assumed, or
    718        a list of Obs, where the dvalues of the Obs are used as xerror for the fit.
    719    y : list
    720        List of Obs, the dvalues of the Obs are used as yerror for the fit.
    721
    722    Returns
    723    -------
    724    fit_parameters : list[Obs]
    725        LIist of fitted observables.
    726    """
    727
    728    def f(a, x):
    729        y = a[0] + a[1] * x
    730        return y
    731
    732    if all(isinstance(n, Obs) for n in x):
    733        out = total_least_squares(x, y, f, **kwargs)
    734        return out.fit_parameters
    735    elif all(isinstance(n, float) or isinstance(n, int) for n in x) or isinstance(x, np.ndarray):
    736        out = least_squares(x, y, f, **kwargs)
    737        return out.fit_parameters
    738    else:
    739        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

    742def qqplot(x, o_y, func, p, title=""):
    743    """Generates a quantile-quantile plot of the fit result which can be used to
    744       check if the residuals of the fit are gaussian distributed.
    745
    746    Returns
    747    -------
    748    None
    749    """
    750
    751    residuals = []
    752    for i_x, i_y in zip(x, o_y):
    753        residuals.append((i_y - func(p, i_x)) / i_y.dvalue)
    754    residuals = sorted(residuals)
    755    my_y = [o.value for o in residuals]
    756    probplot = scipy.stats.probplot(my_y)
    757    my_x = probplot[0][0]
    758    plt.figure(figsize=(8, 8 / 1.618))
    759    plt.errorbar(my_x, my_y, fmt='o')
    760    fit_start = my_x[0]
    761    fit_stop = my_x[-1]
    762    samples = np.arange(fit_start, fit_stop, 0.01)
    763    plt.plot(samples, samples, 'k--', zorder=11, label='Standard normal distribution')
    764    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='-')
    765
    766    plt.xlabel('Theoretical quantiles')
    767    plt.ylabel('Ordered Values')
    768    plt.legend(title=title)
    769    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

    772def residual_plot(x, y, func, fit_res, title=""):
    773    """Generates a plot which compares the fit to the data and displays the corresponding residuals
    774
    775    For uncorrelated data the residuals are expected to be distributed ~N(0,1).
    776
    777    Returns
    778    -------
    779    None
    780    """
    781    sorted_x = sorted(x)
    782    xstart = sorted_x[0] - 0.5 * (sorted_x[1] - sorted_x[0])
    783    xstop = sorted_x[-1] + 0.5 * (sorted_x[-1] - sorted_x[-2])
    784    x_samples = np.arange(xstart, xstop + 0.01, 0.01)
    785
    786    plt.figure(figsize=(8, 8 / 1.618))
    787    gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
    788    ax0 = plt.subplot(gs[0])
    789    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')
    790    ax0.plot(x_samples, func([o.value for o in fit_res], x_samples), label='Fit', zorder=10, ls='-', ms=0)
    791    ax0.set_xticklabels([])
    792    ax0.set_xlim([xstart, xstop])
    793    ax0.set_xticklabels([])
    794    ax0.legend(title=title)
    795
    796    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])
    797    ax1 = plt.subplot(gs[1])
    798    ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
    799    ax1.tick_params(direction='out')
    800    ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
    801    ax1.axhline(y=0.0, ls='--', color='k', marker=" ")
    802    ax1.fill_between(x_samples, -1.0, 1.0, alpha=0.1, facecolor='k')
    803    ax1.set_xlim([xstart, xstop])
    804    ax1.set_ylabel('Residuals')
    805    plt.subplots_adjust(wspace=None, hspace=None)
    806    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

    809def error_band(x, func, beta):
    810    """Calculate the error band for an array of sample values x, for given fit function func with optimized parameters beta.
    811
    812    Returns
    813    -------
    814    err : np.array(Obs)
    815        Error band for an array of sample values x
    816    """
    817    cov = covariance(beta)
    818    if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps):
    819        warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning)
    820
    821    deriv = []
    822    for i, item in enumerate(x):
    823        deriv.append(np.array(egrad(func)([o.value for o in beta], item)))
    824
    825    err = []
    826    for i, item in enumerate(x):
    827        err.append(np.sqrt(deriv[i] @ cov @ deriv[i]))
    828    err = np.array(err)
    829
    830    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

    833def ks_test(objects=None):
    834    """Performs a Kolmogorov–Smirnov test for the p-values of all fit object.
    835
    836    Parameters
    837    ----------
    838    objects : list
    839        List of fit results to include in the analysis (optional).
    840
    841    Returns
    842    -------
    843    None
    844    """
    845
    846    if objects is None:
    847        obs_list = []
    848        for obj in gc.get_objects():
    849            if isinstance(obj, Fit_result):
    850                obs_list.append(obj)
    851    else:
    852        obs_list = objects
    853
    854    p_values = [o.p_value for o in obs_list]
    855
    856    bins = len(p_values)
    857    x = np.arange(0, 1.001, 0.001)
    858    plt.plot(x, x, 'k', zorder=1)
    859    plt.xlim(0, 1)
    860    plt.ylim(0, 1)
    861    plt.xlabel('p-value')
    862    plt.ylabel('Cumulative probability')
    863    plt.title(str(bins) + ' p-values')
    864
    865    n = np.arange(1, bins + 1) / np.float64(bins)
    866    Xs = np.sort(p_values)
    867    plt.step(Xs, n)
    868    diffs = n - Xs
    869    loc_max_diff = np.argmax(np.abs(diffs))
    870    loc = Xs[loc_max_diff]
    871    plt.annotate('', xy=(loc, loc), xytext=(loc, loc + diffs[loc_max_diff]), arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
    872    plt.draw()
    873
    874    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