Stan Math Library  2.20.0
reverse mode automatic differentiation
gaussian_dlm_obs_rng.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_MAT_PROB_GAUSSIAN_DLM_OBS_RNG_HPP
2 #define STAN_MATH_PRIM_MAT_PROB_GAUSSIAN_DLM_OBS_RNG_HPP
3 
13 #include <boost/random/normal_distribution.hpp>
14 #include <boost/random/variate_generator.hpp>
15 #include <vector>
16 
17 namespace stan {
18 namespace math {
19 namespace internal {
20 
34 template <class RNG>
35 inline Eigen::VectorXd multi_normal_semidefinite_rng(
36  const Eigen::VectorXd &mu, const Eigen::LDLT<Eigen::MatrixXd> &S_ldlt,
37  RNG &rng) {
38  using boost::normal_distribution;
39  using boost::variate_generator;
40 
41  variate_generator<RNG &, normal_distribution<>> std_normal_rng(
42  rng, normal_distribution<>(0, 1));
43 
44  Eigen::VectorXd stddev = S_ldlt.vectorD().array().sqrt().matrix();
45  size_t M = S_ldlt.vectorD().size();
46  Eigen::VectorXd z(M);
47  for (int i = 0; i < M; i++)
48  z(i) = stddev(i) * std_normal_rng();
49 
50  Eigen::VectorXd Y
51  = mu + (S_ldlt.transpositionsP().transpose() * (S_ldlt.matrixL() * z));
52  // The inner paranthesis matter, transpositionsP() gives a
53  // permutation matrix from pivoting and matrixL() gives a lower
54  // triangular matrix. The types cannot be directly multiplied.
55 
56  return Y;
57 }
58 
72 template <class RNG>
73 inline Eigen::VectorXd multi_normal_definite_rng(
74  const Eigen::VectorXd &mu, const Eigen::LLT<Eigen::MatrixXd> &S_llt,
75  RNG &rng) {
76  using boost::normal_distribution;
77  using boost::variate_generator;
78 
79  variate_generator<RNG &, normal_distribution<>> std_normal_rng(
80  rng, normal_distribution<>(0, 1));
81 
82  size_t M = S_llt.matrixL().rows();
83  Eigen::VectorXd z(M);
84  for (int i = 0; i < M; i++)
85  z(i) = std_normal_rng();
86 
87  Eigen::VectorXd Y = mu + S_llt.matrixL() * z;
88 
89  return Y;
90 }
91 
92 } // namespace internal
93 
125 template <class RNG>
126 inline Eigen::MatrixXd gaussian_dlm_obs_rng(const Eigen::MatrixXd &F,
127  const Eigen::MatrixXd &G,
128  const Eigen::MatrixXd &V,
129  const Eigen::MatrixXd &W,
130  const Eigen::VectorXd &m0,
131  const Eigen::MatrixXd &C0,
132  const int T, RNG &rng) {
133  static const char *function = "gaussian_dlm_obs_rng";
134 
135  int r = F.cols(); // number of variables
136  int n = G.rows(); // number of states
137 
138  check_size_match(function, "rows of F", F.rows(), "rows of G", G.rows());
139  check_finite(function, "F", F);
140  check_square(function, "G", G);
141  check_finite(function, "G", G);
142  check_size_match(function, "rows of V", V.rows(), "cols of F", F.cols());
143  check_finite(function, "V", V);
144  check_positive(function, "V rows", V.rows());
145  check_symmetric(function, "V", V);
146  check_size_match(function, "rows of W", W.rows(), "rows of G", G.rows());
147  check_finite(function, "W", W);
148  check_positive(function, "W rows", W.rows());
149  check_symmetric(function, "W", W);
150  check_size_match(function, "rows of W", W.rows(), "rows of G", G.rows());
151  check_size_match(function, "size of m0", m0.size(), "rows of G", G.rows());
152  check_finite(function, "m0", m0);
153  check_size_match(function, "rows of C0", C0.rows(), "rows of G", G.rows());
154  check_finite(function, "C0", C0);
155  check_positive(function, "C0 rows", C0.rows());
156  check_symmetric(function, "C0", C0);
157  check_positive(function, "T", T);
158 
159  Eigen::LDLT<Eigen::MatrixXd> V_ldlt = V.ldlt();
160  check_pos_semidefinite(function, "V", V_ldlt);
161  Eigen::LDLT<Eigen::MatrixXd> W_ldlt = W.ldlt();
162  check_pos_semidefinite(function, "W", W_ldlt);
163  // Require strictly positive definite C0 for consistency with lpdf.
164  Eigen::LLT<Eigen::MatrixXd> C0_llt = C0.llt();
165  check_pos_definite(function, "C0", C0_llt);
166 
167  Eigen::MatrixXd y(r, T);
168  Eigen::VectorXd theta_t
169  = internal::multi_normal_definite_rng(m0, C0_llt, rng);
170  for (int t = 0; t < T; ++t) {
172  Eigen::VectorXd(G * theta_t), W_ldlt, rng);
174  Eigen::VectorXd(F.transpose() * theta_t), V_ldlt, rng);
175  }
176  return y;
177 }
178 
179 template <class RNG>
180 inline Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>
181 gaussian_dlm_obs_rng(const Eigen::MatrixXd &F, const Eigen::MatrixXd &G,
182  const Eigen::VectorXd &V, const Eigen::MatrixXd &W,
183  const Eigen::VectorXd &m0, const Eigen::MatrixXd &C0,
184  const int T, RNG &rng) {
185  return gaussian_dlm_obs_rng(F, G, Eigen::MatrixXd(V.asDiagonal()), W, m0, C0,
186  T, rng);
187 }
188 
189 } // namespace math
190 } // namespace stan
191 #endif
void check_finite(const char *function, const char *name, const T_y &y)
Check if y is finite.
Eigen::MatrixXd gaussian_dlm_obs_rng(const Eigen::MatrixXd &F, const Eigen::MatrixXd &G, const Eigen::MatrixXd &V, const Eigen::MatrixXd &W, const Eigen::VectorXd &m0, const Eigen::MatrixXd &C0, const int T, RNG &rng)
Simulate random draw from Gaussian dynamic linear model (GDLM).
void check_square(const char *function, const char *name, const matrix_cl &y)
Check if the matrix_cl is square.
Eigen::VectorXd multi_normal_definite_rng(const Eigen::VectorXd &mu, const Eigen::LLT< Eigen::MatrixXd > &S_llt, RNG &rng)
Return a multivariate normal random variate with the given location and covariance using the specifie...
void check_size_match(const char *function, const char *name_i, T_size1 i, const char *name_j, T_size2 j)
Check if the provided sizes match.
void check_pos_semidefinite(const char *function, const char *name, const Eigen::Matrix< T_y, Eigen::Dynamic, Eigen::Dynamic > &y)
Check if the specified matrix is positive definite.
Eigen::VectorXd multi_normal_semidefinite_rng(const Eigen::VectorXd &mu, const Eigen::LDLT< Eigen::MatrixXd > &S_ldlt, RNG &rng)
Return a multivariate normal random variate with the given location and covariance using the specifie...
void check_positive(const char *function, const char *name, const T_y &y)
Check if y is positive.
void check_pos_definite(const char *function, const char *name, const Eigen::Matrix< T_y, -1, -1 > &y)
Check if the specified square, symmetric matrix is positive definite.
void check_symmetric(const char *function, const char *name, const matrix_cl &y)
Check if the matrix_cl is symmetric.

     [ Stan Home Page ] © 2011–2018, Stan Development Team.