Stan Math Library  2.20.0
reverse mode automatic differentiation
cholesky_corr_constrain.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_MAT_FUN_CHOLESKY_CORR_CONSTRAIN_HPP
2 #define STAN_MATH_PRIM_MAT_FUN_CHOLESKY_CORR_CONSTRAIN_HPP
3 
9 #include <cmath>
10 
11 namespace stan {
12 namespace math {
13 
14 template <typename T>
15 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> cholesky_corr_constrain(
16  const Eigen::Matrix<T, Eigen::Dynamic, 1>& y, int K) {
17  using Eigen::Dynamic;
18  using Eigen::Matrix;
19  using std::sqrt;
20  int k_choose_2 = (K * (K - 1)) / 2;
21  check_size_match("cholesky_corr_constrain", "y.size()", y.size(),
22  "k_choose_2", k_choose_2);
23  Matrix<T, Dynamic, 1> z(k_choose_2);
24  for (int i = 0; i < k_choose_2; ++i)
25  z(i) = corr_constrain(y(i));
26  Matrix<T, Dynamic, Dynamic> x(K, K);
27  if (K == 0)
28  return x;
29  x.setZero();
30  x(0, 0) = 1;
31  int k = 0;
32  for (int i = 1; i < K; ++i) {
33  x(i, 0) = z(k++);
34  T sum_sqs(square(x(i, 0)));
35  for (int j = 1; j < i; ++j) {
36  x(i, j) = z(k++) * sqrt(1.0 - sum_sqs);
37  sum_sqs += square(x(i, j));
38  }
39  x(i, i) = sqrt(1.0 - sum_sqs);
40  }
41  return x;
42 }
43 
44 // FIXME to match above after debugged
45 template <typename T>
46 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> cholesky_corr_constrain(
47  const Eigen::Matrix<T, Eigen::Dynamic, 1>& y, int K, T& lp) {
48  using Eigen::Dynamic;
49  using Eigen::Matrix;
50  using std::sqrt;
51  int k_choose_2 = (K * (K - 1)) / 2;
52  check_size_match("cholesky_corr_constrain", "y.size()", y.size(),
53  "k_choose_2", k_choose_2);
54  Matrix<T, Dynamic, 1> z(k_choose_2);
55  for (int i = 0; i < k_choose_2; ++i)
56  z(i) = corr_constrain(y(i), lp);
57  Matrix<T, Dynamic, Dynamic> x(K, K);
58  if (K == 0)
59  return x;
60  x.setZero();
61  x(0, 0) = 1;
62  int k = 0;
63  for (int i = 1; i < K; ++i) {
64  x(i, 0) = z(k++);
65  T sum_sqs = square(x(i, 0));
66  for (int j = 1; j < i; ++j) {
67  lp += 0.5 * log1m(sum_sqs);
68  x(i, j) = z(k++) * sqrt(1.0 - sum_sqs);
69  sum_sqs += square(x(i, j));
70  }
71  x(i, i) = sqrt(1.0 - sum_sqs);
72  }
73  return x;
74 }
75 
76 } // namespace math
77 } // namespace stan
78 #endif
fvar< T > sqrt(const fvar< T > &x)
Definition: sqrt.hpp:13
T corr_constrain(const T &x)
Return the result of transforming the specified scalar to have a valid correlation value between -1 a...
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.
fvar< T > square(const fvar< T > &x)
Definition: square.hpp:12
Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > cholesky_corr_constrain(const Eigen::Matrix< T, Eigen::Dynamic, 1 > &y, int K)
fvar< T > log1m(const fvar< T > &x)
Definition: log1m.hpp:12

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