pyerrors.fits

  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):
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):
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):
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):
 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):
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):
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=''):
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=''):
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):
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):
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