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