From 009975029522e4688fb78fcf00d84d4739ab59a8 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Thu, 3 Mar 2022 18:29:18 +0000 Subject: [PATCH 1/2] fix: replaced inverse by the pseudo inverse in the calculation of chi_exp --- pyerrors/fits.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyerrors/fits.py b/pyerrors/fits.py index 6a428e65..c023dd29 100644 --- a/pyerrors/fits.py +++ b/pyerrors/fits.py @@ -249,7 +249,7 @@ def total_least_squares(x, y, func, silent=False, **kwargs): new_jac = np.concatenate((fused_row1, fused_row2), axis=1) A = W @ new_jac - P_phi = A @ np.linalg.inv(A.T @ A) @ A.T + P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T expected_chisquare = np.trace((np.identity(P_phi.shape[0]) - P_phi) @ W @ cov @ W) if expected_chisquare <= 0.0: warnings.warn("Negative expected_chisquare.", RuntimeWarning) @@ -528,7 +528,7 @@ def _standard_fit(x, y, func, silent=False, **kwargs): W = np.diag(1 / np.asarray(dy_f)) cov = covariance(y) A = W @ jacobian(func)(fit_result.x, x) - P_phi = A @ np.linalg.inv(A.T @ A) @ A.T + P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T expected_chisquare = np.trace((np.identity(x.shape[-1]) - P_phi) @ W @ cov @ W) output.chisquare_by_expected_chisquare = chisquare / expected_chisquare if not silent: From 99bead518c6c50c60952ebb07d9e88685c739237 Mon Sep 17 00:00:00 2001 From: Fabian Joswig Date: Mon, 7 Mar 2022 11:54:09 +0000 Subject: [PATCH 2/2] tests: additional test for expected chiqsquare and correlated fit based on cov obs added. --- tests/fits_test.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/tests/fits_test.py b/tests/fits_test.py index da43b4ec..cc958dbc 100644 --- a/tests/fits_test.py +++ b/tests/fits_test.py @@ -319,6 +319,23 @@ def test_prior_fit(): fitp = pe.fits.least_squares([0, 1], y, f, priors=y, resplot=True, qqplot=True) +def test_correlated_fit_covobs(): + x1 = pe.cov_Obs(1.01, 0.01 ** 2, 'test1') + x2 = pe.cov_Obs(2.01, 0.01 ** 2, 'test2') + x3 = pe.cov_Obs(2.99, 0.01 ** 2, 'test3') + + [o.gamma_method() for o in [x1, x2, x3]] + + def func(a, x): + return a[0] * x + a[1] + + fit_res = pe.fits.least_squares(np.arange(1, 4), [x1, x2, x3], func, expected_chisquare=True) + assert np.isclose(fit_res.chisquare_by_dof, fit_res.chisquare_by_expected_chisquare) + + fit_res_corr = pe.fits.least_squares(np.arange(1, 4), [x1, x2, x3], func, correlated_fit=True) + assert np.isclose(fit_res.chisquare_by_dof, fit_res_corr.chisquare_by_dof) + + def test_error_band(): def f(a, x): return a[0] + a[1] * x