hmm_hidden_state_prob.hpp

Go to the documentation of this file.

#ifndef STAN_MATH_PRIM_PROB_HMM_HIDDEN_STATE_PROB_HPP
#define STAN_MATH_PRIM_PROB_HMM_HIDDEN_STATE_PROB_HPP
 
#include <stan/math/prim/core.hpp>
#include <stan/math/prim/meta.hpp>
#include <stan/math/prim/err/hmm_check.hpp>
#include <stan/math/prim/fun/Eigen.hpp>
#include <boost/random.hpp>
#include <vector>
 
namespace stan {
namespace math {
 
template <typename T_omega, typename T_Gamma, typename T_rho,
          require_all_eigen_t<T_omega, T_Gamma>* = nullptr,
          require_eigen_col_vector_t<T_rho>* = nullptr>
inline Eigen::MatrixXd hmm_hidden_state_prob(const T_omega& log_omegas,
                                             const T_Gamma& Gamma,
                                             const T_rho& rho) {
  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<decltype(value_of(rho))> rho_dbl = value_of(rho);
  ref_type_t<decltype(value_of(Gamma))> 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);
  alphas.col(0) /= alphas.col(0).maxCoeff();
 
  for (int n = 0; n < n_transitions; ++n)
    alphas.col(n + 1)
        = omegas.col(n + 1).cwiseProduct(Gamma_dbl.transpose() * alphas.col(n));
 
  // Backward pass with running normalization
  Eigen::VectorXd beta = Eigen::VectorXd::Ones(n_states);
 
  alphas.col(n_transitions) /= alphas.col(n_transitions).sum();
 
  for (int n = n_transitions; n-- > 0;) {
    beta = Gamma_dbl * omegas.col(n + 1).cwiseProduct(beta);
    beta /= beta.maxCoeff();
 
    // Reuse alphas to store probabilities
    alphas.col(n) = alphas.col(n).cwiseProduct(beta);
    alphas.col(n) /= alphas.col(n).sum();
  }
 
  return alphas;
}
 
}  // namespace math
}  // namespace stan
#endif