From 0de101c47103ce880e489fbb2f30677f81833806 Mon Sep 17 00:00:00 2001 From: fjosw Date: Fri, 27 May 2022 12:51:52 +0000 Subject: [PATCH] Documentation updated --- docs/pyerrors/fits.html | 775 ++++++++++++++++++++-------------------- 1 file changed, 387 insertions(+), 388 deletions(-) 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')
 

 
 
@@ -1407,30 +1406,30 @@ List of Obs, the dvalues of the Obs are used as yerror for the fit.
- 
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()
 

 
 
@@ -1451,34 +1450,34 @@ check if the residuals of the fit are gaussian distributed.
- 
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()
 

 
 
@@ -1498,22 +1497,22 @@ check if the residuals of the fit are gaussian distributed.
- 
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
 

 
 
@@ -1533,44 +1532,44 @@ check if the residuals of the fit are gaussian distributed.
- 
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'))