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