Better treatment for numerical errors

LAMPSPUC · Oct 9, 2024 · 193f6ae · 193f6ae
1 parent 6a86265
commit 193f6ae
Show file tree

Hide file tree

Showing 5 changed files with 22 additions and 5 deletions.
diff --git a/src/filters/multivariate_kalman_filter.jl b/src/filters/multivariate_kalman_filter.jl
@@ -187,9 +187,14 @@ function update_P!(
 end
 
 function update_llk!(kalman_state::MultivariateKalmanState{Fl}) where Fl
-    kalman_state.llk -=
-        HALF_LOG_2_PI + 0.5 * (logdet(kalman_state.F) +
-        kalman_state.v' * inv(kalman_state.F) * kalman_state.v)
+    detF = det(kalman_state.F)
+    if detF < 0
+        error("Numerical error, the determinant of F $(kalman_state.F) is negative: $detF")
+    else
+        kalman_state.llk -=
+            HALF_LOG_2_PI + 0.5 * (log(detF) +
+            kalman_state.v' * inv(kalman_state.F) * kalman_state.v)
+    end
     return kalman_state
 end
 

diff --git a/src/filters/regression_kalman_filter.jl b/src/filters/regression_kalman_filter.jl
@@ -58,6 +58,9 @@ function regression_update_v!(
 end
 
 function update_llk!(kalman_state::RegressionKalmanState{Fl}) where Fl
+    if kalman_state.F < 0
+        error("Numerical error, F is negative: $(kalman_state.F)")
+    end
     kalman_state.llk -= (
         HALF_LOG_2_PI + 0.5 * (log(kalman_state.F) + kalman_state.v^2 / kalman_state.F)
     )

diff --git a/src/filters/scalar_kalman_filter.jl b/src/filters/scalar_kalman_filter.jl
@@ -125,6 +125,9 @@ function scalar_update_P!(kalman_state::ScalarKalmanState{Fl}, T::Fl, RQR::Fl) w
 end
 
 function update_llk!(kalman_state::ScalarKalmanState{Fl}) where Fl
+    if kalman_state.F < 0
+        error("Numerical error, F is negative: $(kalman_state.F)")
+    end
     kalman_state.llk -= (
         HALF_LOG_2_PI + 0.5 * (log(kalman_state.F) + kalman_state.v^2 / kalman_state.F)
     )

diff --git a/src/filters/sparse_univariate_kalman_filter.jl b/src/filters/sparse_univariate_kalman_filter.jl
@@ -179,8 +179,11 @@ function update_P!(
 end
 
 function update_llk!(kalman_state::SparseUnivariateKalmanState{Fl}) where Fl
+    if kalman_state.F < 0
+        error("Numerical error, F is negative: $(kalman_state.F)")
+    end
     kalman_state.llk -= (
-        HALF_LOG_2_PI + (log(kalman_state.F) + kalman_state.v^2 / kalman_state.F) / 2
+        HALF_LOG_2_PI + 0.5 * (log(kalman_state.F) + kalman_state.v^2 / kalman_state.F)
     )
     return kalman_state
 end

diff --git a/src/filters/univariate_kalman_filter.jl b/src/filters/univariate_kalman_filter.jl
@@ -195,8 +195,11 @@ function update_P!(
 end
 
 function update_llk!(kalman_state::UnivariateKalmanState{Fl}) where Fl
+    if kalman_state.F < 0
+        error("Numerical error, F is negative: $(kalman_state.F)")
+    end
     kalman_state.llk -= (
-        HALF_LOG_2_PI + (log(kalman_state.F) + kalman_state.v^2 / kalman_state.F) / 2
+        HALF_LOG_2_PI + 0.5 * (log(kalman_state.F) + kalman_state.v^2 / kalman_state.F)
     )
     return kalman_state
 end