gwm17/spspy
1
0
Fork 0
mirror of https://github.com/gwm17/spspy.git synced 2025-07-20 03:58:50 -04:00
Code Issues Projects Releases Wiki Activity

Fix bug in printed Chi-Sq. value

Browse Source
This commit is contained in:
gwm17 2025-07-15 07:33:12 -05:00
parent ae5724afbc
commit 97acc1659b
1 changed files with 34 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

53
spspy/Fitter.py
View File

@ -8,6 +8,7 @@ from typing import Optional
INVALID_FIT_RESULT = np.inf INVALID_FIT_RESULT = np.inf
INVALID_NDF = -1 INVALID_NDF = -1
@dataclass @dataclass
class FitPoint: class FitPoint:
x: float = 0.0 x: float = 0.0
@ -15,7 +16,12 @@ class FitPoint:
xError: float = 0.0 xError: float = 0.0
yError: 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)) xArray = np.empty(len(data))
yArray = np.empty(len(data)) yArray = np.empty(len(data))
xErrorArray = 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 yErrorArray[index] = point.yError
return xArray, yArray, xErrorArray, yErrorArray return xArray, yArray, xErrorArray, yErrorArray
@dataclass @dataclass
class FitResidual: class FitResidual:
x: float = 0.0 x: float = 0.0
residual: float = 0.0 residual: float = 0.0
studentizedResidual: 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)) xArray = np.empty(len(data))
residArray = np.empty(len(data)) residArray = np.empty(len(data))
studentResidArray = 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 studentResidArray[index] = point.studentizedResidual
return xArray, residArray, studentResidArray return xArray, residArray, studentResidArray
class Fitter: class Fitter:
def __init__(self, order: int =1): def __init__(self, order: int = 1):
self.polynomialOrder: int = order self.polynomialOrder: int = order
self.fitResults: Optional[odr.Output] = None self.fitResults: Optional[odr.Output] = None
self.fitData: Optional[list[FitPoint]] = None self.fitData: Optional[list[FitPoint]] = None
@ -58,7 +69,9 @@ class Fitter:
self.fitData = data self.fitData = data
if self.fitData is not None: 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) modelData = odr.RealData(xArray, y=yArray, sx=xErrorArray, sy=yErrorArray)
model = odr.polynomial(self.polynomialOrder) model = odr.polynomial(self.polynomialOrder)
self.fitResults = odr.ODR(modelData, model).run() self.fitResults = odr.ODR(modelData, model).run()
@ -66,12 +79,12 @@ class Fitter:
else: else:
print("Cannot run fitter without setting data to be fit!") 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: if self.fitResults is not None:
return self.fitResults.beta return self.fitResults.beta
return np.array({INVALID_FIT_RESULT}) 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: if self.fitResults is not None:
return self.fitResults.sd_beta return self.fitResults.sd_beta
return np.array({INVALID_FIT_RESULT}) return np.array({INVALID_FIT_RESULT})
@ -97,27 +110,23 @@ class Fitter:
return INVALID_FIT_RESULT return INVALID_FIT_RESULT
def get_chisquare(self) -> float: def get_chisquare(self) -> float:
if self.function is not None: if self.fitResults is 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 return INVALID_FIT_RESULT
return self.fitResults.res_var * self.get_ndf()
def get_reduced_chisquare(self) -> float: def get_reduced_chisquare(self) -> float:
ndf = self.get_ndf() if self.fitResults is None:
chisq = self.get_chisquare()
if chisq == INVALID_FIT_RESULT or ndf == INVALID_NDF:
return INVALID_FIT_RESULT return INVALID_FIT_RESULT
else: return self.fitResults.res_var
return chisq/ndf
def get_residuals(self) -> list[FitResidual]: def get_residuals(self) -> list[FitResidual]:
if self.fitData is not None: 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 xMean = 0.0
rmse = 0.0 rmse = 0.0
npoints = len(fitResiduals) npoints = len(fitResiduals)
@ -132,7 +141,9 @@ class Fitter:
meanDiffSq += (resid.x - xMean) ** 2.0 meanDiffSq += (resid.x - xMean) ** 2.0
meanDiffSq /= npoints meanDiffSq /= npoints
for resid in fitResiduals: for resid in fitResiduals:
leverage = 1.0/npoints + (resid.x - xMean)/meanDiffSq leverage = 1.0 / npoints + (resid.x - xMean) / meanDiffSq
resid.studentizedResidual = resid.residual / (rmse * np.sqrt(1.0 - leverage)) resid.studentizedResidual = resid.residual / (
rmse * np.sqrt(1.0 - leverage)
)
return fitResiduals return fitResiduals
return [] return []