Automatic Differentiation
 
Loading...
Searching...
No Matches
hmm_latent_rng.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_PROB_HMM_LATENT_RNG_HPP
2#define STAN_MATH_PRIM_PROB_HMM_LATENT_RNG_HPP
3
8#include <boost/random.hpp>
9#include <vector>
10
11namespace stan {
12namespace math {
13
41template <typename T_omega, typename T_Gamma, typename T_rho, class RNG,
42 require_all_eigen_t<T_omega, T_Gamma>* = nullptr,
43 require_eigen_col_vector_t<T_rho>* = nullptr>
44inline std::vector<int> hmm_latent_rng(const T_omega& log_omegas,
45 const T_Gamma& Gamma, const T_rho& rho,
46 RNG& rng) {
47 int n_states = log_omegas.rows();
48 int n_transitions = log_omegas.cols() - 1;
49
50 Eigen::MatrixXd omegas = value_of(log_omegas).array().exp();
51 ref_type_t<decltype(value_of(rho))> rho_dbl = value_of(rho);
52 ref_type_t<decltype(value_of(Gamma))> Gamma_dbl = value_of(Gamma);
53 hmm_check(log_omegas, Gamma_dbl, rho_dbl, "hmm_latent_rng");
54
55 Eigen::MatrixXd alphas(n_states, n_transitions + 1);
56 alphas.col(0) = omegas.col(0).cwiseProduct(rho_dbl);
57 alphas.col(0) /= alphas.col(0).maxCoeff();
58
59 for (int n = 0; n < n_transitions; ++n) {
60 alphas.col(n + 1)
61 = omegas.col(n + 1).cwiseProduct(Gamma_dbl.transpose() * alphas.col(n));
62 alphas.col(n + 1) /= alphas.col(n + 1).maxCoeff();
63 }
64
65 Eigen::VectorXd beta = Eigen::VectorXd::Ones(n_states);
66
67 // sample the last hidden state
68 std::vector<int> hidden_states(n_transitions + 1);
69 std::vector<double> probs(n_states);
70 Eigen::Map<Eigen::VectorXd> probs_vec(probs.data(), n_states);
71 probs_vec = alphas.col(n_transitions) / alphas.col(n_transitions).sum();
72 boost::random::discrete_distribution<> cat_hidden(probs);
73 hidden_states[n_transitions] = cat_hidden(rng) + stan::error_index::value;
74
75 for (int n = n_transitions; n-- > 0;) {
76 // Sample the nth hidden state conditional on (n + 1)st hidden state.
77 // Subtract error_index in order to use C++ index.
78 int last_hs = hidden_states[n + 1] - stan::error_index::value;
79
80 probs_vec = alphas.col(n).cwiseProduct(Gamma_dbl.col(last_hs))
81 * beta(last_hs) * omegas(last_hs, n + 1);
82
83 probs_vec /= probs_vec.sum();
84
85 // discrete_distribution produces samples in [0, K), so
86 // we need to add 1 to generate over [1, K).
87 boost::random::discrete_distribution<> cat_hidden(probs);
88 hidden_states[n] = cat_hidden(rng) + stan::error_index::value;
89
90 // update backwards state
91 beta = Gamma_dbl * (omegas.col(n + 1).cwiseProduct(beta));
92 beta /= beta.maxCoeff();
93 }
94
95 return hidden_states;
96}
97
98} // namespace math
99} // namespace stan
100#endif
std::vector< int > hmm_latent_rng(const T_omega &log_omegas, const T_Gamma &Gamma, const T_rho &rho, RNG &rng)
For a hidden Markov model with observation y, hidden state x, and parameters theta,...
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
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 ...