diff --git a/docs/pyerrors/fits.html b/docs/pyerrors/fits.html
index e9ebafa7..aabb6154 100644
--- a/docs/pyerrors/fits.html
+++ b/docs/pyerrors/fits.html
@@ -816,35 +816,35 @@
705
706def _combined_fit(x, y, func, silent=False, **kwargs):
707
- 708 if kwargs.get('correlated_fit') is True:
- 709 raise Exception("Correlated fit has not been implemented yet")
+ 708 output = Fit_result()
+ 709 output.fit_function = func
710
- 711 output = Fit_result()
- 712 output.fit_function = func
- 713
- 714 if kwargs.get('num_grad') is True:
- 715 jacobian = num_jacobian
- 716 hessian = num_hessian
- 717 else:
- 718 jacobian = auto_jacobian
- 719 hessian = auto_hessian
- 720
- 721 key_ls = sorted(list(x.keys()))
+ 711 if kwargs.get('num_grad') is True:
+ 712 jacobian = num_jacobian
+ 713 hessian = num_hessian
+ 714 else:
+ 715 jacobian = auto_jacobian
+ 716 hessian = auto_hessian
+ 717
+ 718 key_ls = sorted(list(x.keys()))
+ 719
+ 720 if sorted(list(y.keys())) != key_ls:
+ 721 raise Exception('x and y dictionaries do not contain the same keys.')
722
- 723 if sorted(list(y.keys())) != key_ls:
- 724 raise Exception('x and y dictionaries do not contain the same keys.')
+ 723 if sorted(list(func.keys())) != key_ls:
+ 724 raise Exception('x and func dictionaries do not contain the same keys.')
725
- 726 if sorted(list(func.keys())) != key_ls:
- 727 raise Exception('x and func dictionaries do not contain the same keys.')
+ 726 x_all = np.concatenate([np.array(x[key]) for key in key_ls])
+ 727 y_all = np.concatenate([np.array(y[key]) for key in key_ls])
728
- 729 x_all = np.concatenate([np.array(x[key]) for key in key_ls])
- 730 y_all = np.concatenate([np.array(y[key]) for key in key_ls])
+ 729 y_f = [o.value for o in y_all]
+ 730 dy_f = [o.dvalue for o in y_all]
731
- 732 y_f = [o.value for o in y_all]
- 733 dy_f = [o.dvalue for o in y_all]
+ 732 if len(x_all.shape) > 2:
+ 733 raise Exception('Unknown format for x values')
734
- 735 if len(x_all.shape) > 2:
- 736 raise Exception('Unknown format for x values')
+ 735 if np.any(np.asarray(dy_f) <= 0.0):
+ 736 raise Exception('No y errors available, run the gamma method first.')
737
738 # number of fit parameters
739 n_parms_ls = []
@@ -877,279 +877,326 @@
766 else:
767 x0 = [0.1] * n_parms
768
- 769 def chisqfunc(p):
- 770 func_list = np.concatenate([[func[k]] * len(x[k]) for k in key_ls])
- 771 model = anp.array([func_list[i](p, x_all[i]) for i in range(len(x_all))])
- 772 chisq = anp.sum(((y_f - model) / dy_f) ** 2)
- 773 return chisq
- 774
- 775 output.method = kwargs.get('method', 'Levenberg-Marquardt')
- 776 if not silent:
- 777 print('Method:', output.method)
- 778
- 779 if output.method != 'Levenberg-Marquardt':
- 780 if output.method == 'migrad':
- 781 tolerance = 1e-4 # default value of 1e-1 set by iminuit can be problematic
- 782 if 'tol' in kwargs:
- 783 tolerance = kwargs.get('tol')
- 784 fit_result = iminuit.minimize(chisqfunc, x0, tol=tolerance) # Stopping criterion 0.002 * tol * errordef
- 785 output.iterations = fit_result.nfev
- 786 else:
- 787 tolerance = 1e-12
- 788 if 'tol' in kwargs:
- 789 tolerance = kwargs.get('tol')
- 790 fit_result = scipy.optimize.minimize(chisqfunc, x0, method=kwargs.get('method'), tol=tolerance)
- 791 output.iterations = fit_result.nit
- 792
- 793 chisquare = fit_result.fun
+ 769 if kwargs.get('correlated_fit') is True:
+ 770 corr = covariance(y_all, correlation=True, **kwargs)
+ 771 covdiag = np.diag(1 / np.asarray(dy_f))
+ 772 condn = np.linalg.cond(corr)
+ 773 if condn > 0.1 / np.finfo(float).eps:
+ 774 raise Exception(f"Cannot invert correlation matrix as its condition number exceeds machine precision ({condn:1.2e})")
+ 775 if condn > 1e13:
+ 776 warnings.warn("Correlation matrix may be ill-conditioned, condition number: {%1.2e}" % (condn), RuntimeWarning)
+ 777 chol = np.linalg.cholesky(corr)
+ 778 chol_inv = scipy.linalg.solve_triangular(chol, covdiag, lower=True)
+ 779
+ 780 def chisqfunc_corr(p):
+ 781 model = np.concatenate([np.array(func[key](p, np.asarray(x[key]))) for key in key_ls])
+ 782 chisq = anp.sum(anp.dot(chol_inv, (y_f - model)) ** 2)
+ 783 return chisq
+ 784
+ 785 def chisqfunc(p):
+ 786 func_list = np.concatenate([[func[k]] * len(x[k]) for k in key_ls])
+ 787 model = anp.array([func_list[i](p, x_all[i]) for i in range(len(x_all))])
+ 788 chisq = anp.sum(((y_f - model) / dy_f) ** 2)
+ 789 return chisq
+ 790
+ 791 output.method = kwargs.get('method', 'Levenberg-Marquardt')
+ 792 if not silent:
+ 793 print('Method:', output.method)
794
- 795 else:
- 796 def chisqfunc_residuals(p):
- 797 model = np.concatenate([np.array(func[key](p, np.asarray(x[key]))) for key in key_ls])
- 798 chisq = ((y_f - model) / dy_f)
- 799 return chisq
- 800 if 'tol' in kwargs:
- 801 print('tol cannot be set for Levenberg-Marquardt')
- 802 fit_result = scipy.optimize.least_squares(chisqfunc_residuals, x0, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15)
- 803
- 804 chisquare = np.sum(fit_result.fun ** 2)
- 805 assert np.isclose(chisquare, chisqfunc(fit_result.x), atol=1e-14)
- 806 output.iterations = fit_result.nfev
- 807
- 808 output.message = fit_result.message
- 809
- 810 if not fit_result.success:
- 811 raise Exception('The minimization procedure did not converge.')
+ 795 if output.method != 'Levenberg-Marquardt':
+ 796 if output.method == 'migrad':
+ 797 tolerance = 1e-4 # default value of 1e-1 set by iminuit can be problematic
+ 798 if 'tol' in kwargs:
+ 799 tolerance = kwargs.get('tol')
+ 800 fit_result = iminuit.minimize(chisqfunc, x0, tol=tolerance) # Stopping criterion 0.002 * tol * errordef
+ 801 if kwargs.get('correlated_fit') is True:
+ 802 fit_result = iminuit.minimize(chisqfunc_corr, fit_result.x, tol=tolerance)
+ 803 output.iterations = fit_result.nfev
+ 804 else:
+ 805 tolerance = 1e-12
+ 806 if 'tol' in kwargs:
+ 807 tolerance = kwargs.get('tol')
+ 808 fit_result = scipy.optimize.minimize(chisqfunc, x0, method=kwargs.get('method'), tol=tolerance)
+ 809 if kwargs.get('correlated_fit') is True:
+ 810 fit_result = scipy.optimize.minimize(chisqfunc_corr, fit_result.x, method=kwargs.get('method'), tol=tolerance)
+ 811 output.iterations = fit_result.nit
812
- 813 if x_all.shape[-1] - n_parms > 0:
- 814 output.chisquare = chisquare
- 815 output.dof = x_all.shape[-1] - n_parms
- 816 output.chisquare_by_dof = output.chisquare / output.dof
- 817 output.p_value = 1 - scipy.stats.chi2.cdf(output.chisquare, output.dof)
- 818 else:
- 819 output.chisquare_by_dof = float('nan')
- 820
- 821 if not silent:
- 822 print(fit_result.message)
- 823 print('chisquare/d.o.f.:', output.chisquare_by_dof)
- 824 print('fit parameters', fit_result.x)
- 825
- 826 def chisqfunc_compact(d):
- 827 func_list = np.concatenate([[func[k]] * len(x[k]) for k in key_ls])
- 828 model = anp.array([func_list[i](d[:n_parms], x_all[i]) for i in range(len(x_all))])
- 829 chisq = anp.sum(((d[n_parms:] - model) / dy_f) ** 2)
- 830 return chisq
- 831
- 832 def prepare_hat_matrix():
- 833 hat_vector = []
- 834 for key in key_ls:
- 835 x_array = np.asarray(x[key])
- 836 if (len(x_array) != 0):
- 837 hat_vector.append(jacobian(func[key])(fit_result.x, x_array))
- 838 hat_vector = [item for sublist in hat_vector for item in sublist]
- 839 return hat_vector
- 840
- 841 fitp = fit_result.x
- 842
- 843 if np.any(np.asarray(dy_f) <= 0.0):
- 844 raise Exception('No y errors available, run the gamma method first.')
- 845
- 846 try:
- 847 hess = hessian(chisqfunc)(fitp)
- 848 except TypeError:
- 849 raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None
- 850
- 851 jac_jac_y = hessian(chisqfunc_compact)(np.concatenate((fitp, y_f)))
+ 813 chisquare = fit_result.fun
+ 814
+ 815 else:
+ 816 if kwargs.get('correlated_fit') is True:
+ 817 def chisqfunc_residuals_corr(p):
+ 818 model = np.concatenate([np.array(func[key](p, np.asarray(x[key]))) for key in key_ls])
+ 819 chisq = anp.dot(chol_inv, (y_f - model))
+ 820 return chisq
+ 821
+ 822 def chisqfunc_residuals(p):
+ 823 model = np.concatenate([np.array(func[key](p, np.asarray(x[key]))) for key in key_ls])
+ 824 chisq = ((y_f - model) / dy_f)
+ 825 return chisq
+ 826
+ 827 if 'tol' in kwargs:
+ 828 print('tol cannot be set for Levenberg-Marquardt')
+ 829
+ 830 fit_result = scipy.optimize.least_squares(chisqfunc_residuals, x0, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15)
+ 831 if kwargs.get('correlated_fit') is True:
+ 832 fit_result = scipy.optimize.least_squares(chisqfunc_residuals_corr, fit_result.x, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15)
+ 833
+ 834 chisquare = np.sum(fit_result.fun ** 2)
+ 835 if kwargs.get('correlated_fit') is True:
+ 836 assert np.isclose(chisquare, chisqfunc_corr(fit_result.x), atol=1e-14)
+ 837 else:
+ 838 assert np.isclose(chisquare, chisqfunc(fit_result.x), atol=1e-14)
+ 839
+ 840 output.iterations = fit_result.nfev
+ 841
+ 842 if not fit_result.success:
+ 843 raise Exception('The minimization procedure did not converge.')
+ 844
+ 845 if x_all.shape[-1] - n_parms > 0:
+ 846 output.chisquare = chisquare
+ 847 output.dof = x_all.shape[-1] - n_parms
+ 848 output.chisquare_by_dof = output.chisquare / output.dof
+ 849 output.p_value = 1 - scipy.stats.chi2.cdf(output.chisquare, output.dof)
+ 850 else:
+ 851 output.chisquare_by_dof = float('nan')
852
- 853 # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv
- 854 try:
- 855 deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms, n_parms:])
- 856 except np.linalg.LinAlgError:
- 857 raise Exception("Cannot invert hessian matrix.")
+ 853 output.message = fit_result.message
+ 854 if not silent:
+ 855 print(fit_result.message)
+ 856 print('chisquare/d.o.f.:', output.chisquare_by_dof)
+ 857 print('fit parameters', fit_result.x)
858
- 859 if kwargs.get('expected_chisquare') is True:
- 860 if kwargs.get('correlated_fit') is not True:
- 861 W = np.diag(1 / np.asarray(dy_f))
- 862 cov = covariance(y_all)
- 863 hat_vector = prepare_hat_matrix()
- 864 A = W @ hat_vector # hat_vector = 'jacobian(func)(fit_result.x, x)'
- 865 P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T
- 866 expected_chisquare = np.trace((np.identity(x_all.shape[-1]) - P_phi) @ W @ cov @ W)
- 867 output.chisquare_by_expected_chisquare = output.chisquare / expected_chisquare
- 868 if not silent:
- 869 print('chisquare/expected_chisquare:', output.chisquare_by_expected_chisquare)
- 870
- 871 result = []
- 872 for i in range(n_parms):
- 873 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])))
- 874
- 875 output.fit_parameters = result
- 876
- 877 return output
- 878
+ 859 def prepare_hat_matrix():
+ 860 hat_vector = []
+ 861 for key in key_ls:
+ 862 x_array = np.asarray(x[key])
+ 863 if (len(x_array) != 0):
+ 864 hat_vector.append(jacobian(func[key])(fit_result.x, x_array))
+ 865 hat_vector = [item for sublist in hat_vector for item in sublist]
+ 866 return hat_vector
+ 867
+ 868 if kwargs.get('expected_chisquare') is True:
+ 869 if kwargs.get('correlated_fit') is not True:
+ 870 W = np.diag(1 / np.asarray(dy_f))
+ 871 cov = covariance(y_all)
+ 872 hat_vector = prepare_hat_matrix()
+ 873 A = W @ hat_vector # hat_vector = 'jacobian(func)(fit_result.x, x)'
+ 874 P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T
+ 875 expected_chisquare = np.trace((np.identity(x_all.shape[-1]) - P_phi) @ W @ cov @ W)
+ 876 output.chisquare_by_expected_chisquare = output.chisquare / expected_chisquare
+ 877 if not silent:
+ 878 print('chisquare/expected_chisquare:', output.chisquare_by_expected_chisquare)
879
- 880def fit_lin(x, y, **kwargs):
- 881 """Performs a linear fit to y = n + m * x and returns two Obs n, m.
- 882
- 883 Parameters
- 884 ----------
- 885 x : list
- 886 Can either be a list of floats in which case no xerror is assumed, or
- 887 a list of Obs, where the dvalues of the Obs are used as xerror for the fit.
- 888 y : list
- 889 List of Obs, the dvalues of the Obs are used as yerror for the fit.
- 890
- 891 Returns
- 892 -------
- 893 fit_parameters : list[Obs]
- 894 LIist of fitted observables.
- 895 """
- 896
- 897 def f(a, x):
- 898 y = a[0] + a[1] * x
- 899 return y
- 900
- 901 if all(isinstance(n, Obs) for n in x):
- 902 out = total_least_squares(x, y, f, **kwargs)
- 903 return out.fit_parameters
- 904 elif all(isinstance(n, float) or isinstance(n, int) for n in x) or isinstance(x, np.ndarray):
- 905 out = least_squares(x, y, f, **kwargs)
- 906 return out.fit_parameters
- 907 else:
- 908 raise Exception('Unsupported types for x')
- 909
- 910
- 911def qqplot(x, o_y, func, p):
- 912 """Generates a quantile-quantile plot of the fit result which can be used to
- 913 check if the residuals of the fit are gaussian distributed.
- 914
- 915 Returns
- 916 -------
- 917 None
- 918 """
- 919
- 920 residuals = []
- 921 for i_x, i_y in zip(x, o_y):
- 922 residuals.append((i_y - func(p, i_x)) / i_y.dvalue)
- 923 residuals = sorted(residuals)
- 924 my_y = [o.value for o in residuals]
- 925 probplot = scipy.stats.probplot(my_y)
- 926 my_x = probplot[0][0]
- 927 plt.figure(figsize=(8, 8 / 1.618))
- 928 plt.errorbar(my_x, my_y, fmt='o')
- 929 fit_start = my_x[0]
- 930 fit_stop = my_x[-1]
- 931 samples = np.arange(fit_start, fit_stop, 0.01)
- 932 plt.plot(samples, samples, 'k--', zorder=11, label='Standard normal distribution')
- 933 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='-')
- 934
- 935 plt.xlabel('Theoretical quantiles')
- 936 plt.ylabel('Ordered Values')
- 937 plt.legend()
- 938 plt.draw()
- 939
- 940
- 941def residual_plot(x, y, func, fit_res):
- 942 """Generates a plot which compares the fit to the data and displays the corresponding residuals
+ 880 fitp = fit_result.x
+ 881 if np.any(np.asarray(dy_f) <= 0.0):
+ 882 raise Exception('No y errors available, run the gamma method first.')
+ 883
+ 884 try:
+ 885 if kwargs.get('correlated_fit') is True:
+ 886 hess = hessian(chisqfunc_corr)(fitp)
+ 887 else:
+ 888 hess = hessian(chisqfunc)(fitp)
+ 889 except TypeError:
+ 890 raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None
+ 891
+ 892 if kwargs.get('correlated_fit') is True:
+ 893 def chisqfunc_compact(d):
+ 894 func_list = np.concatenate([[func[k]] * len(x[k]) for k in key_ls])
+ 895 model = anp.array([func_list[i](d[:n_parms], x_all[i]) for i in range(len(x_all))])
+ 896 chisq = anp.sum(anp.dot(chol_inv, (d[n_parms:] - model)) ** 2)
+ 897 return chisq
+ 898 else:
+ 899 def chisqfunc_compact(d):
+ 900 func_list = np.concatenate([[func[k]] * len(x[k]) for k in key_ls])
+ 901 model = anp.array([func_list[i](d[:n_parms], x_all[i]) for i in range(len(x_all))])
+ 902 chisq = anp.sum(((d[n_parms:] - model) / dy_f) ** 2)
+ 903 return chisq
+ 904
+ 905 jac_jac_y = hessian(chisqfunc_compact)(np.concatenate((fitp, y_f)))
+ 906
+ 907 # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv
+ 908 try:
+ 909 deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms, n_parms:])
+ 910 except np.linalg.LinAlgError:
+ 911 raise Exception("Cannot invert hessian matrix.")
+ 912
+ 913 result = []
+ 914 for i in range(n_parms):
+ 915 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])))
+ 916
+ 917 output.fit_parameters = result
+ 918
+ 919 if kwargs.get('correlated_fit') is True:
+ 920 n_cov = np.min(np.vectorize(lambda x_all: x_all.N)(y_all))
+ 921 output.t2_p_value = 1 - scipy.stats.f.cdf((n_cov - output.dof) / (output.dof * (n_cov - 1)) * output.chisquare,
+ 922 output.dof, n_cov - output.dof)
+ 923
+ 924 return output
+ 925
+ 926
+ 927def fit_lin(x, y, **kwargs):
+ 928 """Performs a linear fit to y = n + m * x and returns two Obs n, m.
+ 929
+ 930 Parameters
+ 931 ----------
+ 932 x : list
+ 933 Can either be a list of floats in which case no xerror is assumed, or
+ 934 a list of Obs, where the dvalues of the Obs are used as xerror for the fit.
+ 935 y : list
+ 936 List of Obs, the dvalues of the Obs are used as yerror for the fit.
+ 937
+ 938 Returns
+ 939 -------
+ 940 fit_parameters : list[Obs]
+ 941 LIist of fitted observables.
+ 942 """
943
- 944 Returns
- 945 -------
- 946 None
- 947 """
- 948 sorted_x = sorted(x)
- 949 xstart = sorted_x[0] - 0.5 * (sorted_x[1] - sorted_x[0])
- 950 xstop = sorted_x[-1] + 0.5 * (sorted_x[-1] - sorted_x[-2])
- 951 x_samples = np.arange(xstart, xstop + 0.01, 0.01)
- 952
- 953 plt.figure(figsize=(8, 8 / 1.618))
- 954 gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
- 955 ax0 = plt.subplot(gs[0])
- 956 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')
- 957 ax0.plot(x_samples, func([o.value for o in fit_res], x_samples), label='Fit', zorder=10, ls='-', ms=0)
- 958 ax0.set_xticklabels([])
- 959 ax0.set_xlim([xstart, xstop])
- 960 ax0.set_xticklabels([])
- 961 ax0.legend()
- 962
- 963 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])
- 964 ax1 = plt.subplot(gs[1])
- 965 ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
- 966 ax1.tick_params(direction='out')
- 967 ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
- 968 ax1.axhline(y=0.0, ls='--', color='k', marker=" ")
- 969 ax1.fill_between(x_samples, -1.0, 1.0, alpha=0.1, facecolor='k')
- 970 ax1.set_xlim([xstart, xstop])
- 971 ax1.set_ylabel('Residuals')
- 972 plt.subplots_adjust(wspace=None, hspace=None)
- 973 plt.draw()
- 974
- 975
- 976def error_band(x, func, beta):
- 977 """Calculate the error band for an array of sample values x, for given fit function func with optimized parameters beta.
- 978
- 979 Returns
- 980 -------
- 981 err : np.array(Obs)
- 982 Error band for an array of sample values x
- 983 """
- 984 cov = covariance(beta)
- 985 if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps):
- 986 warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning)
+ 944 def f(a, x):
+ 945 y = a[0] + a[1] * x
+ 946 return y
+ 947
+ 948 if all(isinstance(n, Obs) for n in x):
+ 949 out = total_least_squares(x, y, f, **kwargs)
+ 950 return out.fit_parameters
+ 951 elif all(isinstance(n, float) or isinstance(n, int) for n in x) or isinstance(x, np.ndarray):
+ 952 out = least_squares(x, y, f, **kwargs)
+ 953 return out.fit_parameters
+ 954 else:
+ 955 raise Exception('Unsupported types for x')
+ 956
+ 957
+ 958def qqplot(x, o_y, func, p):
+ 959 """Generates a quantile-quantile plot of the fit result which can be used to
+ 960 check if the residuals of the fit are gaussian distributed.
+ 961
+ 962 Returns
+ 963 -------
+ 964 None
+ 965 """
+ 966
+ 967 residuals = []
+ 968 for i_x, i_y in zip(x, o_y):
+ 969 residuals.append((i_y - func(p, i_x)) / i_y.dvalue)
+ 970 residuals = sorted(residuals)
+ 971 my_y = [o.value for o in residuals]
+ 972 probplot = scipy.stats.probplot(my_y)
+ 973 my_x = probplot[0][0]
+ 974 plt.figure(figsize=(8, 8 / 1.618))
+ 975 plt.errorbar(my_x, my_y, fmt='o')
+ 976 fit_start = my_x[0]
+ 977 fit_stop = my_x[-1]
+ 978 samples = np.arange(fit_start, fit_stop, 0.01)
+ 979 plt.plot(samples, samples, 'k--', zorder=11, label='Standard normal distribution')
+ 980 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='-')
+ 981
+ 982 plt.xlabel('Theoretical quantiles')
+ 983 plt.ylabel('Ordered Values')
+ 984 plt.legend()
+ 985 plt.draw()
+ 986
987
- 988 deriv = []
- 989 for i, item in enumerate(x):
- 990 deriv.append(np.array(egrad(func)([o.value for o in beta], item)))
- 991
- 992 err = []
- 993 for i, item in enumerate(x):
- 994 err.append(np.sqrt(deriv[i] @ cov @ deriv[i]))
- 995 err = np.array(err)
- 996
- 997 return err
- 998
+ 988def residual_plot(x, y, func, fit_res):
+ 989 """Generates a plot which compares the fit to the data and displays the corresponding residuals
+ 990
+ 991 Returns
+ 992 -------
+ 993 None
+ 994 """
+ 995 sorted_x = sorted(x)
+ 996 xstart = sorted_x[0] - 0.5 * (sorted_x[1] - sorted_x[0])
+ 997 xstop = sorted_x[-1] + 0.5 * (sorted_x[-1] - sorted_x[-2])
+ 998 x_samples = np.arange(xstart, xstop + 0.01, 0.01)
999
-1000def ks_test(objects=None):
-1001 """Performs a Kolmogorov–Smirnov test for the p-values of all fit object.
-1002
-1003 Parameters
-1004 ----------
-1005 objects : list
-1006 List of fit results to include in the analysis (optional).
-1007
-1008 Returns
-1009 -------
-1010 None
-1011 """
-1012
-1013 if objects is None:
-1014 obs_list = []
-1015 for obj in gc.get_objects():
-1016 if isinstance(obj, Fit_result):
-1017 obs_list.append(obj)
-1018 else:
-1019 obs_list = objects
-1020
-1021 p_values = [o.p_value for o in obs_list]
+1000 plt.figure(figsize=(8, 8 / 1.618))
+1001 gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
+1002 ax0 = plt.subplot(gs[0])
+1003 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')
+1004 ax0.plot(x_samples, func([o.value for o in fit_res], x_samples), label='Fit', zorder=10, ls='-', ms=0)
+1005 ax0.set_xticklabels([])
+1006 ax0.set_xlim([xstart, xstop])
+1007 ax0.set_xticklabels([])
+1008 ax0.legend()
+1009
+1010 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])
+1011 ax1 = plt.subplot(gs[1])
+1012 ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
+1013 ax1.tick_params(direction='out')
+1014 ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
+1015 ax1.axhline(y=0.0, ls='--', color='k', marker=" ")
+1016 ax1.fill_between(x_samples, -1.0, 1.0, alpha=0.1, facecolor='k')
+1017 ax1.set_xlim([xstart, xstop])
+1018 ax1.set_ylabel('Residuals')
+1019 plt.subplots_adjust(wspace=None, hspace=None)
+1020 plt.draw()
+1021
1022
-1023 bins = len(p_values)
-1024 x = np.arange(0, 1.001, 0.001)
-1025 plt.plot(x, x, 'k', zorder=1)
-1026 plt.xlim(0, 1)
-1027 plt.ylim(0, 1)
-1028 plt.xlabel('p-value')
-1029 plt.ylabel('Cumulative probability')
-1030 plt.title(str(bins) + ' p-values')
-1031
-1032 n = np.arange(1, bins + 1) / np.float64(bins)
-1033 Xs = np.sort(p_values)
-1034 plt.step(Xs, n)
-1035 diffs = n - Xs
-1036 loc_max_diff = np.argmax(np.abs(diffs))
-1037 loc = Xs[loc_max_diff]
-1038 plt.annotate('', xy=(loc, loc), xytext=(loc, loc + diffs[loc_max_diff]), arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
-1039 plt.draw()
-1040
-1041 print(scipy.stats.kstest(p_values, 'uniform'))
+1023def error_band(x, func, beta):
+1024 """Calculate the error band for an array of sample values x, for given fit function func with optimized parameters beta.
+1025
+1026 Returns
+1027 -------
+1028 err : np.array(Obs)
+1029 Error band for an array of sample values x
+1030 """
+1031 cov = covariance(beta)
+1032 if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps):
+1033 warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning)
+1034
+1035 deriv = []
+1036 for i, item in enumerate(x):
+1037 deriv.append(np.array(egrad(func)([o.value for o in beta], item)))
+1038
+1039 err = []
+1040 for i, item in enumerate(x):
+1041 err.append(np.sqrt(deriv[i] @ cov @ deriv[i]))
+1042 err = np.array(err)
+1043
+1044 return err
+1045
+1046
+1047def ks_test(objects=None):
+1048 """Performs a Kolmogorov–Smirnov test for the p-values of all fit object.
+1049
+1050 Parameters
+1051 ----------
+1052 objects : list
+1053 List of fit results to include in the analysis (optional).
+1054
+1055 Returns
+1056 -------
+1057 None
+1058 """
+1059
+1060 if objects is None:
+1061 obs_list = []
+1062 for obj in gc.get_objects():
+1063 if isinstance(obj, Fit_result):
+1064 obs_list.append(obj)
+1065 else:
+1066 obs_list = objects
+1067
+1068 p_values = [o.p_value for o in obs_list]
+1069
+1070 bins = len(p_values)
+1071 x = np.arange(0, 1.001, 0.001)
+1072 plt.plot(x, x, 'k', zorder=1)
+1073 plt.xlim(0, 1)
+1074 plt.ylim(0, 1)
+1075 plt.xlabel('p-value')
+1076 plt.ylabel('Cumulative probability')
+1077 plt.title(str(bins) + ' p-values')
+1078
+1079 n = np.arange(1, bins + 1) / np.float64(bins)
+1080 Xs = np.sort(p_values)
+1081 plt.step(Xs, n)
+1082 diffs = n - Xs
+1083 loc_max_diff = np.argmax(np.abs(diffs))
+1084 loc = Xs[loc_max_diff]
+1085 plt.annotate('', xy=(loc, loc), xytext=(loc, loc + diffs[loc_max_diff]), arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
+1086 plt.draw()
+1087
+1088 print(scipy.stats.kstest(p_values, 'uniform'))
@@ -1816,35 +1863,35 @@ Parameters and information on the fitted result.
- 881def fit_lin(x, y, **kwargs):
-882 """Performs a linear fit to y = n + m * x and returns two Obs n, m.
-883
-884 Parameters
-885 ----------
-886 x : list
-887 Can either be a list of floats in which case no xerror is assumed, or
-888 a list of Obs, where the dvalues of the Obs are used as xerror for the fit.
-889 y : list
-890 List of Obs, the dvalues of the Obs are used as yerror for the fit.
-891
-892 Returns
-893 -------
-894 fit_parameters : list[Obs]
-895 LIist of fitted observables.
-896 """
-897
-898 def f(a, x):
-899 y = a[0] + a[1] * x
-900 return y
-901
-902 if all(isinstance(n, Obs) for n in x):
-903 out = total_least_squares(x, y, f, **kwargs)
-904 return out.fit_parameters
-905 elif all(isinstance(n, float) or isinstance(n, int) for n in x) or isinstance(x, np.ndarray):
-906 out = least_squares(x, y, f, **kwargs)
-907 return out.fit_parameters
-908 else:
-909 raise Exception('Unsupported types for x')
+ 928def fit_lin(x, y, **kwargs):
+929 """Performs a linear fit to y = n + m * x and returns two Obs n, m.
+930
+931 Parameters
+932 ----------
+933 x : list
+934 Can either be a list of floats in which case no xerror is assumed, or
+935 a list of Obs, where the dvalues of the Obs are used as xerror for the fit.
+936 y : list
+937 List of Obs, the dvalues of the Obs are used as yerror for the fit.
+938
+939 Returns
+940 -------
+941 fit_parameters : list[Obs]
+942 LIist of fitted observables.
+943 """
+944
+945 def f(a, x):
+946 y = a[0] + a[1] * x
+947 return y
+948
+949 if all(isinstance(n, Obs) for n in x):
+950 out = total_least_squares(x, y, f, **kwargs)
+951 return out.fit_parameters
+952 elif all(isinstance(n, float) or isinstance(n, int) for n in x) or isinstance(x, np.ndarray):
+953 out = least_squares(x, y, f, **kwargs)
+954 return out.fit_parameters
+955 else:
+956 raise Exception('Unsupported types for x')
912def qqplot(x, o_y, func, p):
-913 """Generates a quantile-quantile plot of the fit result which can be used to
-914 check if the residuals of the fit are gaussian distributed.
-915
-916 Returns
-917 -------
-918 None
-919 """
-920
-921 residuals = []
-922 for i_x, i_y in zip(x, o_y):
-923 residuals.append((i_y - func(p, i_x)) / i_y.dvalue)
-924 residuals = sorted(residuals)
-925 my_y = [o.value for o in residuals]
-926 probplot = scipy.stats.probplot(my_y)
-927 my_x = probplot[0][0]
-928 plt.figure(figsize=(8, 8 / 1.618))
-929 plt.errorbar(my_x, my_y, fmt='o')
-930 fit_start = my_x[0]
-931 fit_stop = my_x[-1]
-932 samples = np.arange(fit_start, fit_stop, 0.01)
-933 plt.plot(samples, samples, 'k--', zorder=11, label='Standard normal distribution')
-934 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='-')
-935
-936 plt.xlabel('Theoretical quantiles')
-937 plt.ylabel('Ordered Values')
-938 plt.legend()
-939 plt.draw()
+ 959def qqplot(x, o_y, func, p):
+960 """Generates a quantile-quantile plot of the fit result which can be used to
+961 check if the residuals of the fit are gaussian distributed.
+962
+963 Returns
+964 -------
+965 None
+966 """
+967
+968 residuals = []
+969 for i_x, i_y in zip(x, o_y):
+970 residuals.append((i_y - func(p, i_x)) / i_y.dvalue)
+971 residuals = sorted(residuals)
+972 my_y = [o.value for o in residuals]
+973 probplot = scipy.stats.probplot(my_y)
+974 my_x = probplot[0][0]
+975 plt.figure(figsize=(8, 8 / 1.618))
+976 plt.errorbar(my_x, my_y, fmt='o')
+977 fit_start = my_x[0]
+978 fit_stop = my_x[-1]
+979 samples = np.arange(fit_start, fit_stop, 0.01)
+980 plt.plot(samples, samples, 'k--', zorder=11, label='Standard normal distribution')
+981 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='-')
+982
+983 plt.xlabel('Theoretical quantiles')
+984 plt.ylabel('Ordered Values')
+985 plt.legend()
+986 plt.draw()
942def residual_plot(x, y, func, fit_res):
-943 """Generates a plot which compares the fit to the data and displays the corresponding residuals
-944
-945 Returns
-946 -------
-947 None
-948 """
-949 sorted_x = sorted(x)
-950 xstart = sorted_x[0] - 0.5 * (sorted_x[1] - sorted_x[0])
-951 xstop = sorted_x[-1] + 0.5 * (sorted_x[-1] - sorted_x[-2])
-952 x_samples = np.arange(xstart, xstop + 0.01, 0.01)
-953
-954 plt.figure(figsize=(8, 8 / 1.618))
-955 gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
-956 ax0 = plt.subplot(gs[0])
-957 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')
-958 ax0.plot(x_samples, func([o.value for o in fit_res], x_samples), label='Fit', zorder=10, ls='-', ms=0)
-959 ax0.set_xticklabels([])
-960 ax0.set_xlim([xstart, xstop])
-961 ax0.set_xticklabels([])
-962 ax0.legend()
-963
-964 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])
-965 ax1 = plt.subplot(gs[1])
-966 ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
-967 ax1.tick_params(direction='out')
-968 ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
-969 ax1.axhline(y=0.0, ls='--', color='k', marker=" ")
-970 ax1.fill_between(x_samples, -1.0, 1.0, alpha=0.1, facecolor='k')
-971 ax1.set_xlim([xstart, xstop])
-972 ax1.set_ylabel('Residuals')
-973 plt.subplots_adjust(wspace=None, hspace=None)
-974 plt.draw()
+ 989def residual_plot(x, y, func, fit_res):
+ 990 """Generates a plot which compares the fit to the data and displays the corresponding residuals
+ 991
+ 992 Returns
+ 993 -------
+ 994 None
+ 995 """
+ 996 sorted_x = sorted(x)
+ 997 xstart = sorted_x[0] - 0.5 * (sorted_x[1] - sorted_x[0])
+ 998 xstop = sorted_x[-1] + 0.5 * (sorted_x[-1] - sorted_x[-2])
+ 999 x_samples = np.arange(xstart, xstop + 0.01, 0.01)
+1000
+1001 plt.figure(figsize=(8, 8 / 1.618))
+1002 gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0)
+1003 ax0 = plt.subplot(gs[0])
+1004 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')
+1005 ax0.plot(x_samples, func([o.value for o in fit_res], x_samples), label='Fit', zorder=10, ls='-', ms=0)
+1006 ax0.set_xticklabels([])
+1007 ax0.set_xlim([xstart, xstop])
+1008 ax0.set_xticklabels([])
+1009 ax0.legend()
+1010
+1011 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])
+1012 ax1 = plt.subplot(gs[1])
+1013 ax1.plot(x, residuals, 'ko', ls='none', markersize=5)
+1014 ax1.tick_params(direction='out')
+1015 ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True)
+1016 ax1.axhline(y=0.0, ls='--', color='k', marker=" ")
+1017 ax1.fill_between(x_samples, -1.0, 1.0, alpha=0.1, facecolor='k')
+1018 ax1.set_xlim([xstart, xstop])
+1019 ax1.set_ylabel('Residuals')
+1020 plt.subplots_adjust(wspace=None, hspace=None)
+1021 plt.draw()
977def error_band(x, func, beta):
-978 """Calculate the error band for an array of sample values x, for given fit function func with optimized parameters beta.
-979
-980 Returns
-981 -------
-982 err : np.array(Obs)
-983 Error band for an array of sample values x
-984 """
-985 cov = covariance(beta)
-986 if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps):
-987 warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning)
-988
-989 deriv = []
-990 for i, item in enumerate(x):
-991 deriv.append(np.array(egrad(func)([o.value for o in beta], item)))
-992
-993 err = []
-994 for i, item in enumerate(x):
-995 err.append(np.sqrt(deriv[i] @ cov @ deriv[i]))
-996 err = np.array(err)
-997
-998 return err
+ 1024def error_band(x, func, beta):
+1025 """Calculate the error band for an array of sample values x, for given fit function func with optimized parameters beta.
+1026
+1027 Returns
+1028 -------
+1029 err : np.array(Obs)
+1030 Error band for an array of sample values x
+1031 """
+1032 cov = covariance(beta)
+1033 if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps):
+1034 warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning)
+1035
+1036 deriv = []
+1037 for i, item in enumerate(x):
+1038 deriv.append(np.array(egrad(func)([o.value for o in beta], item)))
+1039
+1040 err = []
+1041 for i, item in enumerate(x):
+1042 err.append(np.sqrt(deriv[i] @ cov @ deriv[i]))
+1043 err = np.array(err)
+1044
+1045 return err
1001def ks_test(objects=None):
-1002 """Performs a Kolmogorov–Smirnov test for the p-values of all fit object.
-1003
-1004 Parameters
-1005 ----------
-1006 objects : list
-1007 List of fit results to include in the analysis (optional).
-1008
-1009 Returns
-1010 -------
-1011 None
-1012 """
-1013
-1014 if objects is None:
-1015 obs_list = []
-1016 for obj in gc.get_objects():
-1017 if isinstance(obj, Fit_result):
-1018 obs_list.append(obj)
-1019 else:
-1020 obs_list = objects
-1021
-1022 p_values = [o.p_value for o in obs_list]
-1023
-1024 bins = len(p_values)
-1025 x = np.arange(0, 1.001, 0.001)
-1026 plt.plot(x, x, 'k', zorder=1)
-1027 plt.xlim(0, 1)
-1028 plt.ylim(0, 1)
-1029 plt.xlabel('p-value')
-1030 plt.ylabel('Cumulative probability')
-1031 plt.title(str(bins) + ' p-values')
-1032
-1033 n = np.arange(1, bins + 1) / np.float64(bins)
-1034 Xs = np.sort(p_values)
-1035 plt.step(Xs, n)
-1036 diffs = n - Xs
-1037 loc_max_diff = np.argmax(np.abs(diffs))
-1038 loc = Xs[loc_max_diff]
-1039 plt.annotate('', xy=(loc, loc), xytext=(loc, loc + diffs[loc_max_diff]), arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
-1040 plt.draw()
-1041
-1042 print(scipy.stats.kstest(p_values, 'uniform'))
+ 1048def ks_test(objects=None):
+1049 """Performs a Kolmogorov–Smirnov test for the p-values of all fit object.
+1050
+1051 Parameters
+1052 ----------
+1053 objects : list
+1054 List of fit results to include in the analysis (optional).
+1055
+1056 Returns
+1057 -------
+1058 None
+1059 """
+1060
+1061 if objects is None:
+1062 obs_list = []
+1063 for obj in gc.get_objects():
+1064 if isinstance(obj, Fit_result):
+1065 obs_list.append(obj)
+1066 else:
+1067 obs_list = objects
+1068
+1069 p_values = [o.p_value for o in obs_list]
+1070
+1071 bins = len(p_values)
+1072 x = np.arange(0, 1.001, 0.001)
+1073 plt.plot(x, x, 'k', zorder=1)
+1074 plt.xlim(0, 1)
+1075 plt.ylim(0, 1)
+1076 plt.xlabel('p-value')
+1077 plt.ylabel('Cumulative probability')
+1078 plt.title(str(bins) + ' p-values')
+1079
+1080 n = np.arange(1, bins + 1) / np.float64(bins)
+1081 Xs = np.sort(p_values)
+1082 plt.step(Xs, n)
+1083 diffs = n - Xs
+1084 loc_max_diff = np.argmax(np.abs(diffs))
+1085 loc = Xs[loc_max_diff]
+1086 plt.annotate('', xy=(loc, loc), xytext=(loc, loc + diffs[loc_max_diff]), arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
+1087 plt.draw()
+1088
+1089 print(scipy.stats.kstest(p_values, 'uniform'))