diff --git a/docs/pyerrors/fits.html b/docs/pyerrors/fits.html
index e84840a6..f7ef19ee 100644
--- a/docs/pyerrors/fits.html
+++ b/docs/pyerrors/fits.html
@@ -580,269 +580,268 @@
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 = np.linalg.inv(chol)
-481 chol_inv = np.dot(chol_inv, covdiag)
-482
-483 def chisqfunc_corr(p):
-484 model = func(p, x)
-485 chisq = anp.sum(anp.dot(chol_inv, (y_f - model)) ** 2)
-486 return chisq
-487
-488 def chisqfunc(p):
-489 model = func(p, x)
-490 chisq = anp.sum(((y_f - model) / dy_f) ** 2)
-491 return chisq
-492
-493 output.method = kwargs.get('method', 'Levenberg-Marquardt')
-494 if not silent:
-495 print('Method:', output.method)
-496
-497 if output.method != 'Levenberg-Marquardt':
-498 if output.method == 'migrad':
-499 fit_result = iminuit.minimize(chisqfunc, x0, tol=1e-4) # Stopping criterion 0.002 * tol * errordef
-500 if kwargs.get('correlated_fit') is True:
-501 fit_result = iminuit.minimize(chisqfunc_corr, fit_result.x, tol=1e-4) # Stopping criterion 0.002 * tol * errordef
-502 output.iterations = fit_result.nfev
-503 else:
-504 fit_result = scipy.optimize.minimize(chisqfunc, x0, method=kwargs.get('method'), tol=1e-12)
-505 if kwargs.get('correlated_fit') is True:
-506 fit_result = scipy.optimize.minimize(chisqfunc_corr, fit_result.x, method=kwargs.get('method'), tol=1e-12)
-507 output.iterations = fit_result.nit
-508
-509 chisquare = fit_result.fun
-510
-511 else:
-512 if kwargs.get('correlated_fit') is True:
-513 def chisqfunc_residuals_corr(p):
-514 model = func(p, x)
-515 chisq = anp.dot(chol_inv, (y_f - model))
-516 return chisq
-517
-518 def chisqfunc_residuals(p):
-519 model = func(p, x)
-520 chisq = ((y_f - model) / dy_f)
-521 return chisq
-522
-523 fit_result = scipy.optimize.least_squares(chisqfunc_residuals, x0, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15)
-524 if kwargs.get('correlated_fit') is True:
-525 fit_result = scipy.optimize.least_squares(chisqfunc_residuals_corr, fit_result.x, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15)
-526
-527 chisquare = np.sum(fit_result.fun ** 2)
-528 if kwargs.get('correlated_fit') is True:
-529 assert np.isclose(chisquare, chisqfunc_corr(fit_result.x), atol=1e-14)
-530 else:
-531 assert np.isclose(chisquare, chisqfunc(fit_result.x), atol=1e-14)
-532
-533 output.iterations = fit_result.nfev
-534
-535 if not fit_result.success:
-536 raise Exception('The minimization procedure did not converge.')
-537
-538 if x.shape[-1] - n_parms > 0:
-539 output.chisquare_by_dof = chisquare / (x.shape[-1] - n_parms)
-540 else:
-541 output.chisquare_by_dof = float('nan')
-542
-543 output.message = fit_result.message
-544 if not silent:
-545 print(fit_result.message)
-546 print('chisquare/d.o.f.:', output.chisquare_by_dof)
-547
-548 if kwargs.get('expected_chisquare') is True:
-549 if kwargs.get('correlated_fit') is not True:
-550 W = np.diag(1 / np.asarray(dy_f))
-551 cov = covariance(y)
-552 A = W @ jacobian(func)(fit_result.x, x)
-553 P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T
-554 expected_chisquare = np.trace((np.identity(x.shape[-1]) - P_phi) @ W @ cov @ W)
-555 output.chisquare_by_expected_chisquare = chisquare / expected_chisquare
-556 if not silent:
-557 print('chisquare/expected_chisquare:',
-558 output.chisquare_by_expected_chisquare)
-559
-560 fitp = fit_result.x
-561 try:
-562 hess = jacobian(jacobian(chisqfunc))(fitp)
-563 except TypeError:
-564 raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None
-565
-566 if kwargs.get('correlated_fit') is True:
-567 def chisqfunc_compact(d):
-568 model = func(d[:n_parms], x)
-569 chisq = anp.sum(anp.dot(chol_inv, (d[n_parms:] - model)) ** 2)
-570 return chisq
-571
-572 else:
-573 def chisqfunc_compact(d):
-574 model = func(d[:n_parms], x)
-575 chisq = anp.sum(((d[n_parms:] - model) / dy_f) ** 2)
-576 return chisq
-577
-578 jac_jac = jacobian(jacobian(chisqfunc_compact))(np.concatenate((fitp, y_f)))
-579
-580 # Compute hess^{-1} @ jac_jac[:n_parms, n_parms:] using LAPACK dgesv
-581 try:
-582 deriv = -scipy.linalg.solve(hess, jac_jac[:n_parms, n_parms:])
-583 except np.linalg.LinAlgError:
-584 raise Exception("Cannot invert hessian matrix.")
-585
-586 result = []
-587 for i in range(n_parms):
-588 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])))
-589
-590 output.fit_parameters = result
-591
-592 output.chisquare = chisquare
-593 output.dof = x.shape[-1] - n_parms
-594 output.p_value = 1 - chi2.cdf(output.chisquare, output.dof)
-595
-596 if kwargs.get('resplot') is True:
-597 residual_plot(x, y, func, result)
-598
-599 if kwargs.get('qqplot') is True:
-600 qqplot(x, y, func, result)
-601
-602 return output
+480 chol_inv = scipy.linalg.solve(chol, covdiag)
+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 hess = jacobian(jacobian(chisqfunc))(fitp)
+562 except TypeError:
+563 raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None
+564
+565 if kwargs.get('correlated_fit') is True:
+566 def chisqfunc_compact(d):
+567 model = func(d[:n_parms], x)
+568 chisq = anp.sum(anp.dot(chol_inv, (d[n_parms:] - model)) ** 2)
+569 return chisq
+570
+571 else:
+572 def chisqfunc_compact(d):
+573 model = func(d[:n_parms], x)
+574 chisq = anp.sum(((d[n_parms:] - model) / dy_f) ** 2)
+575 return chisq
+576
+577 jac_jac = jacobian(jacobian(chisqfunc_compact))(np.concatenate((fitp, y_f)))
+578
+579 # Compute hess^{-1} @ jac_jac[:n_parms, n_parms:] using LAPACK dgesv
+580 try:
+581 deriv = -scipy.linalg.solve(hess, jac_jac[:n_parms, n_parms:])
+582 except np.linalg.LinAlgError:
+583 raise Exception("Cannot invert hessian matrix.")
+584
+585 result = []
+586 for i in range(n_parms):
+587 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])))
+588
+589 output.fit_parameters = result
+590
+591 output.chisquare = chisquare
+592 output.dof = x.shape[-1] - n_parms
+593 output.p_value = 1 - chi2.cdf(output.chisquare, output.dof)
+594
+595 if kwargs.get('resplot') is True:
+596 residual_plot(x, y, func, result)
+597
+598 if kwargs.get('qqplot') is True:
+599 qqplot(x, y, func, result)
+600
+601 return output
+602
603
-604
-605def fit_lin(x, y, **kwargs):
-606 """Performs a linear fit to y = n + m * x and returns two Obs n, m.
-607
-608 Parameters
-609 ----------
-610 x : list
-611 Can either be a list of floats in which case no xerror is assumed, or
-612 a list of Obs, where the dvalues of the Obs are used as xerror for the fit.
-613 y : list
-614 List of Obs, the dvalues of the Obs are used as yerror for the fit.
-615 """
-616
-617 def f(a, x):
-618 y = a[0] + a[1] * x
-619 return y
-620
-621 if all(isinstance(n, Obs) for n in x):
-622 out = total_least_squares(x, y, f, **kwargs)
-623 return out.fit_parameters
-624 elif all(isinstance(n, float) or isinstance(n, int) for n in x) or isinstance(x, np.ndarray):
-625 out = least_squares(x, y, f, **kwargs)
-626 return out.fit_parameters
-627 else:
-628 raise Exception('Unsupported types for x')
+604def fit_lin(x, y, **kwargs):
+605 """Performs a linear fit to y = n + m * x and returns two Obs n, m.
+606
+607 Parameters
+608 ----------
+609 x : list
+610 Can either be a list of floats in which case no xerror is assumed, or
+611 a list of Obs, where the dvalues of the Obs are used as xerror for the fit.
+612 y : list
+613 List of Obs, the dvalues of the Obs are used as yerror for the fit.
+614 """
+615
+616 def f(a, x):
+617 y = a[0] + a[1] * x
+618 return y
+619
+620 if all(isinstance(n, Obs) for n in x):
+621 out = total_least_squares(x, y, f, **kwargs)
+622 return out.fit_parameters
+623 elif all(isinstance(n, float) or isinstance(n, int) for n in x) or isinstance(x, np.ndarray):
+624 out = least_squares(x, y, f, **kwargs)
+625 return out.fit_parameters
+626 else:
+627 raise Exception('Unsupported types for x')
+628
629
-630
-631def qqplot(x, o_y, func, p):
-632 """Generates a quantile-quantile plot of the fit result which can be used to
-633 check if the residuals of the fit are gaussian distributed.
-634 """
-635
-636 residuals = []
-637 for i_x, i_y in zip(x, o_y):
-638 residuals.append((i_y - func(p, i_x)) / i_y.dvalue)
-639 residuals = sorted(residuals)
-640 my_y = [o.value for o in residuals]
-641 probplot = scipy.stats.probplot(my_y)
-642 my_x = probplot[0][0]
-643 plt.figure(figsize=(8, 8 / 1.618))
-644 plt.errorbar(my_x, my_y, fmt='o')
-645 fit_start = my_x[0]
-646 fit_stop = my_x[-1]
-647 samples = np.arange(fit_start, fit_stop, 0.01)
-648 plt.plot(samples, samples, 'k--', zorder=11, label='Standard normal distribution')
-649 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='-')
-650
-651 plt.xlabel('Theoretical quantiles')
-652 plt.ylabel('Ordered Values')
-653 plt.legend()
-654 plt.draw()
+630def qqplot(x, o_y, func, p):
+631 """Generates a quantile-quantile plot of the fit result which can be used to
+632 check if the residuals of the fit are gaussian distributed.
+633 """
+634
+635 residuals = []
+636 for i_x, i_y in zip(x, o_y):
+637 residuals.append((i_y - func(p, i_x)) / i_y.dvalue)
+638 residuals = sorted(residuals)
+639 my_y = [o.value for o in residuals]
+640 probplot = scipy.stats.probplot(my_y)
+641 my_x = probplot[0][0]
+642 plt.figure(figsize=(8, 8 / 1.618))
+643 plt.errorbar(my_x, my_y, fmt='o')
+644 fit_start = my_x[0]
+645 fit_stop = my_x[-1]
+646 samples = np.arange(fit_start, fit_stop, 0.01)
+647 plt.plot(samples, samples, 'k--', zorder=11, label='Standard normal distribution')
+648 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='-')
+649
+650 plt.xlabel('Theoretical quantiles')
+651 plt.ylabel('Ordered Values')
+652 plt.legend()
+653 plt.draw()
+654
655
-656
-657def residual_plot(x, y, func, fit_res):
-658 """ Generates a plot which compares the fit to the data and displays the corresponding residuals"""
-659 sorted_x = sorted(x)
-660 xstart = sorted_x[0] - 0.5 * (sorted_x[1] - sorted_x[0])
-661 xstop = sorted_x[-1] + 0.5 * (sorted_x[-1] - sorted_x[-2])
-662 x_samples = np.arange(xstart, xstop + 0.01, 0.01)
-663
-664 plt.figure(figsize=(8, 8 / 1.618))
-665 gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
-666 ax0 = plt.subplot(gs[0])
-667 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')
-668 ax0.plot(x_samples, func([o.value for o in fit_res], x_samples), label='Fit', zorder=10, ls='-', ms=0)
-669 ax0.set_xticklabels([])
-670 ax0.set_xlim([xstart, xstop])
-671 ax0.set_xticklabels([])
-672 ax0.legend()
-673
-674 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])
-675 ax1 = plt.subplot(gs[1])
-676 ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
-677 ax1.tick_params(direction='out')
-678 ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
-679 ax1.axhline(y=0.0, ls='--', color='k', marker=" ")
-680 ax1.fill_between(x_samples, -1.0, 1.0, alpha=0.1, facecolor='k')
-681 ax1.set_xlim([xstart, xstop])
-682 ax1.set_ylabel('Residuals')
-683 plt.subplots_adjust(wspace=None, hspace=None)
-684 plt.draw()
+656def residual_plot(x, y, func, fit_res):
+657 """ Generates a plot which compares the fit to the data and displays the corresponding residuals"""
+658 sorted_x = sorted(x)
+659 xstart = sorted_x[0] - 0.5 * (sorted_x[1] - sorted_x[0])
+660 xstop = sorted_x[-1] + 0.5 * (sorted_x[-1] - sorted_x[-2])
+661 x_samples = np.arange(xstart, xstop + 0.01, 0.01)
+662
+663 plt.figure(figsize=(8, 8 / 1.618))
+664 gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
+665 ax0 = plt.subplot(gs[0])
+666 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')
+667 ax0.plot(x_samples, func([o.value for o in fit_res], x_samples), label='Fit', zorder=10, ls='-', ms=0)
+668 ax0.set_xticklabels([])
+669 ax0.set_xlim([xstart, xstop])
+670 ax0.set_xticklabels([])
+671 ax0.legend()
+672
+673 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])
+674 ax1 = plt.subplot(gs[1])
+675 ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
+676 ax1.tick_params(direction='out')
+677 ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
+678 ax1.axhline(y=0.0, ls='--', color='k', marker=" ")
+679 ax1.fill_between(x_samples, -1.0, 1.0, alpha=0.1, facecolor='k')
+680 ax1.set_xlim([xstart, xstop])
+681 ax1.set_ylabel('Residuals')
+682 plt.subplots_adjust(wspace=None, hspace=None)
+683 plt.draw()
+684
685
-686
-687def error_band(x, func, beta):
-688 """Returns the error band for an array of sample values x, for given fit function func with optimized parameters beta."""
-689 cov = covariance(beta)
-690 if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps):
-691 warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning)
-692
-693 deriv = []
-694 for i, item in enumerate(x):
-695 deriv.append(np.array(egrad(func)([o.value for o in beta], item)))
-696
-697 err = []
-698 for i, item in enumerate(x):
-699 err.append(np.sqrt(deriv[i] @ cov @ deriv[i]))
-700 err = np.array(err)
-701
-702 return err
+686def error_band(x, func, beta):
+687 """Returns the error band for an array of sample values x, for given fit function func with optimized parameters beta."""
+688 cov = covariance(beta)
+689 if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps):
+690 warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning)
+691
+692 deriv = []
+693 for i, item in enumerate(x):
+694 deriv.append(np.array(egrad(func)([o.value for o in beta], item)))
+695
+696 err = []
+697 for i, item in enumerate(x):
+698 err.append(np.sqrt(deriv[i] @ cov @ deriv[i]))
+699 err = np.array(err)
+700
+701 return err
+702
703
-704
-705def ks_test(objects=None):
-706 """Performs a Kolmogorov–Smirnov test for the p-values of all fit object.
-707
-708 Parameters
-709 ----------
-710 objects : list
-711 List of fit results to include in the analysis (optional).
-712 """
-713
-714 if objects is None:
-715 obs_list = []
-716 for obj in gc.get_objects():
-717 if isinstance(obj, Fit_result):
-718 obs_list.append(obj)
-719 else:
-720 obs_list = objects
-721
-722 p_values = [o.p_value for o in obs_list]
-723
-724 bins = len(p_values)
-725 x = np.arange(0, 1.001, 0.001)
-726 plt.plot(x, x, 'k', zorder=1)
-727 plt.xlim(0, 1)
-728 plt.ylim(0, 1)
-729 plt.xlabel('p-value')
-730 plt.ylabel('Cumulative probability')
-731 plt.title(str(bins) + ' p-values')
-732
-733 n = np.arange(1, bins + 1) / np.float64(bins)
-734 Xs = np.sort(p_values)
-735 plt.step(Xs, n)
-736 diffs = n - Xs
-737 loc_max_diff = np.argmax(np.abs(diffs))
-738 loc = Xs[loc_max_diff]
-739 plt.annotate('', xy=(loc, loc), xytext=(loc, loc + diffs[loc_max_diff]), arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
-740 plt.draw()
-741
-742 print(scipy.stats.kstest(p_values, 'uniform'))
+704def ks_test(objects=None):
+705 """Performs a Kolmogorov–Smirnov test for the p-values of all fit object.
+706
+707 Parameters
+708 ----------
+709 objects : list
+710 List of fit results to include in the analysis (optional).
+711 """
+712
+713 if objects is None:
+714 obs_list = []
+715 for obj in gc.get_objects():
+716 if isinstance(obj, Fit_result):
+717 obs_list.append(obj)
+718 else:
+719 obs_list = objects
+720
+721 p_values = [o.p_value for o in obs_list]
+722
+723 bins = len(p_values)
+724 x = np.arange(0, 1.001, 0.001)
+725 plt.plot(x, x, 'k', zorder=1)
+726 plt.xlim(0, 1)
+727 plt.ylim(0, 1)
+728 plt.xlabel('p-value')
+729 plt.ylabel('Cumulative probability')
+730 plt.title(str(bins) + ' p-values')
+731
+732 n = np.arange(1, bins + 1) / np.float64(bins)
+733 Xs = np.sort(p_values)
+734 plt.step(Xs, n)
+735 diffs = n - Xs
+736 loc_max_diff = np.argmax(np.abs(diffs))
+737 loc = Xs[loc_max_diff]
+738 plt.annotate('', xy=(loc, loc), xytext=(loc, loc + diffs[loc_max_diff]), arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
+739 plt.draw()
+740
+741 print(scipy.stats.kstest(p_values, 'uniform'))
@@ -1354,30 +1353,30 @@ has to be calculated (default False).
- 606def fit_lin(x, y, **kwargs):
-607 """Performs a linear fit to y = n + m * x and returns two Obs n, m.
-608
-609 Parameters
-610 ----------
-611 x : list
-612 Can either be a list of floats in which case no xerror is assumed, or
-613 a list of Obs, where the dvalues of the Obs are used as xerror for the fit.
-614 y : list
-615 List of Obs, the dvalues of the Obs are used as yerror for the fit.
-616 """
-617
-618 def f(a, x):
-619 y = a[0] + a[1] * x
-620 return y
-621
-622 if all(isinstance(n, Obs) for n in x):
-623 out = total_least_squares(x, y, f, **kwargs)
-624 return out.fit_parameters
-625 elif all(isinstance(n, float) or isinstance(n, int) for n in x) or isinstance(x, np.ndarray):
-626 out = least_squares(x, y, f, **kwargs)
-627 return out.fit_parameters
-628 else:
-629 raise Exception('Unsupported types for x')
+ 605def fit_lin(x, y, **kwargs):
+606 """Performs a linear fit to y = n + m * x and returns two Obs n, m.
+607
+608 Parameters
+609 ----------
+610 x : list
+611 Can either be a list of floats in which case no xerror is assumed, or
+612 a list of Obs, where the dvalues of the Obs are used as xerror for the fit.
+613 y : list
+614 List of Obs, the dvalues of the Obs are used as yerror for the fit.
+615 """
+616
+617 def f(a, x):
+618 y = a[0] + a[1] * x
+619 return y
+620
+621 if all(isinstance(n, Obs) for n in x):
+622 out = total_least_squares(x, y, f, **kwargs)
+623 return out.fit_parameters
+624 elif all(isinstance(n, float) or isinstance(n, int) for n in x) or isinstance(x, np.ndarray):
+625 out = least_squares(x, y, f, **kwargs)
+626 return out.fit_parameters
+627 else:
+628 raise Exception('Unsupported types for x')
632def qqplot(x, o_y, func, p):
-633 """Generates a quantile-quantile plot of the fit result which can be used to
-634 check if the residuals of the fit are gaussian distributed.
-635 """
-636
-637 residuals = []
-638 for i_x, i_y in zip(x, o_y):
-639 residuals.append((i_y - func(p, i_x)) / i_y.dvalue)
-640 residuals = sorted(residuals)
-641 my_y = [o.value for o in residuals]
-642 probplot = scipy.stats.probplot(my_y)
-643 my_x = probplot[0][0]
-644 plt.figure(figsize=(8, 8 / 1.618))
-645 plt.errorbar(my_x, my_y, fmt='o')
-646 fit_start = my_x[0]
-647 fit_stop = my_x[-1]
-648 samples = np.arange(fit_start, fit_stop, 0.01)
-649 plt.plot(samples, samples, 'k--', zorder=11, label='Standard normal distribution')
-650 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='-')
-651
-652 plt.xlabel('Theoretical quantiles')
-653 plt.ylabel('Ordered Values')
-654 plt.legend()
-655 plt.draw()
+ 631def qqplot(x, o_y, func, p):
+632 """Generates a quantile-quantile plot of the fit result which can be used to
+633 check if the residuals of the fit are gaussian distributed.
+634 """
+635
+636 residuals = []
+637 for i_x, i_y in zip(x, o_y):
+638 residuals.append((i_y - func(p, i_x)) / i_y.dvalue)
+639 residuals = sorted(residuals)
+640 my_y = [o.value for o in residuals]
+641 probplot = scipy.stats.probplot(my_y)
+642 my_x = probplot[0][0]
+643 plt.figure(figsize=(8, 8 / 1.618))
+644 plt.errorbar(my_x, my_y, fmt='o')
+645 fit_start = my_x[0]
+646 fit_stop = my_x[-1]
+647 samples = np.arange(fit_start, fit_stop, 0.01)
+648 plt.plot(samples, samples, 'k--', zorder=11, label='Standard normal distribution')
+649 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='-')
+650
+651 plt.xlabel('Theoretical quantiles')
+652 plt.ylabel('Ordered Values')
+653 plt.legend()
+654 plt.draw()
658def residual_plot(x, y, func, fit_res):
-659 """ Generates a plot which compares the fit to the data and displays the corresponding residuals"""
-660 sorted_x = sorted(x)
-661 xstart = sorted_x[0] - 0.5 * (sorted_x[1] - sorted_x[0])
-662 xstop = sorted_x[-1] + 0.5 * (sorted_x[-1] - sorted_x[-2])
-663 x_samples = np.arange(xstart, xstop + 0.01, 0.01)
-664
-665 plt.figure(figsize=(8, 8 / 1.618))
-666 gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
-667 ax0 = plt.subplot(gs[0])
-668 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')
-669 ax0.plot(x_samples, func([o.value for o in fit_res], x_samples), label='Fit', zorder=10, ls='-', ms=0)
-670 ax0.set_xticklabels([])
-671 ax0.set_xlim([xstart, xstop])
-672 ax0.set_xticklabels([])
-673 ax0.legend()
-674
-675 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])
-676 ax1 = plt.subplot(gs[1])
-677 ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
-678 ax1.tick_params(direction='out')
-679 ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
-680 ax1.axhline(y=0.0, ls='--', color='k', marker=" ")
-681 ax1.fill_between(x_samples, -1.0, 1.0, alpha=0.1, facecolor='k')
-682 ax1.set_xlim([xstart, xstop])
-683 ax1.set_ylabel('Residuals')
-684 plt.subplots_adjust(wspace=None, hspace=None)
-685 plt.draw()
+ 657def residual_plot(x, y, func, fit_res):
+658 """ Generates a plot which compares the fit to the data and displays the corresponding residuals"""
+659 sorted_x = sorted(x)
+660 xstart = sorted_x[0] - 0.5 * (sorted_x[1] - sorted_x[0])
+661 xstop = sorted_x[-1] + 0.5 * (sorted_x[-1] - sorted_x[-2])
+662 x_samples = np.arange(xstart, xstop + 0.01, 0.01)
+663
+664 plt.figure(figsize=(8, 8 / 1.618))
+665 gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
+666 ax0 = plt.subplot(gs[0])
+667 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')
+668 ax0.plot(x_samples, func([o.value for o in fit_res], x_samples), label='Fit', zorder=10, ls='-', ms=0)
+669 ax0.set_xticklabels([])
+670 ax0.set_xlim([xstart, xstop])
+671 ax0.set_xticklabels([])
+672 ax0.legend()
+673
+674 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])
+675 ax1 = plt.subplot(gs[1])
+676 ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
+677 ax1.tick_params(direction='out')
+678 ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
+679 ax1.axhline(y=0.0, ls='--', color='k', marker=" ")
+680 ax1.fill_between(x_samples, -1.0, 1.0, alpha=0.1, facecolor='k')
+681 ax1.set_xlim([xstart, xstop])
+682 ax1.set_ylabel('Residuals')
+683 plt.subplots_adjust(wspace=None, hspace=None)
+684 plt.draw()
688def error_band(x, func, beta):
-689 """Returns the error band for an array of sample values x, for given fit function func with optimized parameters beta."""
-690 cov = covariance(beta)
-691 if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps):
-692 warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning)
-693
-694 deriv = []
-695 for i, item in enumerate(x):
-696 deriv.append(np.array(egrad(func)([o.value for o in beta], item)))
-697
-698 err = []
-699 for i, item in enumerate(x):
-700 err.append(np.sqrt(deriv[i] @ cov @ deriv[i]))
-701 err = np.array(err)
-702
-703 return err
+ 687def error_band(x, func, beta):
+688 """Returns the error band for an array of sample values x, for given fit function func with optimized parameters beta."""
+689 cov = covariance(beta)
+690 if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps):
+691 warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning)
+692
+693 deriv = []
+694 for i, item in enumerate(x):
+695 deriv.append(np.array(egrad(func)([o.value for o in beta], item)))
+696
+697 err = []
+698 for i, item in enumerate(x):
+699 err.append(np.sqrt(deriv[i] @ cov @ deriv[i]))
+700 err = np.array(err)
+701
+702 return err
706def ks_test(objects=None):
-707 """Performs a Kolmogorov–Smirnov test for the p-values of all fit object.
-708
-709 Parameters
-710 ----------
-711 objects : list
-712 List of fit results to include in the analysis (optional).
-713 """
-714
-715 if objects is None:
-716 obs_list = []
-717 for obj in gc.get_objects():
-718 if isinstance(obj, Fit_result):
-719 obs_list.append(obj)
-720 else:
-721 obs_list = objects
-722
-723 p_values = [o.p_value for o in obs_list]
-724
-725 bins = len(p_values)
-726 x = np.arange(0, 1.001, 0.001)
-727 plt.plot(x, x, 'k', zorder=1)
-728 plt.xlim(0, 1)
-729 plt.ylim(0, 1)
-730 plt.xlabel('p-value')
-731 plt.ylabel('Cumulative probability')
-732 plt.title(str(bins) + ' p-values')
-733
-734 n = np.arange(1, bins + 1) / np.float64(bins)
-735 Xs = np.sort(p_values)
-736 plt.step(Xs, n)
-737 diffs = n - Xs
-738 loc_max_diff = np.argmax(np.abs(diffs))
-739 loc = Xs[loc_max_diff]
-740 plt.annotate('', xy=(loc, loc), xytext=(loc, loc + diffs[loc_max_diff]), arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
-741 plt.draw()
-742
-743 print(scipy.stats.kstest(p_values, 'uniform'))
+ 705def ks_test(objects=None):
+706 """Performs a Kolmogorov–Smirnov test for the p-values of all fit object.
+707
+708 Parameters
+709 ----------
+710 objects : list
+711 List of fit results to include in the analysis (optional).
+712 """
+713
+714 if objects is None:
+715 obs_list = []
+716 for obj in gc.get_objects():
+717 if isinstance(obj, Fit_result):
+718 obs_list.append(obj)
+719 else:
+720 obs_list = objects
+721
+722 p_values = [o.p_value for o in obs_list]
+723
+724 bins = len(p_values)
+725 x = np.arange(0, 1.001, 0.001)
+726 plt.plot(x, x, 'k', zorder=1)
+727 plt.xlim(0, 1)
+728 plt.ylim(0, 1)
+729 plt.xlabel('p-value')
+730 plt.ylabel('Cumulative probability')
+731 plt.title(str(bins) + ' p-values')
+732
+733 n = np.arange(1, bins + 1) / np.float64(bins)
+734 Xs = np.sort(p_values)
+735 plt.step(Xs, n)
+736 diffs = n - Xs
+737 loc_max_diff = np.argmax(np.abs(diffs))
+738 loc = Xs[loc_max_diff]
+739 plt.annotate('', xy=(loc, loc), xytext=(loc, loc + diffs[loc_max_diff]), arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
+740 plt.draw()
+741
+742 print(scipy.stats.kstest(p_values, 'uniform'))