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