1 #ifndef STAN_MATH_PRIM_FUN_READ_CORR_L_HPP
2 #define STAN_MATH_PRIM_FUN_READ_CORR_L_HPP
3
4 #include <stan/math/prim/fun/Eigen.hpp>
5 #include <stan/math/prim/fun/log1m.hpp>
6 #include <stan/math/prim/fun/sqrt.hpp>
7 #include <stan/math/prim/fun/square.hpp>
8 #include <stan/math/prim/fun/sum.hpp>
9 #include <cstddef>
10 #include <cmath>
11
12 namespace stan {
13 namespace math {
14
15 /**
16 * Return the Cholesky factor of the correlation matrix of the
17 * specified dimensionality corresponding to the specified
18 * canonical partial correlations.
19 *
20 * It is generally better to work with the Cholesky factor rather
21 * than the correlation matrix itself when the determinant,
22 * inverse, etc. of the correlation matrix is needed for some
23 * statistical calculation.
24 *
25 * <p>See <code>read_corr_matrix(Array, size_t, T)</code>
26 * for more information.
27 *
28 * @tparam T type of the array (must be derived from \c Eigen::ArrayBase and
29 * have one compile-time dimension equal to 1)
30 * @param CPCs The (K choose 2) canonical partial correlations in
31 * (-1, 1).
32 * @param K Dimensionality of correlation matrix.
33 * @return Cholesky factor of correlation matrix for specified
34 * canonical partial correlations.
35 */
36 template <typename T, require_eigen_vector_t<T>* = nullptr>
read_corr_L(const T & CPCs,size_t K)37 Eigen::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
75 /**
76 * Return the Cholesky factor of the correlation matrix of the
77 * specified dimensionality corresponding to the specified
78 * canonical partial correlations, incrementing the specified
79 * scalar reference with the log absolute determinant of the
80 * Jacobian of the transformation.
81 *
82 * <p>The implementation is Ben Goodrich's Cholesky
83 * factor-based approach to the C-vine method of:
84 *
85 * <ul><li> Daniel Lewandowski, Dorota Kurowicka, and Harry Joe,
86 * Generating random correlation matrices based on vines and
87 * extended onion method Journal of Multivariate Analysis 100
88 * (2009) 1989–2001 </li></ul>
89 *
90 * @tparam T type of the array (must be derived from \c Eigen::ArrayBase and
91 * have one compile-time dimension equal to 1)
92 * @param CPCs The (K choose 2) canonical partial correlations in
93 * (-1, 1).
94 * @param K Dimensionality of correlation matrix.
95 * @param log_prob Reference to variable to increment with the log
96 * Jacobian determinant.
97 * @return Cholesky factor of correlation matrix for specified
98 * partial correlations.
99 */
100 template <typename T, require_eigen_vector_t<T>* = nullptr>
read_corr_L(const T & CPCs,size_t K,value_type_t<T> & log_prob)101 Eigen::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
133