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 Exception:
184            continue
185        else:
186            break
187    else:
188        raise RuntimeError("Fit function is not valid.")
189
190    n_parms = i
191    if not silent:
192        print('Fit with', n_parms, 'parameter' + 's' * (n_parms > 1))
193
194    x_f = np.vectorize(lambda o: o.value)(x)
195    dx_f = np.vectorize(lambda o: o.dvalue)(x)
196    y_f = np.array([o.value for o in y])
197    dy_f = np.array([o.dvalue for o in y])
198
199    if np.any(np.asarray(dx_f) <= 0.0):
200        raise Exception('No x errors available, run the gamma method first.')
201
202    if np.any(np.asarray(dy_f) <= 0.0):
203        raise Exception('No y errors available, run the gamma method first.')
204
205    if 'initial_guess' in kwargs:
206        x0 = kwargs.get('initial_guess')
207        if len(x0) != n_parms:
208            raise Exception('Initial guess does not have the correct length: %d vs. %d' % (len(x0), n_parms))
209    else:
210        x0 = [1] * n_parms
211
212    data = RealData(x_f, y_f, sx=dx_f, sy=dy_f)
213    model = Model(func)
214    odr = ODR(data, model, x0, partol=np.finfo(np.float64).eps)
215    odr.set_job(fit_type=0, deriv=1)
216    out = odr.run()
217
218    output.residual_variance = out.res_var
219
220    output.method = 'ODR'
221
222    output.message = out.stopreason
223
224    output.xplus = out.xplus
225
226    if not silent:
227        print('Method: ODR')
228        print(*out.stopreason)
229        print('Residual variance:', output.residual_variance)
230
231    if out.info > 3:
232        raise Exception('The minimization procedure did not converge.')
233
234    m = x_f.size
235
236    def odr_chisquare(p):
237        model = func(p[:n_parms], p[n_parms:].reshape(x_shape))
238        chisq = anp.sum(((y_f - model) / dy_f) ** 2) + anp.sum(((x_f - p[n_parms:].reshape(x_shape)) / dx_f) ** 2)
239        return chisq
240
241    if kwargs.get('expected_chisquare') is True:
242        W = np.diag(1 / np.asarray(np.concatenate((dy_f.ravel(), dx_f.ravel()))))
243
244        if kwargs.get('covariance') is not None:
245            cov = kwargs.get('covariance')
246        else:
247            cov = covariance(np.concatenate((y, x.ravel())))
248
249        number_of_x_parameters = int(m / x_f.shape[-1])
250
251        old_jac = jacobian(func)(out.beta, out.xplus)
252        fused_row1 = np.concatenate((old_jac, np.concatenate((number_of_x_parameters * [np.zeros(old_jac.shape)]), axis=0)))
253        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])))
254        new_jac = np.concatenate((fused_row1, fused_row2), axis=1)
255
256        A = W @ new_jac
257        P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T
258        expected_chisquare = np.trace((np.identity(P_phi.shape[0]) - P_phi) @ W @ cov @ W)
259        if expected_chisquare <= 0.0:
260            warnings.warn("Negative expected_chisquare.", RuntimeWarning)
261            expected_chisquare = np.abs(expected_chisquare)
262        output.chisquare_by_expected_chisquare = odr_chisquare(np.concatenate((out.beta, out.xplus.ravel()))) / expected_chisquare
263        if not silent:
264            print('chisquare/expected_chisquare:',
265                  output.chisquare_by_expected_chisquare)
266
267    fitp = out.beta
268    try:
269        hess = jacobian(jacobian(odr_chisquare))(np.concatenate((fitp, out.xplus.ravel())))
270    except TypeError:
271        raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None
272
273    def odr_chisquare_compact_x(d):
274        model = func(d[:n_parms], d[n_parms:n_parms + m].reshape(x_shape))
275        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)
276        return chisq
277
278    jac_jac_x = jacobian(jacobian(odr_chisquare_compact_x))(np.concatenate((fitp, out.xplus.ravel(), x_f.ravel())))
279
280    # Compute hess^{-1} @ jac_jac_x[:n_parms + m, n_parms + m:] using LAPACK dgesv
281    try:
282        deriv_x = -scipy.linalg.solve(hess, jac_jac_x[:n_parms + m, n_parms + m:])
283    except np.linalg.LinAlgError:
284        raise Exception("Cannot invert hessian matrix.")
285
286    def odr_chisquare_compact_y(d):
287        model = func(d[:n_parms], d[n_parms:n_parms + m].reshape(x_shape))
288        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)
289        return chisq
290
291    jac_jac_y = jacobian(jacobian(odr_chisquare_compact_y))(np.concatenate((fitp, out.xplus.ravel(), y_f)))
292
293    # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv
294    try:
295        deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms + m, n_parms + m:])
296    except np.linalg.LinAlgError:
297        raise Exception("Cannot invert hessian matrix.")
298
299    result = []
300    for i in range(n_parms):
301        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])))
302
303    output.fit_parameters = result
304
305    output.odr_chisquare = odr_chisquare(np.concatenate((out.beta, out.xplus.ravel())))
306    output.dof = x.shape[-1] - n_parms
307    output.p_value = 1 - scipy.stats.chi2.cdf(output.odr_chisquare, output.dof)
308
309    return output
310
311
312def _prior_fit(x, y, func, priors, silent=False, **kwargs):
313    output = Fit_result()
314
315    output.fit_function = func
316
317    x = np.asarray(x)
318
319    if not callable(func):
320        raise TypeError('func has to be a function.')
321
322    for i in range(100):
323        try:
324            func(np.arange(i), 0)
325        except Exception:
326            continue
327        else:
328            break
329    else:
330        raise RuntimeError("Fit function is not valid.")
331
332    n_parms = i
333
334    if n_parms != len(priors):
335        raise Exception('Priors does not have the correct length.')
336
337    def extract_val_and_dval(string):
338        split_string = string.split('(')
339        if '.' in split_string[0] and '.' not in split_string[1][:-1]:
340            factor = 10 ** -len(split_string[0].partition('.')[2])
341        else:
342            factor = 1
343        return float(split_string[0]), float(split_string[1][:-1]) * factor
344
345    loc_priors = []
346    for i_n, i_prior in enumerate(priors):
347        if isinstance(i_prior, Obs):
348            loc_priors.append(i_prior)
349        else:
350            loc_val, loc_dval = extract_val_and_dval(i_prior)
351            loc_priors.append(cov_Obs(loc_val, loc_dval ** 2, '#prior' + str(i_n) + f"_{np.random.randint(2147483647):010d}"))
352
353    output.priors = loc_priors
354
355    if not silent:
356        print('Fit with', n_parms, 'parameter' + 's' * (n_parms > 1))
357
358    y_f = [o.value for o in y]
359    dy_f = [o.dvalue for o in y]
360
361    if np.any(np.asarray(dy_f) <= 0.0):
362        raise Exception('No y errors available, run the gamma method first.')
363
364    p_f = [o.value for o in loc_priors]
365    dp_f = [o.dvalue for o in loc_priors]
366
367    if np.any(np.asarray(dp_f) <= 0.0):
368        raise Exception('No prior errors available, run the gamma method first.')
369
370    if 'initial_guess' in kwargs:
371        x0 = kwargs.get('initial_guess')
372        if len(x0) != n_parms:
373            raise Exception('Initial guess does not have the correct length.')
374    else:
375        x0 = p_f
376
377    def chisqfunc(p):
378        model = func(p, x)
379        chisq = anp.sum(((y_f - model) / dy_f) ** 2) + anp.sum(((p_f - p) / dp_f) ** 2)
380        return chisq
381
382    if not silent:
383        print('Method: migrad')
384
385    m = iminuit.Minuit(chisqfunc, x0)
386    m.errordef = 1
387    m.print_level = 0
388    if 'tol' in kwargs:
389        m.tol = kwargs.get('tol')
390    else:
391        m.tol = 1e-4
392    m.migrad()
393    params = np.asarray(m.values)
394
395    output.chisquare_by_dof = m.fval / len(x)
396
397    output.method = 'migrad'
398
399    if not silent:
400        print('chisquare/d.o.f.:', output.chisquare_by_dof)
401
402    if not m.fmin.is_valid:
403        raise Exception('The minimization procedure did not converge.')
404
405    hess_inv = np.linalg.pinv(jacobian(jacobian(chisqfunc))(params))
406
407    def chisqfunc_compact(d):
408        model = func(d[:n_parms], x)
409        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)
410        return chisq
411
412    jac_jac = jacobian(jacobian(chisqfunc_compact))(np.concatenate((params, y_f, p_f)))
413
414    deriv = -hess_inv @ jac_jac[:n_parms, n_parms:]
415
416    result = []
417    for i in range(n_parms):
418        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])))
419
420    output.fit_parameters = result
421    output.chisquare = chisqfunc(np.asarray(params))
422
423    if kwargs.get('resplot') is True:
424        residual_plot(x, y, func, result)
425
426    if kwargs.get('qqplot') is True:
427        qqplot(x, y, func, result)
428
429    return output
430
431
432def _standard_fit(x, y, func, silent=False, **kwargs):
433
434    output = Fit_result()
435
436    output.fit_function = func
437
438    x = np.asarray(x)
439
440    if x.shape[-1] != len(y):
441        raise Exception('x and y input have to have the same length')
442
443    if len(x.shape) > 2:
444        raise Exception('Unknown format for x values')
445
446    if not callable(func):
447        raise TypeError('func has to be a function.')
448
449    for i in range(42):
450        try:
451            func(np.arange(i), x.T[0])
452        except Exception:
453            continue
454        else:
455            break
456    else:
457        raise RuntimeError("Fit function is not valid.")
458
459    n_parms = i
460
461    if not silent:
462        print('Fit with', n_parms, 'parameter' + 's' * (n_parms > 1))
463
464    y_f = [o.value for o in y]
465    dy_f = [o.dvalue for o in y]
466
467    if np.any(np.asarray(dy_f) <= 0.0):
468        raise Exception('No y errors available, run the gamma method first.')
469
470    if 'initial_guess' in kwargs:
471        x0 = kwargs.get('initial_guess')
472        if len(x0) != n_parms:
473            raise Exception('Initial guess does not have the correct length: %d vs. %d' % (len(x0), n_parms))
474    else:
475        x0 = [0.1] * n_parms
476
477    if kwargs.get('correlated_fit') is True:
478        corr = covariance(y, correlation=True, **kwargs)
479        covdiag = np.diag(1 / np.asarray(dy_f))
480        condn = np.linalg.cond(corr)
481        if condn > 0.1 / np.finfo(float).eps:
482            raise Exception(f"Cannot invert correlation matrix as its condition number exceeds machine precision ({condn:1.2e})")
483        if condn > 1 / np.sqrt(np.finfo(float).eps):
484            warnings.warn("Correlation matrix may be ill-conditioned, condition number: {%1.2e}" % (condn), RuntimeWarning)
485        chol = np.linalg.cholesky(corr)
486        chol_inv = scipy.linalg.solve_triangular(chol, covdiag, lower=True)
487
488        def chisqfunc_corr(p):
489            model = func(p, x)
490            chisq = anp.sum(anp.dot(chol_inv, (y_f - model)) ** 2)
491            return chisq
492
493    def chisqfunc(p):
494        model = func(p, x)
495        chisq = anp.sum(((y_f - model) / dy_f) ** 2)
496        return chisq
497
498    output.method = kwargs.get('method', 'Levenberg-Marquardt')
499    if not silent:
500        print('Method:', output.method)
501
502    if output.method != 'Levenberg-Marquardt':
503        if output.method == 'migrad':
504            fit_result = iminuit.minimize(chisqfunc, x0, tol=1e-4)  # Stopping criterion 0.002 * tol * errordef
505            if kwargs.get('correlated_fit') is True:
506                fit_result = iminuit.minimize(chisqfunc_corr, fit_result.x, tol=1e-4)  # Stopping criterion 0.002 * tol * errordef
507            output.iterations = fit_result.nfev
508        else:
509            fit_result = scipy.optimize.minimize(chisqfunc, x0, method=kwargs.get('method'), tol=1e-12)
510            if kwargs.get('correlated_fit') is True:
511                fit_result = scipy.optimize.minimize(chisqfunc_corr, fit_result.x, method=kwargs.get('method'), tol=1e-12)
512            output.iterations = fit_result.nit
513
514        chisquare = fit_result.fun
515
516    else:
517        if kwargs.get('correlated_fit') is True:
518            def chisqfunc_residuals_corr(p):
519                model = func(p, x)
520                chisq = anp.dot(chol_inv, (y_f - model))
521                return chisq
522
523        def chisqfunc_residuals(p):
524            model = func(p, x)
525            chisq = ((y_f - model) / dy_f)
526            return chisq
527
528        fit_result = scipy.optimize.least_squares(chisqfunc_residuals, x0, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15)
529        if kwargs.get('correlated_fit') is True:
530            fit_result = scipy.optimize.least_squares(chisqfunc_residuals_corr, fit_result.x, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15)
531
532        chisquare = np.sum(fit_result.fun ** 2)
533        if kwargs.get('correlated_fit') is True:
534            assert np.isclose(chisquare, chisqfunc_corr(fit_result.x), atol=1e-14)
535        else:
536            assert np.isclose(chisquare, chisqfunc(fit_result.x), atol=1e-14)
537
538        output.iterations = fit_result.nfev
539
540    if not fit_result.success:
541        raise Exception('The minimization procedure did not converge.')
542
543    if x.shape[-1] - n_parms > 0:
544        output.chisquare_by_dof = chisquare / (x.shape[-1] - n_parms)
545    else:
546        output.chisquare_by_dof = float('nan')
547
548    output.message = fit_result.message
549    if not silent:
550        print(fit_result.message)
551        print('chisquare/d.o.f.:', output.chisquare_by_dof)
552
553    if kwargs.get('expected_chisquare') is True:
554        if kwargs.get('correlated_fit') is not True:
555            W = np.diag(1 / np.asarray(dy_f))
556            cov = covariance(y)
557            A = W @ jacobian(func)(fit_result.x, x)
558            P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T
559            expected_chisquare = np.trace((np.identity(x.shape[-1]) - P_phi) @ W @ cov @ W)
560            output.chisquare_by_expected_chisquare = chisquare / expected_chisquare
561            if not silent:
562                print('chisquare/expected_chisquare:',
563                      output.chisquare_by_expected_chisquare)
564
565    fitp = fit_result.x
566    try:
567        if kwargs.get('correlated_fit') is True:
568            hess = jacobian(jacobian(chisqfunc_corr))(fitp)
569        else:
570            hess = jacobian(jacobian(chisqfunc))(fitp)
571    except TypeError:
572        raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None
573
574    if kwargs.get('correlated_fit') is True:
575        def chisqfunc_compact(d):
576            model = func(d[:n_parms], x)
577            chisq = anp.sum(anp.dot(chol_inv, (d[n_parms:] - model)) ** 2)
578            return chisq
579
580    else:
581        def chisqfunc_compact(d):
582            model = func(d[:n_parms], x)
583            chisq = anp.sum(((d[n_parms:] - model) / dy_f) ** 2)
584            return chisq
585
586    jac_jac = jacobian(jacobian(chisqfunc_compact))(np.concatenate((fitp, y_f)))
587
588    # Compute hess^{-1} @ jac_jac[:n_parms, n_parms:] using LAPACK dgesv
589    try:
590        deriv = -scipy.linalg.solve(hess, jac_jac[:n_parms, n_parms:])
591    except np.linalg.LinAlgError:
592        raise Exception("Cannot invert hessian matrix.")
593
594    result = []
595    for i in range(n_parms):
596        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])))
597
598    output.fit_parameters = result
599
600    output.chisquare = chisquare
601    output.dof = x.shape[-1] - n_parms
602    output.p_value = 1 - scipy.stats.chi2.cdf(output.chisquare, output.dof)
603
604    if kwargs.get('resplot') is True:
605        residual_plot(x, y, func, result)
606
607    if kwargs.get('qqplot') is True:
608        qqplot(x, y, func, result)
609
610    return output
611
612
613def fit_lin(x, y, **kwargs):
614    """Performs a linear fit to y = n + m * x and returns two Obs n, m.
615
616    Parameters
617    ----------
618    x : list
619        Can either be a list of floats in which case no xerror is assumed, or
620        a list of Obs, where the dvalues of the Obs are used as xerror for the fit.
621    y : list
622        List of Obs, the dvalues of the Obs are used as yerror for the fit.
623    """
624
625    def f(a, x):
626        y = a[0] + a[1] * x
627        return y
628
629    if all(isinstance(n, Obs) for n in x):
630        out = total_least_squares(x, y, f, **kwargs)
631        return out.fit_parameters
632    elif all(isinstance(n, float) or isinstance(n, int) for n in x) or isinstance(x, np.ndarray):
633        out = least_squares(x, y, f, **kwargs)
634        return out.fit_parameters
635    else:
636        raise Exception('Unsupported types for x')
637
638
639def qqplot(x, o_y, func, p):
640    """Generates a quantile-quantile plot of the fit result which can be used to
641       check if the residuals of the fit are gaussian distributed.
642    """
643
644    residuals = []
645    for i_x, i_y in zip(x, o_y):
646        residuals.append((i_y - func(p, i_x)) / i_y.dvalue)
647    residuals = sorted(residuals)
648    my_y = [o.value for o in residuals]
649    probplot = scipy.stats.probplot(my_y)
650    my_x = probplot[0][0]
651    plt.figure(figsize=(8, 8 / 1.618))
652    plt.errorbar(my_x, my_y, fmt='o')
653    fit_start = my_x[0]
654    fit_stop = my_x[-1]
655    samples = np.arange(fit_start, fit_stop, 0.01)
656    plt.plot(samples, samples, 'k--', zorder=11, label='Standard normal distribution')
657    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='-')
658
659    plt.xlabel('Theoretical quantiles')
660    plt.ylabel('Ordered Values')
661    plt.legend()
662    plt.draw()
663
664
665def residual_plot(x, y, func, fit_res):
666    """ Generates a plot which compares the fit to the data and displays the corresponding residuals"""
667    sorted_x = sorted(x)
668    xstart = sorted_x[0] - 0.5 * (sorted_x[1] - sorted_x[0])
669    xstop = sorted_x[-1] + 0.5 * (sorted_x[-1] - sorted_x[-2])
670    x_samples = np.arange(xstart, xstop + 0.01, 0.01)
671
672    plt.figure(figsize=(8, 8 / 1.618))
673    gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
674    ax0 = plt.subplot(gs[0])
675    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')
676    ax0.plot(x_samples, func([o.value for o in fit_res], x_samples), label='Fit', zorder=10, ls='-', ms=0)
677    ax0.set_xticklabels([])
678    ax0.set_xlim([xstart, xstop])
679    ax0.set_xticklabels([])
680    ax0.legend()
681
682    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])
683    ax1 = plt.subplot(gs[1])
684    ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
685    ax1.tick_params(direction='out')
686    ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
687    ax1.axhline(y=0.0, ls='--', color='k', marker=" ")
688    ax1.fill_between(x_samples, -1.0, 1.0, alpha=0.1, facecolor='k')
689    ax1.set_xlim([xstart, xstop])
690    ax1.set_ylabel('Residuals')
691    plt.subplots_adjust(wspace=None, hspace=None)
692    plt.draw()
693
694
695def error_band(x, func, beta):
696    """Returns the error band for an array of sample values x, for given fit function func with optimized parameters beta."""
697    cov = covariance(beta)
698    if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps):
699        warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning)
700
701    deriv = []
702    for i, item in enumerate(x):
703        deriv.append(np.array(egrad(func)([o.value for o in beta], item)))
704
705    err = []
706    for i, item in enumerate(x):
707        err.append(np.sqrt(deriv[i] @ cov @ deriv[i]))
708    err = np.array(err)
709
710    return err
711
712
713def ks_test(objects=None):
714    """Performs a Kolmogorov–Smirnov test for the p-values of all fit object.
715
716    Parameters
717    ----------
718    objects : list
719        List of fit results to include in the analysis (optional).
720    """
721
722    if objects is None:
723        obs_list = []
724        for obj in gc.get_objects():
725            if isinstance(obj, Fit_result):
726                obs_list.append(obj)
727    else:
728        obs_list = objects
729
730    p_values = [o.p_value for o in obs_list]
731
732    bins = len(p_values)
733    x = np.arange(0, 1.001, 0.001)
734    plt.plot(x, x, 'k', zorder=1)
735    plt.xlim(0, 1)
736    plt.ylim(0, 1)
737    plt.xlabel('p-value')
738    plt.ylabel('Cumulative probability')
739    plt.title(str(bins) + ' p-values')
740
741    n = np.arange(1, bins + 1) / np.float64(bins)
742    Xs = np.sort(p_values)
743    plt.step(Xs, n)
744    diffs = n - Xs
745    loc_max_diff = np.argmax(np.abs(diffs))
746    loc = Xs[loc_max_diff]
747    plt.annotate('', xy=(loc, loc), xytext=(loc, loc + diffs[loc_max_diff]), arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
748    plt.draw()
749
750    print(scipy.stats.kstest(p_values, 'uniform'))