Stan  1.3
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
lkj_corr.hpp
Go to the documentation of this file.
1 #ifndef __STAN__PROB__DISTRIBUTIONS__MULTIVARIATE__CONTINUOUS__LKJ_CORR_HPP__
2 #define __STAN__PROB__DISTRIBUTIONS__MULTIVARIATE__CONTINUOUS__LKJ_CORR_HPP__
3 
4 #include <stan/prob/constants.hpp>
5 #include <stan/math/matrix_error_handling.hpp>
6 #include <stan/math/error_handling.hpp>
7 #include <stan/prob/traits.hpp>
8 #include <stan/prob/distributions/univariate/continuous/beta.hpp>
9 #include <stan/prob/transform.hpp>
10 
11 namespace stan {
12  namespace prob {
13 
14  template <typename T_shape>
15  T_shape do_lkj_constant(const T_shape& eta, const unsigned int& K) {
16 
17  // Lewandowski, Kurowicka, and Joe (2009) equations 15 and 16
18 
19  if (stan::is_constant<typename stan::scalar_type<T_shape> >::value
20  && eta == 1.0) {
21  double sum = 0.0;
22  double constant = 0.0;
23  double beta_arg = 0.0;
24  for (unsigned int k = 1; k < K; k++) { // yes, go from 1 to K - 1
25  beta_arg = 0.5 * (k + 1.0);
26  constant += k * (2.0 * lgamma(beta_arg) - lgamma(2.0 * beta_arg));
27  sum += pow(static_cast<double>(k),2.0);
28  }
29  constant += sum * LOG_TWO;
30  return constant;
31  }
32  T_shape sum = 0.0;
33  T_shape constant = 0.0;
34  T_shape beta_arg;
35  for (unsigned int k = 1; k < K; k++) { // yes, go from 1 to K - 1
36  unsigned int diff = K - k;
37  beta_arg = eta + 0.5 * (diff - 1);
38  constant += diff * (2.0 * lgamma(beta_arg) - lgamma(2.0 * beta_arg));
39  sum += (2.0 * eta - 2.0 + diff) * diff;
40  }
41  constant += sum * LOG_TWO;
42  return constant;
43  }
44 
45  // LKJ_Corr(L|eta) [ L Cholesky factor of correlation matrix
46  // eta > 0; eta == 1 <-> uniform]
47  template <bool propto,
48  typename T_covar, typename T_shape,
49  class Policy>
50  typename boost::math::tools::promote_args<T_covar, T_shape>::type
51  lkj_corr_cholesky_log(
52  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L,
53  const T_shape& eta,
54  const Policy&) {
55 
56  static const char* function
57  = "stan::prob::lkj_corr_cholesky_log(%1%)";
58 
59  using boost::math::tools::promote_args;
60  using stan::math::check_positive;
61 
62  typename promote_args<T_covar,T_shape>::type lp(0.0);
63  if (!check_positive(function, eta, "Shape parameter", &lp, Policy()))
64  return lp;
65 
66  const unsigned int K = L.rows();
67  if (K == 0)
68  return 0.0;
69 
70  if (include_summand<propto,T_shape>::value)
71  lp += do_lkj_constant(eta, K);
72  if (include_summand<propto,T_covar,T_shape>::value) {
73  if ( (eta == 1.0) &&
74  stan::is_constant<typename stan::scalar_type<T_shape> >::value)
75  return lp;
76  lp += (eta - 1.0) * 2.0 * L.diagonal().array().log().sum();
77  }
78 
79  return lp;
80  }
81 
82  template <bool propto,
83  typename T_covar, typename T_shape>
84  inline
85  typename boost::math::tools::promote_args<T_covar, T_shape>::type
86  lkj_corr_cholesky_log(
87  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L,
88  const T_shape& eta) {
89  return lkj_corr_cholesky_log<propto>(L,eta,stan::math::default_policy());
90  }
91 
92 
93  template <typename T_covar, typename T_shape,
94  class Policy>
95  inline
96  typename boost::math::tools::promote_args<T_covar, T_shape>::type
97  lkj_corr_cholesky_log(
98  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L,
99  const T_shape& eta,
100  const Policy&) {
101  return lkj_corr_cholesky_log<false>(L,eta,Policy());
102  }
103 
104  template <typename T_covar, typename T_shape>
105  inline
106  typename boost::math::tools::promote_args<T_covar, T_shape>::type
107  lkj_corr_cholesky_log(
108  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L,
109  const T_shape& eta) {
110  return lkj_corr_cholesky_log<false>(L,eta,stan::math::default_policy());
111  }
112 
113 
114 
115  // LKJ_Corr(y|eta) [ y correlation matrix (not covariance matrix)
116  // eta > 0; eta == 1 <-> uniform]
117  template <bool propto,
118  typename T_y, typename T_shape,
119  class Policy>
120  typename boost::math::tools::promote_args<T_y, T_shape>::type
121  lkj_corr_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
122  const T_shape& eta,
123  const Policy&) {
124  static const char* function = "stan::prob::lkj_corr_log(%1%)";
125 
126  using stan::math::check_size_match;
127  using stan::math::check_not_nan;
128  using stan::math::check_positive;
129  using stan::math::check_corr_matrix;
130  using boost::math::tools::promote_args;
131 
132  typename promote_args<T_y,T_shape>::type lp;
133  if (!check_positive(function, eta, "Shape parameter", &lp, Policy()))
134  return lp;
135  if (!check_size_match(function,
136  y.rows(), "Rows of correlation matrix",
137  y.cols(), "columns of correlation matrix",
138  &lp, Policy()))
139  return lp;
140  if (!check_not_nan(function, y, "Correlation matrix", &lp, Policy()))
141  return lp;
142  if (!check_corr_matrix(function, y, "Correlation matrix", &lp, Policy())) {
143  return lp;
144  }
145 
146  const unsigned int K = y.rows();
147  if (K == 0)
148  return 0.0;
149 
150  Eigen::LLT< Eigen::Matrix<T_y, Eigen::Dynamic, Eigen::Dynamic> > Cholesky = y.llt();
151  // FIXME: check_numerical_issue function?
152  if (Cholesky.info() == Eigen::NumericalIssue)
153  return lp;
154 
155  Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic> L = Cholesky.matrixL();
156  return lkj_corr_cholesky_log<propto>(L, eta, Policy());
157  }
158 
159 
160 
161 
162  template <bool propto,
163  typename T_y, typename T_shape>
164  inline
165  typename boost::math::tools::promote_args<T_y, T_shape>::type
166  lkj_corr_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
167  const T_shape& eta) {
168  return lkj_corr_log<propto>(y,eta,stan::math::default_policy());
169  }
170 
171 
172  template <typename T_y, typename T_shape,
173  class Policy>
174  inline
175  typename boost::math::tools::promote_args<T_y, T_shape>::type
176  lkj_corr_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
177  const T_shape& eta,
178  const Policy&) {
179  return lkj_corr_log<false>(y,eta,Policy());
180  }
181 
182 
183  template <typename T_y, typename T_shape>
184  inline
185  typename boost::math::tools::promote_args<T_y, T_shape>::type
186  lkj_corr_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
187  const T_shape& eta) {
188  return lkj_corr_log<false>(y,eta,stan::math::default_policy());
189  }
190 
191  template <class RNG>
192  inline Eigen::MatrixXd
193  lkj_corr_cholesky_rng(const size_t K,
194  const double eta,
195  RNG& rng) {
196  // Need checks
197  Eigen::ArrayXd CPCs( (K * (K - 1)) / 2 );
198  double alpha = eta + 0.5 * (K - 1);
199  unsigned int count = 0;
200  for (size_t i = 0; i < (K - 1); i++) {
201  alpha -= 0.5;
202  for (size_t j = i + 1; j < K; j++) {
203  CPCs(count) = 2.0 * stan::prob::beta_rng(alpha,alpha,rng) - 1.0;
204  count++;
205  }
206  }
207  return stan::prob::read_corr_L(CPCs, K);
208  }
209 
210  template <class RNG>
211  inline Eigen::MatrixXd
212  lkj_corr_rng(const size_t K,
213  const double eta,
214  RNG& rng) {
215 
216  using stan::math::multiply_lower_tri_self_transpose;
217  return multiply_lower_tri_self_transpose(
218  lkj_corr_cholesky_rng(K, eta, rng) );
219  }
220 
221  }
222 }
223 #endif
     [ Stan Home Page ] © 2011–2013, Stan Development Team.