Automatic Differentiation
 
Loading...
Searching...
No Matches
gaussian_dlm_obs_rng.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_PROB_GAUSSIAN_DLM_OBS_RNG_HPP
2#define STAN_MATH_PRIM_PROB_GAUSSIAN_DLM_OBS_RNG_HPP
3
7#include <boost/random/normal_distribution.hpp>
8#include <boost/random/variate_generator.hpp>
9#include <vector>
10
11namespace stan {
12namespace math {
13namespace internal {
14
28template <class RNG>
29inline Eigen::VectorXd multi_normal_semidefinite_rng(
30 const Eigen::VectorXd &mu, const Eigen::LDLT<Eigen::MatrixXd> &S_ldlt,
31 RNG &rng) {
32 using boost::normal_distribution;
33 using boost::variate_generator;
34
35 variate_generator<RNG &, normal_distribution<>> std_normal_rng(
36 rng, normal_distribution<>(0, 1));
37
38 Eigen::VectorXd stddev = S_ldlt.vectorD().array().sqrt().matrix();
39 size_t M = S_ldlt.vectorD().size();
40 Eigen::VectorXd z(M);
41 for (int i = 0; i < M; i++) {
42 z(i) = stddev(i) * std_normal_rng();
43 }
44
45 Eigen::VectorXd Y
46 = mu + (S_ldlt.transpositionsP().transpose() * (S_ldlt.matrixL() * z));
47 // The inner paranthesis matter, transpositionsP() gives a
48 // permutation matrix from pivoting and matrixL() gives a lower
49 // triangular matrix. The types cannot be directly multiplied.
50
51 return Y;
52}
53
54} // namespace internal
55
87template <class RNG>
88inline Eigen::MatrixXd gaussian_dlm_obs_rng(const Eigen::MatrixXd &F,
89 const Eigen::MatrixXd &G,
90 const Eigen::MatrixXd &V,
91 const Eigen::MatrixXd &W,
92 const Eigen::VectorXd &m0,
93 const Eigen::MatrixXd &C0,
94 const int T, RNG &rng) {
95 static constexpr const char *function = "gaussian_dlm_obs_rng";
96
97 int r = F.cols(); // number of variables
98 int n = G.rows(); // number of states
99
100 check_size_match(function, "rows of F", F.rows(), "rows of G", n);
101 check_finite(function, "F", F);
102 check_square(function, "G", G);
103 check_finite(function, "G", G);
104 check_size_match(function, "rows of V", V.rows(), "cols of F", r);
105 check_finite(function, "V", V);
106 check_positive(function, "V rows", V.rows());
107 check_symmetric(function, "V", V);
108 check_size_match(function, "rows of W", W.rows(), "rows of G", n);
109 check_finite(function, "W", W);
110 check_positive(function, "W rows", W.rows());
111 check_symmetric(function, "W", W);
112 check_size_match(function, "rows of W", W.rows(), "rows of G", n);
113 check_size_match(function, "size of m0", m0.size(), "rows of G", n);
114 check_finite(function, "m0", m0);
115 check_size_match(function, "rows of C0", C0.rows(), "rows of G", n);
116 check_finite(function, "C0", C0);
117 check_positive(function, "C0 rows", C0.rows());
118 check_symmetric(function, "C0", C0);
119 check_positive(function, "T", T);
120
121 Eigen::LDLT<Eigen::MatrixXd> V_ldlt = V.ldlt();
122 check_pos_semidefinite(function, "V", V_ldlt);
123 Eigen::LDLT<Eigen::MatrixXd> W_ldlt = W.ldlt();
124 check_pos_semidefinite(function, "W", W_ldlt);
125 Eigen::LDLT<Eigen::MatrixXd> C0_ldlt = C0.ldlt();
126 check_pos_semidefinite(function, "C0", C0_ldlt);
127
128 Eigen::MatrixXd y(r, T);
129 Eigen::VectorXd theta_t
130 = internal::multi_normal_semidefinite_rng(m0, C0_ldlt, rng);
131 for (int t = 0; t < T; ++t) {
133 Eigen::VectorXd(G * theta_t), W_ldlt, rng);
135 Eigen::VectorXd(F.transpose() * theta_t), V_ldlt, rng);
136 }
137 return y;
138}
139
140template <class RNG>
141inline Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>
142gaussian_dlm_obs_rng(const Eigen::MatrixXd &F, const Eigen::MatrixXd &G,
143 const Eigen::VectorXd &V, const Eigen::MatrixXd &W,
144 const Eigen::VectorXd &m0, const Eigen::MatrixXd &C0,
145 const int T, RNG &rng) {
146 return gaussian_dlm_obs_rng(F, G, Eigen::MatrixXd(V.asDiagonal()), W, m0, C0,
147 T, rng);
148}
149
150} // namespace math
151} // namespace stan
152#endif
void check_symmetric(const char *function, const char *name, const matrix_cl< T > &y)
Check if the matrix_cl is symmetric.
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).
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...
double std_normal_rng(RNG &rng)
Return a standard Normal random variate using the specified random number generator.
void check_square(const char *function, const char *name, const T_y &y)
Check if the specified matrix is square.
void check_pos_semidefinite(const char *function, const char *name, const EigMat &y)
Check if the specified matrix is positive definite.
void check_finite(const char *function, const char *name, const T_y &y)
Return true if all values in y are finite.
void check_positive(const char *function, const char *name, const T_y &y)
Check if y is positive.
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.
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...