diff --git a/docs/pyerrors/fits.html b/docs/pyerrors/fits.html
index c2d1aba5..0df35391 100644
--- a/docs/pyerrors/fits.html
+++ b/docs/pyerrors/fits.html
@@ -688,337 +688,323 @@
580 else:
581 x0 = [0.1] * n_parms
582
-583 if kwargs.get('correlated_fit') is True:
-584 corr = covariance(y_all, correlation=True, **kwargs)
-585 covdiag = np.diag(1 / np.asarray(dy_f))
-586 condn = np.linalg.cond(corr)
-587 if condn > 0.1 / np.finfo(float).eps:
-588 raise Exception(f"Cannot invert correlation matrix as its condition number exceeds machine precision ({condn:1.2e})")
-589 if condn > 1e13:
-590 warnings.warn("Correlation matrix may be ill-conditioned, condition number: {%1.2e}" % (condn), RuntimeWarning)
-591 chol = np.linalg.cholesky(corr)
-592 chol_inv = scipy.linalg.solve_triangular(chol, covdiag, lower=True)
-593
-594 def chisqfunc_corr(p):
-595 model = np.concatenate([np.array(funcd[key](p, np.asarray(xd[key]))).reshape(-1) for key in key_ls])
-596 chisq = anp.sum(anp.dot(chol_inv, (y_f - model)) ** 2)
-597 return chisq
-598
-599 def chisqfunc(p):
-600 func_list = np.concatenate([[funcd[k]] * len(xd[k]) for k in key_ls])
-601 model = anp.array([func_list[i](p, x_all[i]) for i in range(len(x_all))])
-602 chisq = anp.sum(((y_f - model) / dy_f) ** 2)
-603 return chisq
+583 def general_chisqfunc_uncorr(p, ivars):
+584 model = anp.concatenate([anp.array(funcd[key](p, anp.asarray(xd[key]))).reshape(-1) for key in key_ls])
+585 return ((ivars - model) / dy_f)
+586
+587 def chisqfunc_uncorr(p):
+588 return anp.sum(general_chisqfunc_uncorr(p, y_f) ** 2)
+589
+590 if kwargs.get('correlated_fit') is True:
+591 corr = covariance(y_all, correlation=True, **kwargs)
+592 covdiag = np.diag(1 / np.asarray(dy_f))
+593 condn = np.linalg.cond(corr)
+594 if condn > 0.1 / np.finfo(float).eps:
+595 raise Exception(f"Cannot invert correlation matrix as its condition number exceeds machine precision ({condn:1.2e})")
+596 if condn > 1e13:
+597 warnings.warn("Correlation matrix may be ill-conditioned, condition number: {%1.2e}" % (condn), RuntimeWarning)
+598 chol = np.linalg.cholesky(corr)
+599 chol_inv = scipy.linalg.solve_triangular(chol, covdiag, lower=True)
+600
+601 def general_chisqfunc(p, ivars):
+602 model = anp.concatenate([anp.array(funcd[key](p, anp.asarray(xd[key]))).reshape(-1) for key in key_ls])
+603 return anp.dot(chol_inv, (ivars - model))
604
-605 output.method = kwargs.get('method', 'Levenberg-Marquardt')
-606 if not silent:
-607 print('Method:', output.method)
-608
-609 if output.method != 'Levenberg-Marquardt':
-610 if output.method == 'migrad':
-611 tolerance = 1e-4 # default value of 1e-1 set by iminuit can be problematic
-612 if 'tol' in kwargs:
-613 tolerance = kwargs.get('tol')
-614 fit_result = iminuit.minimize(chisqfunc, x0, tol=tolerance) # Stopping criterion 0.002 * tol * errordef
-615 if kwargs.get('correlated_fit') is True:
-616 fit_result = iminuit.minimize(chisqfunc_corr, fit_result.x, tol=tolerance)
-617 output.iterations = fit_result.nfev
-618 else:
-619 tolerance = 1e-12
-620 if 'tol' in kwargs:
-621 tolerance = kwargs.get('tol')
-622 fit_result = scipy.optimize.minimize(chisqfunc, x0, method=kwargs.get('method'), tol=tolerance)
-623 if kwargs.get('correlated_fit') is True:
-624 fit_result = scipy.optimize.minimize(chisqfunc_corr, fit_result.x, method=kwargs.get('method'), tol=tolerance)
-625 output.iterations = fit_result.nit
-626
-627 chisquare = fit_result.fun
-628
-629 else:
-630 if kwargs.get('correlated_fit') is True:
-631 def chisqfunc_residuals_corr(p):
-632 model = np.concatenate([np.array(funcd[key](p, np.asarray(xd[key]))).reshape(-1) for key in key_ls])
-633 chisq = anp.dot(chol_inv, (y_f - model))
-634 return chisq
-635
-636 def chisqfunc_residuals(p):
-637 model = np.concatenate([np.array(funcd[key](p, np.asarray(xd[key]))).reshape(-1) for key in key_ls])
-638 chisq = ((y_f - model) / dy_f)
-639 return chisq
-640
-641 if 'tol' in kwargs:
-642 print('tol cannot be set for Levenberg-Marquardt')
-643
-644 fit_result = scipy.optimize.least_squares(chisqfunc_residuals, x0, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15)
-645 if kwargs.get('correlated_fit') is True:
-646 fit_result = scipy.optimize.least_squares(chisqfunc_residuals_corr, fit_result.x, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15)
+605 def chisqfunc(p):
+606 return anp.sum(general_chisqfunc(p, y_f) ** 2)
+607 else:
+608 general_chisqfunc = general_chisqfunc_uncorr
+609 chisqfunc = chisqfunc_uncorr
+610
+611 output.method = kwargs.get('method', 'Levenberg-Marquardt')
+612 if not silent:
+613 print('Method:', output.method)
+614
+615 if output.method != 'Levenberg-Marquardt':
+616 if output.method == 'migrad':
+617 tolerance = 1e-4 # default value of 1e-1 set by iminuit can be problematic
+618 if 'tol' in kwargs:
+619 tolerance = kwargs.get('tol')
+620 fit_result = iminuit.minimize(chisqfunc_uncorr, x0, tol=tolerance) # Stopping criterion 0.002 * tol * errordef
+621 if kwargs.get('correlated_fit') is True:
+622 fit_result = iminuit.minimize(chisqfunc, fit_result.x, tol=tolerance)
+623 output.iterations = fit_result.nfev
+624 else:
+625 tolerance = 1e-12
+626 if 'tol' in kwargs:
+627 tolerance = kwargs.get('tol')
+628 fit_result = scipy.optimize.minimize(chisqfunc_uncorr, x0, method=kwargs.get('method'), tol=tolerance)
+629 if kwargs.get('correlated_fit') is True:
+630 fit_result = scipy.optimize.minimize(chisqfunc, fit_result.x, method=kwargs.get('method'), tol=tolerance)
+631 output.iterations = fit_result.nit
+632
+633 chisquare = fit_result.fun
+634
+635 else:
+636 if 'tol' in kwargs:
+637 print('tol cannot be set for Levenberg-Marquardt')
+638
+639 def chisqfunc_residuals_uncorr(p):
+640 return general_chisqfunc_uncorr(p, y_f)
+641
+642 fit_result = scipy.optimize.least_squares(chisqfunc_residuals_uncorr, x0, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15)
+643 if kwargs.get('correlated_fit') is True:
+644
+645 def chisqfunc_residuals(p):
+646 return general_chisqfunc(p, y_f)
647
-648 chisquare = np.sum(fit_result.fun ** 2)
-649 if kwargs.get('correlated_fit') is True:
-650 assert np.isclose(chisquare, chisqfunc_corr(fit_result.x), atol=1e-14)
-651 else:
-652 assert np.isclose(chisquare, chisqfunc(fit_result.x), atol=1e-14)
-653
-654 output.iterations = fit_result.nfev
-655
-656 if not fit_result.success:
-657 raise Exception('The minimization procedure did not converge.')
-658
-659 if x_all.shape[-1] - n_parms > 0:
-660 output.chisquare = chisquare
-661 output.dof = x_all.shape[-1] - n_parms
-662 output.chisquare_by_dof = output.chisquare / output.dof
-663 output.p_value = 1 - scipy.stats.chi2.cdf(output.chisquare, output.dof)
-664 else:
-665 output.chisquare_by_dof = float('nan')
-666
-667 output.message = fit_result.message
-668 if not silent:
-669 print(fit_result.message)
-670 print('chisquare/d.o.f.:', output.chisquare_by_dof)
-671 print('fit parameters', fit_result.x)
-672
-673 def prepare_hat_matrix():
-674 hat_vector = []
-675 for key in key_ls:
-676 x_array = np.asarray(xd[key])
-677 if (len(x_array) != 0):
-678 hat_vector.append(jacobian(funcd[key])(fit_result.x, x_array))
-679 hat_vector = [item for sublist in hat_vector for item in sublist]
-680 return hat_vector
-681
-682 if kwargs.get('expected_chisquare') is True:
-683 if kwargs.get('correlated_fit') is not True:
-684 W = np.diag(1 / np.asarray(dy_f))
-685 cov = covariance(y_all)
-686 hat_vector = prepare_hat_matrix()
-687 A = W @ hat_vector
-688 P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T
-689 expected_chisquare = np.trace((np.identity(x_all.shape[-1]) - P_phi) @ W @ cov @ W)
-690 output.chisquare_by_expected_chisquare = output.chisquare / expected_chisquare
-691 if not silent:
-692 print('chisquare/expected_chisquare:', output.chisquare_by_expected_chisquare)
-693
-694 fitp = fit_result.x
-695 if np.any(np.asarray(dy_f) <= 0.0):
-696 raise Exception('No y errors available, run the gamma method first.')
-697
-698 try:
-699 if kwargs.get('correlated_fit') is True:
-700 hess = hessian(chisqfunc_corr)(fitp)
-701 else:
-702 hess = hessian(chisqfunc)(fitp)
-703 except TypeError:
-704 raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None
-705
-706 if kwargs.get('correlated_fit') is True:
-707 def chisqfunc_compact(d):
-708 func_list = np.concatenate([[funcd[k]] * len(xd[k]) for k in key_ls])
-709 model = anp.array([func_list[i](d[:n_parms], x_all[i]) for i in range(len(x_all))])
-710 chisq = anp.sum(anp.dot(chol_inv, (d[n_parms:] - model)) ** 2)
-711 return chisq
-712 else:
-713 def chisqfunc_compact(d):
-714 func_list = np.concatenate([[funcd[k]] * len(xd[k]) for k in key_ls])
-715 model = anp.array([func_list[i](d[:n_parms], x_all[i]) for i in range(len(x_all))])
-716 chisq = anp.sum(((d[n_parms:] - model) / dy_f) ** 2)
-717 return chisq
+648 fit_result = scipy.optimize.least_squares(chisqfunc_residuals, fit_result.x, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15)
+649
+650 chisquare = np.sum(fit_result.fun ** 2)
+651 assert np.isclose(chisquare, chisqfunc(fit_result.x), atol=1e-14)
+652
+653 output.iterations = fit_result.nfev
+654
+655 if not fit_result.success:
+656 raise Exception('The minimization procedure did not converge.')
+657
+658 if x_all.shape[-1] - n_parms > 0:
+659 output.chisquare = chisquare
+660 output.dof = x_all.shape[-1] - n_parms
+661 output.chisquare_by_dof = output.chisquare / output.dof
+662 output.p_value = 1 - scipy.stats.chi2.cdf(output.chisquare, output.dof)
+663 else:
+664 output.chisquare_by_dof = float('nan')
+665
+666 output.message = fit_result.message
+667 if not silent:
+668 print(fit_result.message)
+669 print('chisquare/d.o.f.:', output.chisquare_by_dof)
+670 print('fit parameters', fit_result.x)
+671
+672 def prepare_hat_matrix():
+673 hat_vector = []
+674 for key in key_ls:
+675 x_array = np.asarray(xd[key])
+676 if (len(x_array) != 0):
+677 hat_vector.append(jacobian(funcd[key])(fit_result.x, x_array))
+678 hat_vector = [item for sublist in hat_vector for item in sublist]
+679 return hat_vector
+680
+681 if kwargs.get('expected_chisquare') is True:
+682 if kwargs.get('correlated_fit') is not True:
+683 W = np.diag(1 / np.asarray(dy_f))
+684 cov = covariance(y_all)
+685 hat_vector = prepare_hat_matrix()
+686 A = W @ hat_vector
+687 P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T
+688 expected_chisquare = np.trace((np.identity(x_all.shape[-1]) - P_phi) @ W @ cov @ W)
+689 output.chisquare_by_expected_chisquare = output.chisquare / expected_chisquare
+690 if not silent:
+691 print('chisquare/expected_chisquare:', output.chisquare_by_expected_chisquare)
+692
+693 fitp = fit_result.x
+694 if np.any(np.asarray(dy_f) <= 0.0):
+695 raise Exception('No y errors available, run the gamma method first.')
+696
+697 try:
+698 hess = hessian(chisqfunc)(fitp)
+699 except TypeError:
+700 raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None
+701
+702 def chisqfunc_compact(d):
+703 return anp.sum(general_chisqfunc(d[:n_parms], d[n_parms:]) ** 2)
+704
+705 jac_jac_y = hessian(chisqfunc_compact)(np.concatenate((fitp, y_f)))
+706
+707 # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv
+708 try:
+709 deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms, n_parms:])
+710 except np.linalg.LinAlgError:
+711 raise Exception("Cannot invert hessian matrix.")
+712
+713 result = []
+714 for i in range(n_parms):
+715 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), man_grad=list(deriv_y[i])))
+716
+717 output.fit_parameters = result
718
-719 jac_jac_y = hessian(chisqfunc_compact)(np.concatenate((fitp, y_f)))
-720
-721 # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv
-722 try:
-723 deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms, n_parms:])
-724 except np.linalg.LinAlgError:
-725 raise Exception("Cannot invert hessian matrix.")
-726
-727 result = []
-728 for i in range(n_parms):
-729 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), man_grad=list(deriv_y[i])))
-730
-731 output.fit_parameters = result
+719 # Hotelling t-squared p-value for correlated fits.
+720 if kwargs.get('correlated_fit') is True:
+721 n_cov = np.min(np.vectorize(lambda x_all: x_all.N)(y_all))
+722 output.t2_p_value = 1 - scipy.stats.f.cdf((n_cov - output.dof) / (output.dof * (n_cov - 1)) * output.chisquare,
+723 output.dof, n_cov - output.dof)
+724
+725 if kwargs.get('resplot') is True:
+726 for key in key_ls:
+727 residual_plot(xd[key], yd[key], funcd[key], result, title=key)
+728
+729 if kwargs.get('qqplot') is True:
+730 for key in key_ls:
+731 qqplot(xd[key], yd[key], funcd[key], result, title=key)
732
-733 # Hotelling t-squared p-value for correlated fits.
-734 if kwargs.get('correlated_fit') is True:
-735 n_cov = np.min(np.vectorize(lambda x_all: x_all.N)(y_all))
-736 output.t2_p_value = 1 - scipy.stats.f.cdf((n_cov - output.dof) / (output.dof * (n_cov - 1)) * output.chisquare,
-737 output.dof, n_cov - output.dof)
+733 return output
+734
+735
+736def fit_lin(x, y, **kwargs):
+737 """Performs a linear fit to y = n + m * x and returns two Obs n, m.
738
-739 if kwargs.get('resplot') is True:
-740 for key in key_ls:
-741 residual_plot(xd[key], yd[key], funcd[key], result, title=key)
-742
-743 if kwargs.get('qqplot') is True:
-744 for key in key_ls:
-745 qqplot(xd[key], yd[key], funcd[key], result, title=key)
+739 Parameters
+740 ----------
+741 x : list
+742 Can either be a list of floats in which case no xerror is assumed, or
+743 a list of Obs, where the dvalues of the Obs are used as xerror for the fit.
+744 y : list
+745 List of Obs, the dvalues of the Obs are used as yerror for the fit.
746
-747 return output
-748
-749
-750def fit_lin(x, y, **kwargs):
-751 """Performs a linear fit to y = n + m * x and returns two Obs n, m.
+747 Returns
+748 -------
+749 fit_parameters : list[Obs]
+750 LIist of fitted observables.
+751 """
752
-753 Parameters
-754 ----------
-755 x : list
-756 Can either be a list of floats in which case no xerror is assumed, or
-757 a list of Obs, where the dvalues of the Obs are used as xerror for the fit.
-758 y : list
-759 List of Obs, the dvalues of the Obs are used as yerror for the fit.
-760
-761 Returns
-762 -------
-763 fit_parameters : list[Obs]
-764 LIist of fitted observables.
-765 """
+753 def f(a, x):
+754 y = a[0] + a[1] * x
+755 return y
+756
+757 if all(isinstance(n, Obs) for n in x):
+758 out = total_least_squares(x, y, f, **kwargs)
+759 return out.fit_parameters
+760 elif all(isinstance(n, float) or isinstance(n, int) for n in x) or isinstance(x, np.ndarray):
+761 out = least_squares(x, y, f, **kwargs)
+762 return out.fit_parameters
+763 else:
+764 raise Exception('Unsupported types for x')
+765
766
-767 def f(a, x):
-768 y = a[0] + a[1] * x
-769 return y
+767def qqplot(x, o_y, func, p, title=""):
+768 """Generates a quantile-quantile plot of the fit result which can be used to
+769 check if the residuals of the fit are gaussian distributed.
770
-771 if all(isinstance(n, Obs) for n in x):
-772 out = total_least_squares(x, y, f, **kwargs)
-773 return out.fit_parameters
-774 elif all(isinstance(n, float) or isinstance(n, int) for n in x) or isinstance(x, np.ndarray):
-775 out = least_squares(x, y, f, **kwargs)
-776 return out.fit_parameters
-777 else:
-778 raise Exception('Unsupported types for x')
-779
-780
-781def qqplot(x, o_y, func, p, title=""):
-782 """Generates a quantile-quantile plot of the fit result which can be used to
-783 check if the residuals of the fit are gaussian distributed.
-784
-785 Returns
-786 -------
-787 None
-788 """
-789
-790 residuals = []
-791 for i_x, i_y in zip(x, o_y):
-792 residuals.append((i_y - func(p, i_x)) / i_y.dvalue)
-793 residuals = sorted(residuals)
-794 my_y = [o.value for o in residuals]
-795 probplot = scipy.stats.probplot(my_y)
-796 my_x = probplot[0][0]
-797 plt.figure(figsize=(8, 8 / 1.618))
-798 plt.errorbar(my_x, my_y, fmt='o')
-799 fit_start = my_x[0]
-800 fit_stop = my_x[-1]
-801 samples = np.arange(fit_start, fit_stop, 0.01)
-802 plt.plot(samples, samples, 'k--', zorder=11, label='Standard normal distribution')
-803 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='-')
-804
-805 plt.xlabel('Theoretical quantiles')
-806 plt.ylabel('Ordered Values')
-807 plt.legend(title=title)
-808 plt.draw()
-809
+771 Returns
+772 -------
+773 None
+774 """
+775
+776 residuals = []
+777 for i_x, i_y in zip(x, o_y):
+778 residuals.append((i_y - func(p, i_x)) / i_y.dvalue)
+779 residuals = sorted(residuals)
+780 my_y = [o.value for o in residuals]
+781 probplot = scipy.stats.probplot(my_y)
+782 my_x = probplot[0][0]
+783 plt.figure(figsize=(8, 8 / 1.618))
+784 plt.errorbar(my_x, my_y, fmt='o')
+785 fit_start = my_x[0]
+786 fit_stop = my_x[-1]
+787 samples = np.arange(fit_start, fit_stop, 0.01)
+788 plt.plot(samples, samples, 'k--', zorder=11, label='Standard normal distribution')
+789 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='-')
+790
+791 plt.xlabel('Theoretical quantiles')
+792 plt.ylabel('Ordered Values')
+793 plt.legend(title=title)
+794 plt.draw()
+795
+796
+797def residual_plot(x, y, func, fit_res, title=""):
+798 """Generates a plot which compares the fit to the data and displays the corresponding residuals
+799
+800 For uncorrelated data the residuals are expected to be distributed ~N(0,1).
+801
+802 Returns
+803 -------
+804 None
+805 """
+806 sorted_x = sorted(x)
+807 xstart = sorted_x[0] - 0.5 * (sorted_x[1] - sorted_x[0])
+808 xstop = sorted_x[-1] + 0.5 * (sorted_x[-1] - sorted_x[-2])
+809 x_samples = np.arange(xstart, xstop + 0.01, 0.01)
810
-811def residual_plot(x, y, func, fit_res, title=""):
-812 """Generates a plot which compares the fit to the data and displays the corresponding residuals
-813
-814 For uncorrelated data the residuals are expected to be distributed ~N(0,1).
-815
-816 Returns
-817 -------
-818 None
-819 """
-820 sorted_x = sorted(x)
-821 xstart = sorted_x[0] - 0.5 * (sorted_x[1] - sorted_x[0])
-822 xstop = sorted_x[-1] + 0.5 * (sorted_x[-1] - sorted_x[-2])
-823 x_samples = np.arange(xstart, xstop + 0.01, 0.01)
-824
-825 plt.figure(figsize=(8, 8 / 1.618))
-826 gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
-827 ax0 = plt.subplot(gs[0])
-828 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')
-829 ax0.plot(x_samples, func([o.value for o in fit_res], x_samples), label='Fit', zorder=10, ls='-', ms=0)
-830 ax0.set_xticklabels([])
-831 ax0.set_xlim([xstart, xstop])
-832 ax0.set_xticklabels([])
-833 ax0.legend(title=title)
-834
-835 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])
-836 ax1 = plt.subplot(gs[1])
-837 ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
-838 ax1.tick_params(direction='out')
-839 ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
-840 ax1.axhline(y=0.0, ls='--', color='k', marker=" ")
-841 ax1.fill_between(x_samples, -1.0, 1.0, alpha=0.1, facecolor='k')
-842 ax1.set_xlim([xstart, xstop])
-843 ax1.set_ylabel('Residuals')
-844 plt.subplots_adjust(wspace=None, hspace=None)
-845 plt.draw()
-846
-847
-848def error_band(x, func, beta):
-849 """Calculate the error band for an array of sample values x, for given fit function func with optimized parameters beta.
-850
-851 Returns
-852 -------
-853 err : np.array(Obs)
-854 Error band for an array of sample values x
-855 """
-856 cov = covariance(beta)
-857 if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps):
-858 warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning)
-859
-860 deriv = []
-861 for i, item in enumerate(x):
-862 deriv.append(np.array(egrad(func)([o.value for o in beta], item)))
-863
-864 err = []
-865 for i, item in enumerate(x):
-866 err.append(np.sqrt(deriv[i] @ cov @ deriv[i]))
-867 err = np.array(err)
-868
-869 return err
+811 plt.figure(figsize=(8, 8 / 1.618))
+812 gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
+813 ax0 = plt.subplot(gs[0])
+814 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')
+815 ax0.plot(x_samples, func([o.value for o in fit_res], x_samples), label='Fit', zorder=10, ls='-', ms=0)
+816 ax0.set_xticklabels([])
+817 ax0.set_xlim([xstart, xstop])
+818 ax0.set_xticklabels([])
+819 ax0.legend(title=title)
+820
+821 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])
+822 ax1 = plt.subplot(gs[1])
+823 ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
+824 ax1.tick_params(direction='out')
+825 ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
+826 ax1.axhline(y=0.0, ls='--', color='k', marker=" ")
+827 ax1.fill_between(x_samples, -1.0, 1.0, alpha=0.1, facecolor='k')
+828 ax1.set_xlim([xstart, xstop])
+829 ax1.set_ylabel('Residuals')
+830 plt.subplots_adjust(wspace=None, hspace=None)
+831 plt.draw()
+832
+833
+834def error_band(x, func, beta):
+835 """Calculate the error band for an array of sample values x, for given fit function func with optimized parameters beta.
+836
+837 Returns
+838 -------
+839 err : np.array(Obs)
+840 Error band for an array of sample values x
+841 """
+842 cov = covariance(beta)
+843 if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps):
+844 warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning)
+845
+846 deriv = []
+847 for i, item in enumerate(x):
+848 deriv.append(np.array(egrad(func)([o.value for o in beta], item)))
+849
+850 err = []
+851 for i, item in enumerate(x):
+852 err.append(np.sqrt(deriv[i] @ cov @ deriv[i]))
+853 err = np.array(err)
+854
+855 return err
+856
+857
+858def ks_test(objects=None):
+859 """Performs a Kolmogorov–Smirnov test for the p-values of all fit object.
+860
+861 Parameters
+862 ----------
+863 objects : list
+864 List of fit results to include in the analysis (optional).
+865
+866 Returns
+867 -------
+868 None
+869 """
870
-871
-872def ks_test(objects=None):
-873 """Performs a Kolmogorov–Smirnov test for the p-values of all fit object.
-874
-875 Parameters
-876 ----------
-877 objects : list
-878 List of fit results to include in the analysis (optional).
-879
-880 Returns
-881 -------
-882 None
-883 """
-884
-885 if objects is None:
-886 obs_list = []
-887 for obj in gc.get_objects():
-888 if isinstance(obj, Fit_result):
-889 obs_list.append(obj)
-890 else:
-891 obs_list = objects
-892
-893 p_values = [o.p_value for o in obs_list]
-894
-895 bins = len(p_values)
-896 x = np.arange(0, 1.001, 0.001)
-897 plt.plot(x, x, 'k', zorder=1)
-898 plt.xlim(0, 1)
-899 plt.ylim(0, 1)
-900 plt.xlabel('p-value')
-901 plt.ylabel('Cumulative probability')
-902 plt.title(str(bins) + ' p-values')
-903
-904 n = np.arange(1, bins + 1) / np.float64(bins)
-905 Xs = np.sort(p_values)
-906 plt.step(Xs, n)
-907 diffs = n - Xs
-908 loc_max_diff = np.argmax(np.abs(diffs))
-909 loc = Xs[loc_max_diff]
-910 plt.annotate('', xy=(loc, loc), xytext=(loc, loc + diffs[loc_max_diff]), arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
-911 plt.draw()
-912
-913 print(scipy.stats.kstest(p_values, 'uniform'))
+871 if objects is None:
+872 obs_list = []
+873 for obj in gc.get_objects():
+874 if isinstance(obj, Fit_result):
+875 obs_list.append(obj)
+876 else:
+877 obs_list = objects
+878
+879 p_values = [o.p_value for o in obs_list]
+880
+881 bins = len(p_values)
+882 x = np.arange(0, 1.001, 0.001)
+883 plt.plot(x, x, 'k', zorder=1)
+884 plt.xlim(0, 1)
+885 plt.ylim(0, 1)
+886 plt.xlabel('p-value')
+887 plt.ylabel('Cumulative probability')
+888 plt.title(str(bins) + ' p-values')
+889
+890 n = np.arange(1, bins + 1) / np.float64(bins)
+891 Xs = np.sort(p_values)
+892 plt.step(Xs, n)
+893 diffs = n - Xs
+894 loc_max_diff = np.argmax(np.abs(diffs))
+895 loc = Xs[loc_max_diff]
+896 plt.annotate('', xy=(loc, loc), xytext=(loc, loc + diffs[loc_max_diff]), arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
+897 plt.draw()
+898
+899 print(scipy.stats.kstest(p_values, 'uniform'))
@@ -1662,35 +1648,35 @@ Parameters and information on the fitted result.
- 751def fit_lin(x, y, **kwargs):
-752 """Performs a linear fit to y = n + m * x and returns two Obs n, m.
+ 737def fit_lin(x, y, **kwargs):
+738 """Performs a linear fit to y = n + m * x and returns two Obs n, m.
+739
+740 Parameters
+741 ----------
+742 x : list
+743 Can either be a list of floats in which case no xerror is assumed, or
+744 a list of Obs, where the dvalues of the Obs are used as xerror for the fit.
+745 y : list
+746 List of Obs, the dvalues of the Obs are used as yerror for the fit.
+747
+748 Returns
+749 -------
+750 fit_parameters : list[Obs]
+751 LIist of fitted observables.
+752 """
753
-754 Parameters
-755 ----------
-756 x : list
-757 Can either be a list of floats in which case no xerror is assumed, or
-758 a list of Obs, where the dvalues of the Obs are used as xerror for the fit.
-759 y : list
-760 List of Obs, the dvalues of the Obs are used as yerror for the fit.
-761
-762 Returns
-763 -------
-764 fit_parameters : list[Obs]
-765 LIist of fitted observables.
-766 """
-767
-768 def f(a, x):
-769 y = a[0] + a[1] * x
-770 return y
-771
-772 if all(isinstance(n, Obs) for n in x):
-773 out = total_least_squares(x, y, f, **kwargs)
-774 return out.fit_parameters
-775 elif all(isinstance(n, float) or isinstance(n, int) for n in x) or isinstance(x, np.ndarray):
-776 out = least_squares(x, y, f, **kwargs)
-777 return out.fit_parameters
-778 else:
-779 raise Exception('Unsupported types for x')
+754 def f(a, x):
+755 y = a[0] + a[1] * x
+756 return y
+757
+758 if all(isinstance(n, Obs) for n in x):
+759 out = total_least_squares(x, y, f, **kwargs)
+760 return out.fit_parameters
+761 elif all(isinstance(n, float) or isinstance(n, int) for n in x) or isinstance(x, np.ndarray):
+762 out = least_squares(x, y, f, **kwargs)
+763 return out.fit_parameters
+764 else:
+765 raise Exception('Unsupported types for x')
782def qqplot(x, o_y, func, p, title=""):
-783 """Generates a quantile-quantile plot of the fit result which can be used to
-784 check if the residuals of the fit are gaussian distributed.
-785
-786 Returns
-787 -------
-788 None
-789 """
-790
-791 residuals = []
-792 for i_x, i_y in zip(x, o_y):
-793 residuals.append((i_y - func(p, i_x)) / i_y.dvalue)
-794 residuals = sorted(residuals)
-795 my_y = [o.value for o in residuals]
-796 probplot = scipy.stats.probplot(my_y)
-797 my_x = probplot[0][0]
-798 plt.figure(figsize=(8, 8 / 1.618))
-799 plt.errorbar(my_x, my_y, fmt='o')
-800 fit_start = my_x[0]
-801 fit_stop = my_x[-1]
-802 samples = np.arange(fit_start, fit_stop, 0.01)
-803 plt.plot(samples, samples, 'k--', zorder=11, label='Standard normal distribution')
-804 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='-')
-805
-806 plt.xlabel('Theoretical quantiles')
-807 plt.ylabel('Ordered Values')
-808 plt.legend(title=title)
-809 plt.draw()
+ 768def qqplot(x, o_y, func, p, title=""):
+769 """Generates a quantile-quantile plot of the fit result which can be used to
+770 check if the residuals of the fit are gaussian distributed.
+771
+772 Returns
+773 -------
+774 None
+775 """
+776
+777 residuals = []
+778 for i_x, i_y in zip(x, o_y):
+779 residuals.append((i_y - func(p, i_x)) / i_y.dvalue)
+780 residuals = sorted(residuals)
+781 my_y = [o.value for o in residuals]
+782 probplot = scipy.stats.probplot(my_y)
+783 my_x = probplot[0][0]
+784 plt.figure(figsize=(8, 8 / 1.618))
+785 plt.errorbar(my_x, my_y, fmt='o')
+786 fit_start = my_x[0]
+787 fit_stop = my_x[-1]
+788 samples = np.arange(fit_start, fit_stop, 0.01)
+789 plt.plot(samples, samples, 'k--', zorder=11, label='Standard normal distribution')
+790 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='-')
+791
+792 plt.xlabel('Theoretical quantiles')
+793 plt.ylabel('Ordered Values')
+794 plt.legend(title=title)
+795 plt.draw()
812def residual_plot(x, y, func, fit_res, title=""):
-813 """Generates a plot which compares the fit to the data and displays the corresponding residuals
-814
-815 For uncorrelated data the residuals are expected to be distributed ~N(0,1).
-816
-817 Returns
-818 -------
-819 None
-820 """
-821 sorted_x = sorted(x)
-822 xstart = sorted_x[0] - 0.5 * (sorted_x[1] - sorted_x[0])
-823 xstop = sorted_x[-1] + 0.5 * (sorted_x[-1] - sorted_x[-2])
-824 x_samples = np.arange(xstart, xstop + 0.01, 0.01)
-825
-826 plt.figure(figsize=(8, 8 / 1.618))
-827 gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
-828 ax0 = plt.subplot(gs[0])
-829 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')
-830 ax0.plot(x_samples, func([o.value for o in fit_res], x_samples), label='Fit', zorder=10, ls='-', ms=0)
-831 ax0.set_xticklabels([])
-832 ax0.set_xlim([xstart, xstop])
-833 ax0.set_xticklabels([])
-834 ax0.legend(title=title)
-835
-836 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])
-837 ax1 = plt.subplot(gs[1])
-838 ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
-839 ax1.tick_params(direction='out')
-840 ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
-841 ax1.axhline(y=0.0, ls='--', color='k', marker=" ")
-842 ax1.fill_between(x_samples, -1.0, 1.0, alpha=0.1, facecolor='k')
-843 ax1.set_xlim([xstart, xstop])
-844 ax1.set_ylabel('Residuals')
-845 plt.subplots_adjust(wspace=None, hspace=None)
-846 plt.draw()
+ 798def residual_plot(x, y, func, fit_res, title=""):
+799 """Generates a plot which compares the fit to the data and displays the corresponding residuals
+800
+801 For uncorrelated data the residuals are expected to be distributed ~N(0,1).
+802
+803 Returns
+804 -------
+805 None
+806 """
+807 sorted_x = sorted(x)
+808 xstart = sorted_x[0] - 0.5 * (sorted_x[1] - sorted_x[0])
+809 xstop = sorted_x[-1] + 0.5 * (sorted_x[-1] - sorted_x[-2])
+810 x_samples = np.arange(xstart, xstop + 0.01, 0.01)
+811
+812 plt.figure(figsize=(8, 8 / 1.618))
+813 gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
+814 ax0 = plt.subplot(gs[0])
+815 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')
+816 ax0.plot(x_samples, func([o.value for o in fit_res], x_samples), label='Fit', zorder=10, ls='-', ms=0)
+817 ax0.set_xticklabels([])
+818 ax0.set_xlim([xstart, xstop])
+819 ax0.set_xticklabels([])
+820 ax0.legend(title=title)
+821
+822 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])
+823 ax1 = plt.subplot(gs[1])
+824 ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
+825 ax1.tick_params(direction='out')
+826 ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
+827 ax1.axhline(y=0.0, ls='--', color='k', marker=" ")
+828 ax1.fill_between(x_samples, -1.0, 1.0, alpha=0.1, facecolor='k')
+829 ax1.set_xlim([xstart, xstop])
+830 ax1.set_ylabel('Residuals')
+831 plt.subplots_adjust(wspace=None, hspace=None)
+832 plt.draw()
849def error_band(x, func, beta):
-850 """Calculate the error band for an array of sample values x, for given fit function func with optimized parameters beta.
-851
-852 Returns
-853 -------
-854 err : np.array(Obs)
-855 Error band for an array of sample values x
-856 """
-857 cov = covariance(beta)
-858 if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps):
-859 warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning)
-860
-861 deriv = []
-862 for i, item in enumerate(x):
-863 deriv.append(np.array(egrad(func)([o.value for o in beta], item)))
-864
-865 err = []
-866 for i, item in enumerate(x):
-867 err.append(np.sqrt(deriv[i] @ cov @ deriv[i]))
-868 err = np.array(err)
-869
-870 return err
+ 835def error_band(x, func, beta):
+836 """Calculate the error band for an array of sample values x, for given fit function func with optimized parameters beta.
+837
+838 Returns
+839 -------
+840 err : np.array(Obs)
+841 Error band for an array of sample values x
+842 """
+843 cov = covariance(beta)
+844 if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps):
+845 warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning)
+846
+847 deriv = []
+848 for i, item in enumerate(x):
+849 deriv.append(np.array(egrad(func)([o.value for o in beta], item)))
+850
+851 err = []
+852 for i, item in enumerate(x):
+853 err.append(np.sqrt(deriv[i] @ cov @ deriv[i]))
+854 err = np.array(err)
+855
+856 return err
873def ks_test(objects=None):
-874 """Performs a Kolmogorov–Smirnov test for the p-values of all fit object.
-875
-876 Parameters
-877 ----------
-878 objects : list
-879 List of fit results to include in the analysis (optional).
-880
-881 Returns
-882 -------
-883 None
-884 """
-885
-886 if objects is None:
-887 obs_list = []
-888 for obj in gc.get_objects():
-889 if isinstance(obj, Fit_result):
-890 obs_list.append(obj)
-891 else:
-892 obs_list = objects
-893
-894 p_values = [o.p_value for o in obs_list]
-895
-896 bins = len(p_values)
-897 x = np.arange(0, 1.001, 0.001)
-898 plt.plot(x, x, 'k', zorder=1)
-899 plt.xlim(0, 1)
-900 plt.ylim(0, 1)
-901 plt.xlabel('p-value')
-902 plt.ylabel('Cumulative probability')
-903 plt.title(str(bins) + ' p-values')
-904
-905 n = np.arange(1, bins + 1) / np.float64(bins)
-906 Xs = np.sort(p_values)
-907 plt.step(Xs, n)
-908 diffs = n - Xs
-909 loc_max_diff = np.argmax(np.abs(diffs))
-910 loc = Xs[loc_max_diff]
-911 plt.annotate('', xy=(loc, loc), xytext=(loc, loc + diffs[loc_max_diff]), arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
-912 plt.draw()
-913
-914 print(scipy.stats.kstest(p_values, 'uniform'))
+ 859def ks_test(objects=None):
+860 """Performs a Kolmogorov–Smirnov test for the p-values of all fit object.
+861
+862 Parameters
+863 ----------
+864 objects : list
+865 List of fit results to include in the analysis (optional).
+866
+867 Returns
+868 -------
+869 None
+870 """
+871
+872 if objects is None:
+873 obs_list = []
+874 for obj in gc.get_objects():
+875 if isinstance(obj, Fit_result):
+876 obs_list.append(obj)
+877 else:
+878 obs_list = objects
+879
+880 p_values = [o.p_value for o in obs_list]
+881
+882 bins = len(p_values)
+883 x = np.arange(0, 1.001, 0.001)
+884 plt.plot(x, x, 'k', zorder=1)
+885 plt.xlim(0, 1)
+886 plt.ylim(0, 1)
+887 plt.xlabel('p-value')
+888 plt.ylabel('Cumulative probability')
+889 plt.title(str(bins) + ' p-values')
+890
+891 n = np.arange(1, bins + 1) / np.float64(bins)
+892 Xs = np.sort(p_values)
+893 plt.step(Xs, n)
+894 diffs = n - Xs
+895 loc_max_diff = np.argmax(np.abs(diffs))
+896 loc = Xs[loc_max_diff]
+897 plt.annotate('', xy=(loc, loc), xytext=(loc, loc + diffs[loc_max_diff]), arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
+898 plt.draw()
+899
+900 print(scipy.stats.kstest(p_values, 'uniform'))