Stan Math Library  2.20.0
reverse mode automatic differentiation
read_corr_L.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_MAT_FUN_READ_CORR_L_HPP
2 #define STAN_MATH_PRIM_MAT_FUN_READ_CORR_L_HPP
3 
8 #include <cstddef>
9 
10 namespace stan {
11 namespace math {
12 
34 template <typename T>
35 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> read_corr_L(
36  const Eigen::Array<T, Eigen::Dynamic, 1>& CPCs, // on (-1, 1)
37  size_t K) {
38  Eigen::Array<T, Eigen::Dynamic, 1> temp;
39  Eigen::Array<T, Eigen::Dynamic, 1> acc(K - 1);
40  acc.setOnes();
41  // Cholesky factor of correlation matrix
42  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic> L(K, K);
43  L.setZero();
44 
45  size_t position = 0;
46  size_t pull = K - 1;
47 
48  L(0, 0) = 1.0;
49  L.col(0).tail(pull) = temp = CPCs.head(pull);
50  acc.tail(pull) = T(1.0) - temp.square();
51  for (size_t i = 1; i < (K - 1); i++) {
52  position += pull;
53  pull--;
54  temp = CPCs.segment(position, pull);
55  L(i, i) = sqrt(acc(i - 1));
56  L.col(i).tail(pull) = temp * acc.tail(pull).sqrt();
57  acc.tail(pull) *= T(1.0) - temp.square();
58  }
59  L(K - 1, K - 1) = sqrt(acc(K - 2));
60  return L.matrix();
61 }
62 
87 template <typename T>
88 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> read_corr_L(
89  const Eigen::Array<T, Eigen::Dynamic, 1>& CPCs, size_t K, T& log_prob) {
90  Eigen::Matrix<T, Eigen::Dynamic, 1> values(CPCs.rows() - 1);
91  size_t pos = 0;
92  // no need to abs() because this Jacobian determinant
93  // is strictly positive (and triangular)
94  // see inverse of Jacobian in equation 11 of LKJ paper
95  for (size_t k = 1; k <= (K - 2); k++)
96  for (size_t i = k + 1; i <= K; i++) {
97  values(pos) = (K - k - 1) * log1m(square(CPCs(pos)));
98  pos++;
99  }
100 
101  log_prob += 0.5 * sum(values);
102  return read_corr_L(CPCs, K);
103 }
104 
105 } // namespace math
106 } // namespace stan
107 #endif
fvar< T > sum(const std::vector< fvar< T > > &m)
Return the sum of the entries of the specified standard vector.
Definition: sum.hpp:20
fvar< T > sqrt(const fvar< T > &x)
Definition: sqrt.hpp:13
fvar< T > square(const fvar< T > &x)
Definition: square.hpp:12
fvar< T > log1m(const fvar< T > &x)
Definition: log1m.hpp:12
Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > read_corr_L(const Eigen::Array< T, Eigen::Dynamic, 1 > &CPCs, size_t K)
Return the Cholesky factor of the correlation matrix of the specified dimensionality corresponding to...
Definition: read_corr_L.hpp:35

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