| Stan Math Library 2.11.0 reverse mode automatic differentiation

- Main Page
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

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
4 #include <stan/math/prim/mat/fun/Eigen.hpp>
5 #include <stan/math/prim/scal/fun/log1m.hpp>
6 #include <stan/math/prim/scal/fun/square.hpp>
7 #include <stan/math/prim/scal/fun/corr_constrain.hpp>
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
stan::math::sqrt
fvar< T > sqrt(const fvar< T > &x)
Definition: sqrt.hpp:15
stan
Definition: log_sum_exp.hpp:8
square.hpp
stan::math::square
fvar< T > square(const fvar< T > &x)
Definition: square.hpp:15
corr_constrain.hpp
log1m.hpp
stan::math::cholesky_corr_constrain
Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > cholesky_corr_constrain(const Eigen::Matrix< T, Eigen::Dynamic, 1 > &y, int K)
Definition: cholesky_corr_constrain.hpp:20
stan::math::domain_error
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.
Definition: domain_error.hpp:32
Eigen.hpp
stan::math::log1m
fvar< T > log1m(const fvar< T > &x)
Definition: log1m.hpp:16
stan::math::corr_constrain
T corr_constrain(const T x)
Return the result of transforming the specified scalar to have a valid correlation value between -1 a...
Definition: corr_constrain.hpp:25

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