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