From 09e20b4183d10a3a3c7cd42ba50eaa7bcfe77479 Mon Sep 17 00:00:00 2001 From: fjosw Date: Fri, 3 Mar 2023 16:35:37 +0000 Subject: [PATCH] Documentation updated --- docs/pyerrors/fits.html | 936 ++++++++++++++++++++-------------------- 1 file changed, 461 insertions(+), 475 deletions(-) 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')
 

 
 
@@ -1727,34 +1713,34 @@ LIist of fitted observables.
- 
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()
 

 
 
@@ -1781,41 +1767,41 @@ LIist of fitted observables.
- 
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()
 

 
 
@@ -1843,28 +1829,28 @@ LIist of fitted observables.
- 
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
 

 
 
@@ -1891,48 +1877,48 @@ Error band for an array of sample values x
- 
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'))