pyerrors.fits

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