Automatic Differentiation
 
Loading...
Searching...
No Matches
read_corr_L.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_FUN_READ_CORR_L_HPP
2#define STAN_MATH_PRIM_FUN_READ_CORR_L_HPP
3
9#include <cstddef>
10#include <cmath>
11
12namespace stan {
13namespace math {
14
36template <typename T, require_eigen_vector_t<T>* = nullptr>
37Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic> read_corr_L(
38 const T& CPCs, // on (-1, 1)
39 size_t K) {
40 using T_scalar = value_type_t<T>;
41 if (K == 0) {
42 return {};
43 }
44 if (K == 1) {
45 return Eigen::Matrix<T_scalar, Eigen::Dynamic, Eigen::Dynamic>::Identity(1,
46 1);
47 }
48
49 using std::sqrt;
50 Eigen::Array<T_scalar, Eigen::Dynamic, 1> temp;
51 Eigen::Array<T_scalar, Eigen::Dynamic, 1> acc(K - 1);
52 acc.setOnes();
53 // Cholesky factor of correlation matrix
54 Eigen::Matrix<T_scalar, Eigen::Dynamic, Eigen::Dynamic> L(K, K);
55 L.setZero();
56
57 size_t position = 0;
58 size_t pull = K - 1;
59
60 L(0, 0) = 1.0;
61 L.col(0).tail(pull) = temp = CPCs.head(pull);
62 acc.tail(pull) = T_scalar(1.0) - temp.square();
63 for (size_t i = 1; i < (K - 1); i++) {
64 position += pull;
65 pull--;
66 temp = CPCs.segment(position, pull);
67 L(i, i) = sqrt(acc(i - 1));
68 L.col(i).tail(pull) = temp * acc.tail(pull).sqrt();
69 acc.tail(pull) *= T_scalar(1.0) - temp.square();
70 }
71 L(K - 1, K - 1) = sqrt(acc(K - 2));
72 return L;
73}
74
100template <typename T, require_eigen_vector_t<T>* = nullptr>
101Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic> read_corr_L(
102 const T& CPCs, size_t K, value_type_t<T>& log_prob) {
103 using T_scalar = value_type_t<T>;
104 if (K == 0) {
105 return {};
106 }
107 if (K == 1) {
108 return Eigen::Matrix<T_scalar, Eigen::Dynamic, Eigen::Dynamic>::Identity(1,
109 1);
110 }
111
112 const Eigen::Ref<const plain_type_t<T>>& CPCs_ref = CPCs;
113 size_t pos = 0;
114 T_scalar acc = 0;
115 // no need to abs() because this Jacobian determinant
116 // is strictly positive (and triangular)
117 // see inverse of Jacobian in equation 11 of LKJ paper
118 for (size_t k = 1; k <= (K - 2); k++) {
119 for (size_t i = k + 1; i <= K; i++) {
120 acc += (K - k - 1) * log1m(square(CPCs_ref(pos)));
121 pos++;
122 }
123 }
124
125 log_prob += 0.5 * acc;
126 return read_corr_L(CPCs_ref, K);
127}
128
129} // namespace math
130} // namespace stan
131
132#endif
typename value_type< T >::type value_type_t
Helper function for accessing underlying type.
fvar< T > sqrt(const fvar< T > &x)
Definition sqrt.hpp:17
fvar< T > log1m(const fvar< T > &x)
Definition log1m.hpp:12
Eigen::Matrix< value_type_t< T >, Eigen::Dynamic, Eigen::Dynamic > read_corr_L(const T &CPCs, size_t K)
Return the Cholesky factor of the correlation matrix of the specified dimensionality corresponding to...
fvar< T > square(const fvar< T > &x)
Definition square.hpp:12
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...