Automatic Differentiation
 
Loading...
Searching...
No Matches
hmm_hidden_state_prob.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_PROB_HMM_HIDDEN_STATE_PROB_HPP
2#define STAN_MATH_PRIM_PROB_HMM_HIDDEN_STATE_PROB_HPP
3
8#include <boost/random.hpp>
9#include <vector>
10
11namespace stan {
12namespace math {
13
45template <typename T_omega, typename T_Gamma, typename T_rho,
46 require_all_eigen_t<T_omega, T_Gamma>* = nullptr,
47 require_eigen_col_vector_t<T_rho>* = nullptr>
48inline Eigen::MatrixXd hmm_hidden_state_prob(const T_omega& log_omegas,
49 const T_Gamma& Gamma,
50 const T_rho& rho) {
51 int n_states = log_omegas.rows();
52 int n_transitions = log_omegas.cols() - 1;
53
54 Eigen::MatrixXd omegas = value_of(log_omegas).array().exp();
55 ref_type_t<decltype(value_of(rho))> rho_dbl = value_of(rho);
56 ref_type_t<decltype(value_of(Gamma))> Gamma_dbl = value_of(Gamma);
57 hmm_check(log_omegas, Gamma_dbl, rho_dbl, "hmm_hidden_state_prob");
58
59 Eigen::MatrixXd alphas(n_states, n_transitions + 1);
60 alphas.col(0) = omegas.col(0).cwiseProduct(rho_dbl);
61 alphas.col(0) /= alphas.col(0).maxCoeff();
62
63 for (int n = 0; n < n_transitions; ++n)
64 alphas.col(n + 1)
65 = omegas.col(n + 1).cwiseProduct(Gamma_dbl.transpose() * alphas.col(n));
66
67 // Backward pass with running normalization
68 Eigen::VectorXd beta = Eigen::VectorXd::Ones(n_states);
69
70 alphas.col(n_transitions) /= alphas.col(n_transitions).sum();
71
72 for (int n = n_transitions; n-- > 0;) {
73 beta = Gamma_dbl * omegas.col(n + 1).cwiseProduct(beta);
74 beta /= beta.maxCoeff();
75
76 // Reuse alphas to store probabilities
77 alphas.col(n) = alphas.col(n).cwiseProduct(beta);
78 alphas.col(n) /= alphas.col(n).sum();
79 }
80
81 return alphas;
82}
83
84} // namespace math
85} // namespace stan
86#endif
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition value_of.hpp:18
void hmm_check(const T_omega &log_omegas, const T_Gamma &Gamma, const T_rho &rho, const char *function)
Check arguments for hidden Markov model functions with a discrete latent state (lpdf,...
Definition hmm_check.hpp:30
Eigen::MatrixXd hmm_hidden_state_prob(const T_omega &log_omegas, const T_Gamma &Gamma, const T_rho &rho)
For a hidden Markov model with observation y, hidden state x, and parameters theta,...
fvar< T > beta(const fvar< T > &x1, const fvar< T > &x2)
Return fvar with the beta function applied to the specified arguments and its gradient.
Definition beta.hpp:51
typename ref_type_if< true, T >::type ref_type_t
Definition ref_type.hpp:55
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...