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