From 9c86df074f3f771a2a95102b95ef70f5dc4e96da Mon Sep 17 00:00:00 2001 From: Aleksandr Tarutin Date: Fri, 9 Oct 2026 19:28:59 +0000 Subject: [PATCH 1/3] Fix NaN in hmm_hidden_state_prob for long sequences (#2677) The forward pass rescaled only its first column, so on a long sequence the forward probabilities underflowed to 0 and the final normalization computed 0/0. Rescale every forward step, as hmm_latent_rng and hmm_marginal already do. The hidden state probabilities are invariant to a positive factor on each column, so the scale does not need to be kept. log_omegas was also exponentiated directly, so log densities below about -745 underflowed to 0 and large ones overflowed to inf. Subtract the maximum of each column before exponentiating; the hidden state probabilities are invariant to that shift as well. The test hidden_state_prob1 is taken from Charles Margossian's branch bugfix/issue-2677-hmm_hidden_state (e7a5efa). --- stan/math/prim/prob/hmm_hidden_state_prob.hpp | 19 ++- .../prim/prob/hmm_hidden_state_prob_test.cpp | 141 ++++++++++++++++++ 2 files changed, 158 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/prob/hmm_hidden_state_prob.hpp b/stan/math/prim/prob/hmm_hidden_state_prob.hpp index d90731d23ad..768e27b4d8e 100644 --- a/stan/math/prim/prob/hmm_hidden_state_prob.hpp +++ b/stan/math/prim/prob/hmm_hidden_state_prob.hpp @@ -29,6 +29,10 @@ namespace math { * of a state sequence given parameters and an observation sequence, * p(x | y, theta), * because it only computes marginals on a state-by-state basis. + * Each column of log_omegas is shifted by its maximum before it is + * exponentiated, and the forward and backward passes are normalized at + * every step, so the result neither underflows nor overflows for long + * sequences or for log densities of large magnitude. * * @tparam T_omega type of the log likelihood matrix * @tparam T_Gamma type of the transition matrix @@ -51,7 +55,15 @@ inline Eigen::MatrixXd hmm_hidden_state_prob(const T_omega& log_omegas, int n_states = log_omegas.rows(); int n_transitions = log_omegas.cols() - 1; - Eigen::MatrixXd omegas = value_of(log_omegas).array().exp(); + // Subtract the maximum of each column before exponentiating, so that + // omegas neither underflows nor overflows. The hidden state probabilities + // are invariant to a positive factor on each column of omegas. + ref_type_t log_omegas_dbl + = value_of(log_omegas); + Eigen::MatrixXd omegas + = (log_omegas_dbl.rowwise() - log_omegas_dbl.colwise().maxCoeff()) + .array() + .exp(); ref_type_t rho_dbl = value_of(rho); ref_type_t Gamma_dbl = value_of(Gamma); hmm_check(log_omegas, Gamma_dbl, rho_dbl, "hmm_hidden_state_prob"); @@ -60,9 +72,12 @@ inline Eigen::MatrixXd hmm_hidden_state_prob(const T_omega& log_omegas, alphas.col(0) = omegas.col(0).cwiseProduct(rho_dbl); alphas.col(0) /= alphas.col(0).maxCoeff(); - for (int n = 0; n < n_transitions; ++n) + // Forward pass with running normalization + for (int n = 0; n < n_transitions; ++n) { alphas.col(n + 1) = omegas.col(n + 1).cwiseProduct(Gamma_dbl.transpose() * alphas.col(n)); + alphas.col(n + 1) /= alphas.col(n + 1).maxCoeff(); + } // Backward pass with running normalization Eigen::VectorXd beta = Eigen::VectorXd::Ones(n_states); diff --git a/test/unit/math/prim/prob/hmm_hidden_state_prob_test.cpp b/test/unit/math/prim/prob/hmm_hidden_state_prob_test.cpp index 747784ab5b9..122bd626f09 100644 --- a/test/unit/math/prim/prob/hmm_hidden_state_prob_test.cpp +++ b/test/unit/math/prim/prob/hmm_hidden_state_prob_test.cpp @@ -1,13 +1,56 @@ #include #include +#include #include #include #include #include #include +#include #include #include +namespace { +/** + * Reference implementation of hmm_hidden_state_prob: the forward-backward + * algorithm computed entirely on the log scale, with log_sum_exp in place + * of every sum, so that neither the forward nor the backward pass can + * underflow or overflow. + */ +Eigen::MatrixXd hmm_hidden_state_prob_log_space( + const Eigen::MatrixXd& log_omegas, const Eigen::MatrixXd& Gamma, + const Eigen::VectorXd& rho) { + using stan::math::log_sum_exp; + const int n_states = log_omegas.rows(); + const int n_cols = log_omegas.cols(); + const Eigen::MatrixXd log_Gamma = Gamma.array().log(); + + Eigen::MatrixXd log_alphas(n_states, n_cols); + log_alphas.col(0) = rho.array().log() + log_omegas.col(0).array(); + for (int n = 1; n < n_cols; ++n) { + for (int j = 0; j < n_states; ++j) { + Eigen::VectorXd terms = log_alphas.col(n - 1) + log_Gamma.col(j); + log_alphas(j, n) = log_omegas(j, n) + log_sum_exp(terms); + } + } + + Eigen::MatrixXd prob(n_states, n_cols); + Eigen::VectorXd log_beta = Eigen::VectorXd::Zero(n_states); + for (int n = n_cols - 1; n >= 0; --n) { + if (n < n_cols - 1) { + Eigen::VectorXd next = log_omegas.col(n + 1) + log_beta; + for (int i = 0; i < n_states; ++i) { + Eigen::VectorXd terms = log_Gamma.row(i).transpose() + next; + log_beta(i) = log_sum_exp(terms); + } + } + Eigen::VectorXd log_prob = log_alphas.col(n) + log_beta; + prob.col(n) = (log_prob.array() - log_sum_exp(log_prob)).exp(); + } + return prob; +} +} // namespace + TEST_F(hmm_test, hidden_state_single_outcome) { using stan::math::hmm_hidden_state_prob; @@ -69,3 +112,101 @@ TEST(hmm_test_nonstandard, hidden_state_symmetry) { EXPECT_FLOAT_EQ(prob(1, i), 0.5); } } + +TEST(hmm_test_nonstandard, hidden_state_prob1) { + // This time, the transition matrix forces states to transition + // to the first state with probability 1. The large log densities + // overflow if they are exponentiated directly (#2677). + using stan::math::hmm_hidden_state_prob; + int n_states = 2; + int n_transitions = 2; + Eigen::MatrixXd Gamma(n_states, n_states); + Gamma << 1, 0, 1, 0; + Eigen::VectorXd rho(n_states); + rho << 0.3, 0.7; + Eigen::MatrixXd log_omegas + = 1000 * Eigen::MatrixXd::Ones(n_states, n_transitions + 1); + + Eigen::MatrixXd prob = hmm_hidden_state_prob(log_omegas, Gamma, rho); + EXPECT_FLOAT_EQ(prob(0, 0), 0.3); + EXPECT_FLOAT_EQ(prob(1, 0), 0.7); + + for (int i = 1; i < n_transitions; i++) { + EXPECT_FLOAT_EQ(prob(0, i), 1); + EXPECT_FLOAT_EQ(prob(1, i), 0); + } +} + +TEST_F(hmm_test, hidden_state_very_negative_log_omegas) { + // Regression test for #2677. Every log density is below -1000, so + // exponentiating it directly underflows to 0. Adding a constant to a + // column of log_omegas does not change the hidden state probabilities. + using stan::math::hmm_hidden_state_prob; + Eigen::MatrixXd log_omegas = log_omegas_; + for (int n = 0; n < log_omegas.cols(); ++n) { + log_omegas.col(n).array() -= 1000 + 10 * n; + } + + Eigen::MatrixXd prob = hmm_hidden_state_prob(log_omegas, Gamma_, rho_); + Eigen::MatrixXd prob_unshifted + = hmm_hidden_state_prob(log_omegas_, Gamma_, rho_); + Eigen::MatrixXd prob_ref + = hmm_hidden_state_prob_log_space(log_omegas, Gamma_, rho_); + + for (int n = 0; n < prob.cols(); ++n) { + for (int k = 0; k < prob.rows(); ++k) { + EXPECT_TRUE(std::isfinite(prob(k, n))); + EXPECT_NEAR(prob(k, n), prob_unshifted(k, n), 1e-12); + // The log scale reference loses about machine epsilon times the + // magnitude of the log densities, hence the looser tolerance. + EXPECT_NEAR(prob(k, n), prob_ref(k, n), 1e-10); + } + } +} + +TEST(hmm_test_nonstandard, hidden_state_long_sequence) { + // Regression test for #2677. The log density of a long sequence is far + // below the smallest double, so the forward pass underflows to 0 unless + // it is rescaled at every step. + using stan::math::hmm_hidden_state_prob; + using stan::math::hmm_marginal; + int n_states = 2; + int n_obs = 5000; + Eigen::MatrixXd Gamma(n_states, n_states); + Gamma << 0.95, 0.05, 0.10, 0.90; + Eigen::VectorXd rho(n_states); + rho << 2.0 / 3, 1.0 / 3; + + // Simulate a two-state Gaussian hidden Markov model, with observations + // normal(1, 1) in state 0 and normal(-1, 1) in state 1. + boost::random::mt19937 rng(2677); + boost::random::uniform_01<> uniform; + boost::random::normal_distribution<> normal(0, 1); + Eigen::MatrixXd log_omegas(n_states, n_obs); + int state = 0; + for (int n = 0; n < n_obs; ++n) { + state = uniform(rng) < Gamma(state, 0) ? 0 : 1; + double y = (state == 0 ? 1 : -1) + normal(rng); + for (int k = 0; k < n_states; ++k) { + log_omegas(k, n) = state_lpdf(y, 1, 1, k); + } + } + // The joint density underflows the double range many times over. + EXPECT_LT(hmm_marginal(log_omegas, Gamma, rho), -2000); + + Eigen::MatrixXd prob = hmm_hidden_state_prob(log_omegas, Gamma, rho); + Eigen::MatrixXd prob_ref + = hmm_hidden_state_prob_log_space(log_omegas, Gamma, rho); + + int first_non_finite = -1; + for (int n = 0; n < n_obs && first_non_finite < 0; ++n) { + if (!prob.col(n).allFinite()) { + first_non_finite = n; + } + } + ASSERT_EQ(first_non_finite, -1) << "first non-finite column"; + EXPECT_LT((prob.colwise().sum().array() - 1).abs().maxCoeff(), 1e-12); + // The log scale reference loses about machine epsilon times the + // magnitude of the log densities, hence the looser tolerance. + EXPECT_LT((prob - prob_ref).cwiseAbs().maxCoeff(), 1e-10); +} From 7ceb00425eac475033745b308b1e4cb459111d15 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Wed, 28 Oct 2026 13:44:20 -0400 Subject: [PATCH 2/3] Shift log_omegas in hmm_latent_rng and hmm_marginal (#2677) hmm_latent_rng and hmm_marginal exponentiated log_omegas directly, so log densities below about -745 underflowed to 0 and large ones overflowed to inf, as in hmm_hidden_state_prob before. Subtract the maximum of each column before exponentiating. The sampling probabilities are invariant to the shift. hmm_marginal adds the shifts back to the log density, and its gradients are unchanged because the shift cancels in the adjoint. In hmm_hidden_state_prob, call hmm_check before taking the column maxima, so that a log_omegas with no rows throws instead of hitting an Eigen assert. Fix the docblocks of hmm_hidden_state_prob and hmm_latent_rng, which described log_omegas transposed, and state that the result can still be NaN when the largest density in a column belongs to a state the chain cannot reach. The loops in four hidden state tests skipped the last column. --- stan/math/prim/prob/hmm_hidden_state_prob.hpp | 17 +++++---- stan/math/prim/prob/hmm_latent_rng.hpp | 15 ++++++-- stan/math/prim/prob/hmm_marginal.hpp | 13 +++++-- .../prim/prob/hmm_hidden_state_prob_test.cpp | 35 +++++++++++++++-- .../math/prim/prob/hmm_latent_rng_test.cpp | 21 ++++++++++ .../unit/math/prim/prob/hmm_marginal_test.cpp | 38 +++++++++++++++++++ 6 files changed, 122 insertions(+), 17 deletions(-) diff --git a/stan/math/prim/prob/hmm_hidden_state_prob.hpp b/stan/math/prim/prob/hmm_hidden_state_prob.hpp index 768e27b4d8e..0c44751a4f2 100644 --- a/stan/math/prim/prob/hmm_hidden_state_prob.hpp +++ b/stan/math/prim/prob/hmm_hidden_state_prob.hpp @@ -20,8 +20,8 @@ namespace math { * Hence for each hidden variable x, we compute a simplex with K elements. * The final result is stored in a K by N matrix, where N is the length of x. * log_omegas is a matrix of observational densities, where - * the (i, j)th entry corresponds to the density of the ith observation, y_i, - * given x_i = j. + * the (i, j)th entry corresponds to the density of the jth observation, y_j, + * given x_j = i. * The transition matrix Gamma is such that the (i, j)th entry is the * probability that x_n = j given x_{n - 1} = i. The rows of Gamma are * simplexes. @@ -31,8 +31,10 @@ namespace math { * because it only computes marginals on a state-by-state basis. * Each column of log_omegas is shifted by its maximum before it is * exponentiated, and the forward and backward passes are normalized at - * every step, so the result neither underflows nor overflows for long - * sequences or for log densities of large magnitude. + * every step, so long sequences and log densities of large magnitude do not + * underflow or overflow. The result is still NaN if the largest density in + * a column belongs to a state that the chain cannot reach, and the log + * densities of the reachable states are more than about 745 below it. * * @tparam T_omega type of the log likelihood matrix * @tparam T_Gamma type of the transition matrix @@ -55,6 +57,10 @@ inline Eigen::MatrixXd hmm_hidden_state_prob(const T_omega& log_omegas, int n_states = log_omegas.rows(); int n_transitions = log_omegas.cols() - 1; + ref_type_t rho_dbl = value_of(rho); + ref_type_t Gamma_dbl = value_of(Gamma); + hmm_check(log_omegas, Gamma_dbl, rho_dbl, "hmm_hidden_state_prob"); + // Subtract the maximum of each column before exponentiating, so that // omegas neither underflows nor overflows. The hidden state probabilities // are invariant to a positive factor on each column of omegas. @@ -64,9 +70,6 @@ inline Eigen::MatrixXd hmm_hidden_state_prob(const T_omega& log_omegas, = (log_omegas_dbl.rowwise() - log_omegas_dbl.colwise().maxCoeff()) .array() .exp(); - ref_type_t rho_dbl = value_of(rho); - ref_type_t Gamma_dbl = value_of(Gamma); - hmm_check(log_omegas, Gamma_dbl, rho_dbl, "hmm_hidden_state_prob"); Eigen::MatrixXd alphas(n_states, n_transitions + 1); alphas.col(0) = omegas.col(0).cwiseProduct(rho_dbl); diff --git a/stan/math/prim/prob/hmm_latent_rng.hpp b/stan/math/prim/prob/hmm_latent_rng.hpp index 24037e37df3..c9a13251115 100644 --- a/stan/math/prim/prob/hmm_latent_rng.hpp +++ b/stan/math/prim/prob/hmm_latent_rng.hpp @@ -18,8 +18,8 @@ namespace math { * In this setting, the hidden states are discrete * and takes values over the finite space {1, ..., K}. * log_omegas is a matrix of observational densities, where - * the (i, j)th entry corresponds to the density of the ith observation, y_i, - * given x_i = j. + * the (i, j)th entry corresponds to the density of the jth observation, y_j, + * given x_j = i. * The transition matrix Gamma is such that the (i, j)th entry is the * probability that x_n = j given x_{n - 1} = i. The rows of Gamma are * simplexes. @@ -47,11 +47,20 @@ inline std::vector hmm_latent_rng(const T_omega& log_omegas, int n_states = log_omegas.rows(); int n_transitions = log_omegas.cols() - 1; - Eigen::MatrixXd omegas = value_of(log_omegas).array().exp(); ref_type_t rho_dbl = value_of(rho); ref_type_t Gamma_dbl = value_of(Gamma); hmm_check(log_omegas, Gamma_dbl, rho_dbl, "hmm_latent_rng"); + // Subtract the maximum of each column before exponentiating, so that + // omegas neither underflows nor overflows. The sampling probabilities + // are invariant to a positive factor on each column of omegas. + ref_type_t log_omegas_dbl + = value_of(log_omegas); + Eigen::MatrixXd omegas + = (log_omegas_dbl.rowwise() - log_omegas_dbl.colwise().maxCoeff()) + .array() + .exp(); + Eigen::MatrixXd alphas(n_states, n_transitions + 1); alphas.col(0) = omegas.col(0).cwiseProduct(rho_dbl); alphas.col(0) /= alphas.col(0).maxCoeff(); diff --git a/stan/math/prim/prob/hmm_marginal.hpp b/stan/math/prim/prob/hmm_marginal.hpp index 8ad907f0da7..921c116a9f1 100644 --- a/stan/math/prim/prob/hmm_marginal.hpp +++ b/stan/math/prim/prob/hmm_marginal.hpp @@ -98,8 +98,15 @@ inline auto hmm_marginal(const T_omega& log_omegas, const T_Gamma& Gamma, eig_matrix_partial alphas(n_states, n_transitions + 1); eig_vector_partial alpha_log_norms(n_transitions + 1); + // Subtract the maximum of each column before exponentiating, so that + // omegas neither underflows nor overflows. The shifts are added back to + // the log density. The gradients are invariant to the shifts. + const auto& log_omegas_val = to_ref(value_of(log_omegas_ref)); + Eigen::Matrix log_omegas_max + = log_omegas_val.colwise().maxCoeff(); + eig_matrix_partial omegas + = (log_omegas_val.rowwise() - log_omegas_max).array().exp(); // compute the density using the forward algorithm. - eig_matrix_partial omegas = value_of(log_omegas_ref).array().exp(); T_partial_type norm_norm; auto log_marginal_density = hmm_marginal_val( omegas, Gamma_val, rho_val, alphas, alpha_log_norms, norm_norm); @@ -148,7 +155,7 @@ inline auto hmm_marginal(const T_omega& log_omegas, const T_Gamma& Gamma, edge<2>(ops_partials).partials_ = omegas.col(0) / exp(log_marginal_density); } - return ops_partials.build(log_marginal_density); + return ops_partials.build(log_marginal_density + log_omegas_max.sum()); } else { auto grad_corr_boundary = exp(kappa_log_norms(0) - norm_norm); eig_vector_partial C = Gamma_val * omegas.col(1).cwiseProduct(kappa[0]); @@ -176,7 +183,7 @@ inline auto hmm_marginal(const T_omega& log_omegas, const T_Gamma& Gamma, } } - return ops_partials.build(log_marginal_density); + return ops_partials.build(log_marginal_density + log_omegas_max.sum()); } } // namespace math diff --git a/test/unit/math/prim/prob/hmm_hidden_state_prob_test.cpp b/test/unit/math/prim/prob/hmm_hidden_state_prob_test.cpp index 122bd626f09..e53e01463d4 100644 --- a/test/unit/math/prim/prob/hmm_hidden_state_prob_test.cpp +++ b/test/unit/math/prim/prob/hmm_hidden_state_prob_test.cpp @@ -62,7 +62,7 @@ TEST_F(hmm_test, hidden_state_single_outcome) { Eigen::MatrixXd prob = hmm_hidden_state_prob(log_omegas_, Gamma, rho); - for (int i = 0; i < n_transitions_; i++) { + for (int i = 0; i <= n_transitions_; i++) { EXPECT_EQ(prob(0, i), 1); EXPECT_EQ(prob(1, i), 0); } @@ -80,7 +80,7 @@ TEST_F(hmm_test, hidden_state_identity_transition) { Eigen::MatrixXd prob = hmm_hidden_state_prob(log_omegas, Gamma, rho_); - for (int i = 0; i < n_transitions_; i++) { + for (int i = 0; i <= n_transitions_; i++) { EXPECT_FLOAT_EQ(prob(0, i), rho_(0)); EXPECT_FLOAT_EQ(prob(1, i), rho_(1)); } @@ -107,7 +107,7 @@ TEST(hmm_test_nonstandard, hidden_state_symmetry) { EXPECT_FLOAT_EQ(prob(0, 0), 0.3); EXPECT_FLOAT_EQ(prob(1, 0), 0.7); - for (int i = 1; i < n_transitions; i++) { + for (int i = 1; i <= n_transitions; i++) { EXPECT_FLOAT_EQ(prob(0, i), 0.5); EXPECT_FLOAT_EQ(prob(1, i), 0.5); } @@ -131,7 +131,7 @@ TEST(hmm_test_nonstandard, hidden_state_prob1) { EXPECT_FLOAT_EQ(prob(0, 0), 0.3); EXPECT_FLOAT_EQ(prob(1, 0), 0.7); - for (int i = 1; i < n_transitions; i++) { + for (int i = 1; i <= n_transitions; i++) { EXPECT_FLOAT_EQ(prob(0, i), 1); EXPECT_FLOAT_EQ(prob(1, i), 0); } @@ -210,3 +210,30 @@ TEST(hmm_test_nonstandard, hidden_state_long_sequence) { // magnitude of the log densities, hence the looser tolerance. EXPECT_LT((prob - prob_ref).cwiseAbs().maxCoeff(), 1e-10); } + +TEST_F(hmm_test, hidden_state_neg_inf_log_omega) { + // A log density of -inf rules out that state for that observation. + using stan::math::hmm_hidden_state_prob; + Eigen::MatrixXd log_omegas = log_omegas_; + log_omegas(0, 4) = -std::numeric_limits::infinity(); + + Eigen::MatrixXd prob = hmm_hidden_state_prob(log_omegas, Gamma_, rho_); + Eigen::MatrixXd prob_ref + = hmm_hidden_state_prob_log_space(log_omegas, Gamma_, rho_); + + EXPECT_EQ(prob(0, 4), 0); + EXPECT_TRUE(prob.allFinite()); + EXPECT_LT((prob - prob_ref).cwiseAbs().maxCoeff(), 1e-12); +} + +TEST(hmm_test_nonstandard, hidden_state_no_states_throws) { + // The arguments must be checked before log_omegas is shifted by the + // maximum of each column, which is undefined for an empty column. + using stan::math::hmm_hidden_state_prob; + Eigen::MatrixXd Gamma(2, 2); + Gamma << 0.5, 0.5, 0.5, 0.5; + Eigen::VectorXd rho(0); + Eigen::MatrixXd log_omegas(0, 3); + EXPECT_THROW(hmm_hidden_state_prob(log_omegas, Gamma, rho), + std::invalid_argument); +} diff --git a/test/unit/math/prim/prob/hmm_latent_rng_test.cpp b/test/unit/math/prim/prob/hmm_latent_rng_test.cpp index 0a5606e9abd..44e31d25f42 100644 --- a/test/unit/math/prim/prob/hmm_latent_rng_test.cpp +++ b/test/unit/math/prim/prob/hmm_latent_rng_test.cpp @@ -133,3 +133,24 @@ TEST(hmm_rng_test, chiSquareGoodnessFitTest_symmetric) { // double threshold = 0.1; // CHECK -- what is an appropriate threshold? // EXPECT_TRUE(p_value > threshold); } + +TEST_F(hmm_test, latent_rng_shifted_log_omegas) { + // Adding a constant to a column of log_omegas does not change the + // posterior of the hidden states, so with the same seed the draws match + // the draws for the unshifted log_omegas. Shifts of -1000 and 1000 + // underflow or overflow if log_omegas is exponentiated directly (#2677). + using stan::math::hmm_latent_rng; + for (double shift : {-1000.0, 1000.0}) { + Eigen::MatrixXd log_omegas = log_omegas_; + for (int n = 0; n < log_omegas.cols(); ++n) { + log_omegas.col(n).array() += shift + 10 * n; + } + boost::random::mt19937 rng(2677); + boost::random::mt19937 rng_ref(2677); + for (int i = 0; i < 100; ++i) { + EXPECT_EQ(hmm_latent_rng(log_omegas, Gamma_, rho_, rng), + hmm_latent_rng(log_omegas_, Gamma_, rho_, rng_ref)) + << "shift " << shift << ", draw " << i; + } + } +} diff --git a/test/unit/math/prim/prob/hmm_marginal_test.cpp b/test/unit/math/prim/prob/hmm_marginal_test.cpp index 3f9b653eb31..6c56d8b5eac 100644 --- a/test/unit/math/prim/prob/hmm_marginal_test.cpp +++ b/test/unit/math/prim/prob/hmm_marginal_test.cpp @@ -140,3 +140,41 @@ TEST(hmm_marginal, exceptions) { " all arguments must be scalars or multidimensional values of" " the same shape.") } + +TEST_F(hmm_test, shifted_log_omegas) { + // Adding a constant to a column of log_omegas adds it to the log marginal + // density and leaves the gradients unchanged. Shifts of -1000 and 1000 + // underflow or overflow if log_omegas is exponentiated directly (#2677). + using stan::math::hmm_marginal; + using stan::math::var; + auto gradients = [&](const Eigen::MatrixXd& log_omegas) { + Eigen::Matrix log_omegas_v + = log_omegas; + Eigen::Matrix Gamma_v = Gamma_; + Eigen::Matrix rho_v = rho_; + hmm_marginal(log_omegas_v, Gamma_v, rho_v).grad(); + std::vector adjs{log_omegas_v.adj(), Gamma_v.adj(), + rho_v.adj()}; + stan::math::recover_memory(); + return adjs; + }; + std::vector adjs_ref = gradients(log_omegas_); + + for (double shift : {-1000.0, 1000.0}) { + Eigen::MatrixXd log_omegas = log_omegas_; + double total_shift = 0; + for (int n = 0; n < log_omegas.cols(); ++n) { + log_omegas.col(n).array() += shift + 10 * n; + total_shift += shift + 10 * n; + } + EXPECT_FLOAT_EQ(hmm_marginal(log_omegas, Gamma_, rho_) - total_shift, + -18.37417) + << "shift " << shift; + + std::vector adjs = gradients(log_omegas); + for (size_t i = 0; i < adjs.size(); ++i) { + EXPECT_LT((adjs[i] - adjs_ref[i]).cwiseAbs().maxCoeff(), 1e-10) + << "shift " << shift << ", argument " << i; + } + } +} From 2a468d6f39399665537f00605b4f34a6155b20d8 Mon Sep 17 00:00:00 2001 From: Steve Bronder Date: Wed, 28 Oct 2026 14:08:37 -0400 Subject: [PATCH 3/3] Reject log_omegas with no columns in hmm_check hmm_check did not check that log_omegas has at least one column. With no columns, n_transitions is -1 and hmm_hidden_state_prob, hmm_latent_rng and hmm_marginal index alphas.col(0) out of bounds. Throw std::invalid_argument instead. --- stan/math/prim/err/hmm_check.hpp | 6 ++++-- test/unit/math/prim/err/hmm_check_test.cpp | 7 +++++++ .../unit/math/prim/prob/hmm_hidden_state_prob_test.cpp | 9 +++++++++ test/unit/math/prim/prob/hmm_latent_rng_test.cpp | 10 ++++++++++ test/unit/math/prim/prob/hmm_marginal_test.cpp | 6 ++++++ 5 files changed, 36 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/err/hmm_check.hpp b/stan/math/prim/err/hmm_check.hpp index 32b39ed325a..bc01dc71e9b 100644 --- a/stan/math/prim/err/hmm_check.hpp +++ b/stan/math/prim/err/hmm_check.hpp @@ -19,8 +19,9 @@ namespace math { * @param[in] Gamma transition density between hidden states. * @param[in] rho initial state * @param[in] function the name of the function using the arguments. - * @throw `std::invalid_argument` if Gamma is not square - * or if the size of rho is not the number of rows of log_omegas. + * @throw `std::invalid_argument` if Gamma is not square, if log_omegas + * has size 0, or if the size of rho is not the number of rows of + * log_omegas. * @throw `std::domain_error` if rho is not a simplex or the rows * of Gamma are not a simplex. */ @@ -34,6 +35,7 @@ inline void hmm_check(const T_omega& log_omegas, const T_Gamma& Gamma, check_consistent_size(function, "rho", rho, n_states); check_square(function, "Gamma", Gamma); check_nonzero_size(function, "Gamma", Gamma); + check_nonzero_size(function, "log_omegas", log_omegas); check_multiplicable(function, "Gamma", Gamma, "log_omegas", log_omegas); check_simplex(function, "rho", rho); diff --git a/test/unit/math/prim/err/hmm_check_test.cpp b/test/unit/math/prim/err/hmm_check_test.cpp index f2862746199..45ee3d371a8 100644 --- a/test/unit/math/prim/err/hmm_check_test.cpp +++ b/test/unit/math/prim/err/hmm_check_test.cpp @@ -53,6 +53,13 @@ TEST(err, hmm_check) { "hmm_marginal_lpdf: Columns of Gamma (3)" " and Rows of log_omegas (2) must match in size") + // log_omegas has no columns, so there are no observations. + MatrixXd log_omegas_empty(n_states, 0); + EXPECT_THROW_MSG( + hmm_check(log_omegas_empty, Gamma, rho, "hmm_marginal_lpdf"), + std::invalid_argument, + "hmm_marginal_lpdf: log_omegas has size 0, but must have a non-zero size") + // rho is not a simplex. VectorXd rho_bad = rho; rho_bad(0) = rho(0) + 1; diff --git a/test/unit/math/prim/prob/hmm_hidden_state_prob_test.cpp b/test/unit/math/prim/prob/hmm_hidden_state_prob_test.cpp index e53e01463d4..2f7fbd0c217 100644 --- a/test/unit/math/prim/prob/hmm_hidden_state_prob_test.cpp +++ b/test/unit/math/prim/prob/hmm_hidden_state_prob_test.cpp @@ -237,3 +237,12 @@ TEST(hmm_test_nonstandard, hidden_state_no_states_throws) { EXPECT_THROW(hmm_hidden_state_prob(log_omegas, Gamma, rho), std::invalid_argument); } + +TEST_F(hmm_test, hidden_state_no_observations_throws) { + using stan::math::hmm_hidden_state_prob; + Eigen::MatrixXd log_omegas(n_states_, 0); + EXPECT_THROW_MSG(hmm_hidden_state_prob(log_omegas, Gamma_, rho_), + std::invalid_argument, + "hmm_hidden_state_prob: log_omegas has size 0, but must " + "have a non-zero size"); +} diff --git a/test/unit/math/prim/prob/hmm_latent_rng_test.cpp b/test/unit/math/prim/prob/hmm_latent_rng_test.cpp index 44e31d25f42..8dfdfb08d92 100644 --- a/test/unit/math/prim/prob/hmm_latent_rng_test.cpp +++ b/test/unit/math/prim/prob/hmm_latent_rng_test.cpp @@ -154,3 +154,13 @@ TEST_F(hmm_test, latent_rng_shifted_log_omegas) { } } } + +TEST_F(hmm_test, latent_rng_no_observations_throws) { + using stan::math::hmm_latent_rng; + Eigen::MatrixXd log_omegas(n_states_, 0); + boost::random::mt19937 rng; + EXPECT_THROW_MSG(hmm_latent_rng(log_omegas, Gamma_, rho_, rng), + std::invalid_argument, + "hmm_latent_rng: log_omegas has size 0, but must have a " + "non-zero size"); +} diff --git a/test/unit/math/prim/prob/hmm_marginal_test.cpp b/test/unit/math/prim/prob/hmm_marginal_test.cpp index 6c56d8b5eac..04418a5ad4d 100644 --- a/test/unit/math/prim/prob/hmm_marginal_test.cpp +++ b/test/unit/math/prim/prob/hmm_marginal_test.cpp @@ -123,6 +123,12 @@ TEST(hmm_marginal, exceptions) { "hmm_marginal: Columns of Gamma (3)" " and Rows of log_omegas (2) must match in size") + // log_omegas has no columns, so there are no observations. + MatrixXd log_omegas_empty(n_states, 0); + EXPECT_THROW_MSG( + hmm_marginal(log_omegas_empty, Gamma, rho), std::invalid_argument, + "hmm_marginal: log_omegas has size 0, but must have a non-zero size") + // rho is not a simplex. VectorXd rho_bad = rho; rho_bad(0) = rho(0) + 1;