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