Stan Math Library  2.20.0
reverse mode automatic differentiation
matrix_normal_prec_rng.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_MAT_PROB_MATRIX_NORMAL_PREC_RNG_HPP
2 #define STAN_MATH_PRIM_MAT_PROB_MATRIX_NORMAL_PREC_RNG_HPP
3 
11 #include <boost/random/normal_distribution.hpp>
12 #include <boost/random/variate_generator.hpp>
13 
14 namespace stan {
15 namespace math {
16 
35 template <class RNG>
36 inline Eigen::MatrixXd matrix_normal_prec_rng(const Eigen::MatrixXd &Mu,
37  const Eigen::MatrixXd &Sigma,
38  const Eigen::MatrixXd &D,
39  RNG &rng) {
40  using boost::normal_distribution;
41  using boost::variate_generator;
42 
43  static const char *function = "matrix_normal_prec_rng";
44 
45  check_positive(function, "Sigma rows", Sigma.rows());
46  check_finite(function, "Sigma", Sigma);
47  check_symmetric(function, "Sigma", Sigma);
48 
49  check_positive(function, "D rows", D.rows());
50  check_finite(function, "D", D);
51  check_symmetric(function, "D", D);
52 
53  check_size_match(function, "Rows of location parameter", Mu.rows(),
54  "Rows of Sigma", Sigma.rows());
55  check_size_match(function, "Columns of location parameter", Mu.cols(),
56  "Rows of D", D.rows());
57 
58  check_finite(function, "Location parameter", Mu);
59 
60  Eigen::LDLT<Eigen::MatrixXd> Sigma_ldlt(Sigma);
61  // Sigma = PS^T LS DS LS^T PS
62  // PS a permutation matrix.
63  // LS lower triangular with unit diagonal.
64  // DS diagonal.
65  Eigen::LDLT<Eigen::MatrixXd> D_ldlt(D);
66  // D = PD^T LD DD LD^T PD
67 
68  check_pos_semidefinite(function, "Sigma", Sigma_ldlt);
69  check_pos_semidefinite(function, "D", D_ldlt);
70 
71  // If
72  // C ~ N[0, I, I]
73  // Then
74  // A C B ~ N[0, A A^T, B^T B]
75  // So to get
76  // Y - Mu ~ N[0, Sigma^(-1), D^(-1)]
77  // We need to do
78  // Y - Mu = Q^T^(-1) C R^(-1)
79  // Where Q^T^(-1) and R^(-1) are such that
80  // Q^(-1) Q^(-1)^T = Sigma^(-1)
81  // R^(-1)^T R^(-1) = D^(-1)
82  // We choose:
83  // Q^(-1)^T = PS^T LS^T^(-1) sqrt[DS]^(-1)
84  // R^(-1) = sqrt[DD]^(-1) LD^(-1) PD
85  // And therefore
86  // Y - Mu = (PS^T LS^T^(-1) sqrt[DS]^(-1)) C (sqrt[DD]^(-1) LD^(-1) PD)
87 
88  int m = Sigma.rows();
89  int n = D.rows();
90 
91  variate_generator<RNG &, normal_distribution<>> std_normal_rng(
92  rng, normal_distribution<>(0, 1));
93 
94  // X = sqrt[DS]^(-1) C sqrt[DD]^(-1)
95  // X ~ N[0, DS, DD]
96  Eigen::MatrixXd X(m, n);
97  Eigen::VectorXd row_stddev
98  = Sigma_ldlt.vectorD().array().inverse().sqrt().matrix();
99  Eigen::VectorXd col_stddev
100  = D_ldlt.vectorD().array().inverse().sqrt().matrix();
101  for (int row = 0; row < m; ++row) {
102  for (int col = 0; col < n; ++col) {
103  double stddev = row_stddev(row) * col_stddev(col);
104  // C(row, col) = std_normal_rng();
105  X(row, col) = stddev * std_normal_rng();
106  }
107  }
108 
109  // Y - Mu = PS^T (LS^T^(-1) X LD^(-1)) PD
110  // Y' = LS^T^(-1) X LD^(-1)
111  // Y' = LS^T.solve(X) LD^(-1)
112  // Y' = (LD^(-1)^T (LS^T.solve(X))^T)^T
113  // Y' = (LD^T.solve((LS^T.solve(X))^T))^T
114  // Y = Mu + PS^T Y' PD
115  Eigen::MatrixXd Y = Mu
116  + (Sigma_ldlt.transpositionsP().transpose()
117  * (D_ldlt.matrixU().solve(
118  (Sigma_ldlt.matrixU().solve(X)).transpose()))
119  .transpose()
120  * D_ldlt.transpositionsP());
121 
122  return Y;
123 }
124 } // namespace math
125 } // namespace stan
126 #endif
void check_finite(const char *function, const char *name, const T_y &y)
Check if y is finite.
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::Matrix< T, 1, Eigen::Dynamic > row(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m, size_t i)
Return the specified row of the specified matrix, using start-at-1 indexing.
Definition: row.hpp:24
Eigen::MatrixXd matrix_normal_prec_rng(const Eigen::MatrixXd &Mu, const Eigen::MatrixXd &Sigma, const Eigen::MatrixXd &D, RNG &rng)
Sample from the the matrix normal distribution for the given Mu, Sigma and D where Sigma and D are gi...
Eigen::Matrix< T, Eigen::Dynamic, 1 > col(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m, size_t j)
Return the specified column of the specified matrix using start-at-1 indexing.
Definition: col.hpp:23
matrix_cl transpose(const matrix_cl &src)
Takes the transpose of the matrix on the OpenCL device.
Definition: transpose.hpp:20
void check_positive(const char *function, const char *name, const T_y &y)
Check if y is positive.
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.