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

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