From 97acc1659b546b999deea79a9a7c4bb7b66a66ef Mon Sep 17 00:00:00 2001 From: gwm17 Date: Tue, 15 Jul 2025 07:33:12 -0500 Subject: [PATCH] Fix bug in printed Chi-Sq. value --- spspy/Fitter.py | 57 +++++++++++++++++++++++++++++-------------------- 1 file changed, 34 insertions(+), 23 deletions(-) diff --git a/spspy/Fitter.py b/spspy/Fitter.py index a252d23..1b77975 100644 --- a/spspy/Fitter.py +++ b/spspy/Fitter.py @@ -8,6 +8,7 @@ from typing import Optional INVALID_FIT_RESULT = np.inf INVALID_NDF = -1 + @dataclass class FitPoint: x: float = 0.0 @@ -15,7 +16,12 @@ class FitPoint: xError: float = 0.0 yError: float = 0.0 -def convert_fit_points_to_arrays(data: list[FitPoint]) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]: + +def convert_fit_points_to_arrays( + data: list[FitPoint], +) -> tuple[ + NDArray[np.float64], NDArray[np.float64], NDArray[np.float64], NDArray[np.float64] +]: xArray = np.empty(len(data)) yArray = np.empty(len(data)) xErrorArray = np.empty(len(data)) @@ -27,13 +33,17 @@ def convert_fit_points_to_arrays(data: list[FitPoint]) -> tuple[NDArray[np.float yErrorArray[index] = point.yError return xArray, yArray, xErrorArray, yErrorArray + @dataclass class FitResidual: x: float = 0.0 residual: float = 0.0 studentizedResidual: float = 0.0 -def convert_resid_points_to_arrays(data: list[FitResidual]) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]: + +def convert_resid_points_to_arrays( + data: list[FitResidual], +) -> tuple[NDArray[np.float64], NDArray[np.float64], NDArray[np.float64]]: xArray = np.empty(len(data)) residArray = np.empty(len(data)) studentResidArray = np.empty(len(data)) @@ -43,8 +53,9 @@ def convert_resid_points_to_arrays(data: list[FitResidual]) -> tuple[NDArray[np. studentResidArray[index] = point.studentizedResidual return xArray, residArray, studentResidArray + class Fitter: - def __init__(self, order: int =1): + def __init__(self, order: int = 1): self.polynomialOrder: int = order self.fitResults: Optional[odr.Output] = None self.fitData: Optional[list[FitPoint]] = None @@ -56,9 +67,11 @@ class Fitter: def run(self, data: list[FitPoint] = None) -> None: if data is not None: self.fitData = data - + if self.fitData is not None: - xArray, yArray, xErrorArray, yErrorArray = convert_fit_points_to_arrays(self.fitData) + xArray, yArray, xErrorArray, yErrorArray = convert_fit_points_to_arrays( + self.fitData + ) modelData = odr.RealData(xArray, y=yArray, sx=xErrorArray, sy=yErrorArray) model = odr.polynomial(self.polynomialOrder) self.fitResults = odr.ODR(modelData, model).run() @@ -66,12 +79,12 @@ class Fitter: else: print("Cannot run fitter without setting data to be fit!") - def get_parameters(self) -> NDArray[np.float64] : + def get_parameters(self) -> NDArray[np.float64]: if self.fitResults is not None: return self.fitResults.beta return np.array({INVALID_FIT_RESULT}) - def get_parameter_errors(self) -> NDArray[np.float64] : + def get_parameter_errors(self) -> NDArray[np.float64]: if self.fitResults is not None: return self.fitResults.sd_beta return np.array({INVALID_FIT_RESULT}) @@ -97,27 +110,23 @@ class Fitter: return INVALID_FIT_RESULT def get_chisquare(self) -> float: - if self.function is not None: - yEffErrorArray = [np.sqrt(point.yError**2.0 + (point.xError * self.evaluate_derivative(point.x))**2.0) for point in self.fitData] - chisq = 0.0 - for index, point in enumerate(self.fitData): - chisq = ((point.y - self.evaluate(point.x)) / yEffErrorArray[index])**2.0 - return chisq - return INVALID_FIT_RESULT + if self.fitResults is None: + return INVALID_FIT_RESULT + return self.fitResults.res_var * self.get_ndf() def get_reduced_chisquare(self) -> float: - ndf = self.get_ndf() - chisq = self.get_chisquare() - if chisq == INVALID_FIT_RESULT or ndf == INVALID_NDF: + if self.fitResults is None: return INVALID_FIT_RESULT - else: - return chisq/ndf + return self.fitResults.res_var def get_residuals(self) -> list[FitResidual]: if self.fitData is not None: - fitResiduals = [FitResidual(point.x, point.y - self.evaluate(point.x), 0.0) for point in self.fitData] + fitResiduals = [ + FitResidual(point.x, point.y - self.evaluate(point.x), 0.0) + for point in self.fitData + ] - #compute the leverage and studentize + # compute the leverage and studentize xMean = 0.0 rmse = 0.0 npoints = len(fitResiduals) @@ -132,7 +141,9 @@ class Fitter: meanDiffSq += (resid.x - xMean) ** 2.0 meanDiffSq /= npoints for resid in fitResiduals: - leverage = 1.0/npoints + (resid.x - xMean)/meanDiffSq - resid.studentizedResidual = resid.residual / (rmse * np.sqrt(1.0 - leverage)) + leverage = 1.0 / npoints + (resid.x - xMean) / meanDiffSq + resid.studentizedResidual = resid.residual / ( + rmse * np.sqrt(1.0 - leverage) + ) return fitResiduals return []