diff --git a/particula_beta/data/process/chamber_rate_fitting.py b/particula_beta/data/process/chamber_rate_fitting.py index 9d68d7f..1847b14 100644 --- a/particula_beta/data/process/chamber_rate_fitting.py +++ b/particula_beta/data/process/chamber_rate_fitting.py @@ -47,17 +47,19 @@ def chi_squared_error( def hessian_standard_error( objective_function: Callable, - parameters: NDArray[np.float64], + parameters: np.ndarray, epsilon: float = 1e-4, -) -> Tuple[float, float, float]: + regularization: float = 1e-8, +) -> np.ndarray: """ Calculate the standard error for the optimized parameters using Hessian - estimation. + estimation with regularization and fallback to diagonal approximation. Arguments: - func: The objective function. - params: Array of optimized parameters. + objective_function: The objective function. + parameters: Array of optimized parameters. epsilon: Step size for numerical differentiation. + regularization: Regularization term for Hessian inversion. Returns: standard_errors: Standard errors of the optimized parameters. @@ -69,7 +71,6 @@ def hessian_standard_error( for i in range(n): for j in range(n): params_ij = parameters.copy() - if ( i == j ): # Diagonal elements: second derivative w.r.t. one parameter @@ -99,10 +100,18 @@ def hessian_standard_error( 4 * epsilon**2 ) - # Compute covariance matrix and standard errors - cov_matrix = np.linalg.inv(hessian) - standard_errors = np.sqrt(np.diag(cov_matrix)) - return standard_errors[0], standard_errors[1], standard_errors[2] + # Regularization + hessian += np.eye(n) * regularization + + try: + cov_matrix = np.linalg.inv(hessian) + standard_errors = np.sqrt(np.diag(cov_matrix)) + except np.linalg.LinAlgError: + # Fallback to diagonal approximation + print("Hessian inversion failed; using diagonal approximation.") + standard_errors = np.sqrt(1.0 / np.diag(hessian)) + + return standard_errors # pylint: disable=too-many-positional-arguments, too-many-arguments