From bc190e3bb16ef5fb082fe352d8fbe545503dc618 Mon Sep 17 00:00:00 2001 From: Iahn Cajigas Date: Mon, 16 Mar 2026 21:09:42 -0400 Subject: [PATCH] Fix hybrid decoder likelihood: remove 1e-15 floor that killed high-dim models The determinant ratio sqrt(det(W_u))/sqrt(det(W_p)) in the Laplace approximation used max(..., 1e-15) as a denominator floor. For 6D state models the absolute determinant (~1e-46) falls far below this floor, artificially deflating the likelihood by ~8 orders of magnitude and preventing the filter from ever selecting the reach model. Compute the ratio as sqrt(det(W_u)/det(W_p)) instead, which cancels the absolute scale and matches MATLAB behavior exactly. Co-Authored-By: Claude Opus 4.6 --- nstat/decoding_algorithms.py | 26 ++++++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/nstat/decoding_algorithms.py b/nstat/decoding_algorithms.py index a37c6e3e..a8245eaa 100644 --- a/nstat/decoding_algorithms.py +++ b/nstat/decoding_algorithms.py @@ -2043,7 +2043,18 @@ def PPHybridFilterLinear( X_u[model_index][:, time_index] = upd_x W_u[model_index][:, :, time_index] = upd_W - det_ratio = np.sqrt(max(np.linalg.det(upd_W), 0.0)) / max(np.sqrt(max(np.linalg.det(pred_W), 0.0)), 1e-15) + # Likelihood via Laplace approximation (Srinivasan et al. 2007). + # Compute sqrt(det(W_u)) / sqrt(det(W_p)) as sqrt(det(W_u)/det(W_p)) + # to avoid a fixed floor that destroys the ratio for + # high-dimensional models with tiny absolute determinants. + det_upd = np.linalg.det(upd_W) + det_pred = np.linalg.det(pred_W) + if det_pred > 0.0 and det_upd >= 0.0: + det_ratio = np.sqrt(det_upd / det_pred) + elif det_upd == 0.0 and det_pred == 0.0: + det_ratio = 1.0 + else: + det_ratio = 0.0 log_term = np.sum(obs[:, time_index] * np.log(np.clip(lambda_delta.reshape(-1), 1e-12, np.inf)) - lambda_delta.reshape(-1)) likelihoods[model_index] = float(det_ratio * np.exp(np.clip(log_term, -200.0, 50.0))) @@ -2212,7 +2223,18 @@ def PPHybridFilter(A, Q, p_ij, Mu0, dN, lambdaCIFColl, binwidth=0.001, x0=None, X_u[model_index][:, time_index] = upd_x W_u[model_index][:, :, time_index] = upd_W - det_ratio = np.sqrt(max(np.linalg.det(upd_W), 0.0)) / max(np.sqrt(max(np.linalg.det(pred_W), 0.0)), 1e-15) + # Likelihood via Laplace approximation (Srinivasan et al. 2007). + # Compute sqrt(det(W_u)) / sqrt(det(W_p)) as sqrt(det(W_u)/det(W_p)) + # to avoid a fixed floor that destroys the ratio for + # high-dimensional models with tiny absolute determinants. + det_upd = np.linalg.det(upd_W) + det_pred = np.linalg.det(pred_W) + if det_pred > 0.0 and det_upd >= 0.0: + det_ratio = np.sqrt(det_upd / det_pred) + elif det_upd == 0.0 and det_pred == 0.0: + det_ratio = 1.0 + else: + det_ratio = 0.0 log_term = np.sum(obs[:, time_index] * np.log(np.clip(lambda_delta.reshape(-1), 1e-12, np.inf)) - lambda_delta.reshape(-1)) likelihoods[model_index] = float(det_ratio * np.exp(np.clip(log_term, -200.0, 50.0)))