Stan Math Library  2.6.3
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros
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 
8 #include <cmath>
9 #include <iostream>
10 #include <stdexcept>
11 
12 namespace stan {
13 
14  namespace math {
15 
16  // CHOLESKY CORRELATION MATRIX
17 
18  template <typename T>
19  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>
20  cholesky_corr_constrain(const Eigen::Matrix<T, Eigen::Dynamic, 1>& y,
21  int K) {
22  using std::sqrt;
23  using Eigen::Matrix;
24  using Eigen::Dynamic;
25  using stan::math::square;
26  int k_choose_2 = (K * (K - 1)) / 2;
27  if (k_choose_2 != y.size()) {
28  throw std::domain_error("y is not a valid unconstrained cholesky "
29  "correlation matrix."
30  "Require (K choose 2) elements in y.");
31  }
32  Matrix<T, Dynamic, 1> z(k_choose_2);
33  for (int i = 0; i < k_choose_2; ++i)
34  z(i) = corr_constrain(y(i));
35  Matrix<T, Dynamic, Dynamic> x(K, K);
36  if (K == 0) return x;
37  T zero(0);
38  for (int j = 1; j < K; ++j)
39  for (int i = 0; i < j; ++i)
40  x(i, j) = zero;
41  x(0, 0) = 1;
42  int k = 0;
43  for (int i = 1; i < K; ++i) {
44  x(i, 0) = z(k++);
45  T sum_sqs(square(x(i, 0)));
46  for (int j = 1; j < i; ++j) {
47  x(i, j) = z(k++) * sqrt(1.0 - sum_sqs);
48  sum_sqs += square(x(i, j));
49  }
50  x(i, i) = sqrt(1.0 - sum_sqs);
51  }
52  return x;
53  }
54 
55  // FIXME to match above after debugged
56  template <typename T>
57  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>
58  cholesky_corr_constrain(const Eigen::Matrix<T, Eigen::Dynamic, 1>& y,
59  int K,
60  T& lp) {
61  using std::sqrt;
62  using Eigen::Matrix;
63  using Eigen::Dynamic;
64  using stan::math::log1m;
65  using stan::math::square;
66  int k_choose_2 = (K * (K - 1)) / 2;
67  if (k_choose_2 != y.size()) {
68  throw std::domain_error("y is not a valid unconstrained cholesky "
69  "correlation matrix."
70  " Require (K choose 2) elements in y.");
71  }
72  Matrix<T, Dynamic, 1> z(k_choose_2);
73  for (int i = 0; i < k_choose_2; ++i)
74  z(i) = corr_constrain(y(i), lp);
75  Matrix<T, Dynamic, Dynamic> x(K, K);
76  if (K == 0) return x;
77  T zero(0);
78  for (int j = 1; j < K; ++j)
79  for (int i = 0; i < j; ++i)
80  x(i, j) = zero;
81  x(0, 0) = 1;
82  int k = 0;
83  for (int i = 1; i < K; ++i) {
84  x(i, 0) = z(k++);
85  T sum_sqs = square(x(i, 0));
86  for (int j = 1; j < i; ++j) {
87  lp += 0.5 * log1m(sum_sqs);
88  x(i, j) = z(k++) * sqrt(1.0 - sum_sqs);
89  sum_sqs += square(x(i, j));
90  }
91  x(i, i) = sqrt(1.0 - sum_sqs);
92  }
93  return x;
94  }
95 
96  }
97 
98 }
99 
100 #endif
fvar< T > sqrt(const fvar< T > &x)
Definition: sqrt.hpp:15
fvar< T > square(const fvar< T > &x)
Definition: square.hpp:15
Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > cholesky_corr_constrain(const Eigen::Matrix< T, Eigen::Dynamic, 1 > &y, int K)
void domain_error(const char *function, const char *name, const T &y, const char *msg1, const char *msg2)
Throw a domain error with a consistently formatted message.
fvar< T > log1m(const fvar< T > &x)
Definition: log1m.hpp:16
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–2015, Stan Development Team.