Stan  1.0
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/math/special_functions.hpp>
8 #include <stan/prob/traits.hpp>
9 
10 namespace stan {
11  namespace prob {
12 
13  using Eigen::Matrix;
14  using Eigen::Dynamic;
15  using Eigen::LLT;
16  using Eigen::NumericalIssue;
17 
18  template <typename T_shape>
19  T_shape do_lkj_constant(const T_shape& eta, const unsigned int& K) {
20  // Lewandowski, Kurowicka, and Joe (2009) equations 15 and 16
21 
22  if (stan::is_constant<typename stan::scalar_type<T_shape> >::value
23  && eta == 1.0) {
24  double sum = 0.0;
25  double constant = 0.0;
26  double beta_arg = 0.0;
27  for (unsigned int k = 1; k < K; k++) { // yes, go from 1 to K - 1
28  beta_arg = 0.5 * (k + 1.0);
29  constant += k * (2.0 * lgamma(beta_arg) - lgamma(2.0 * beta_arg));
30  sum += pow(static_cast<double>(k),2.0);
31  }
32  constant += sum * LOG_TWO;
33  return constant;
34  }
35  T_shape sum = 0.0;
36  T_shape constant = 0.0;
37  T_shape beta_arg;
38  for (unsigned int k = 1; k < K; k++) { // yes, go from 1 to K - 1
39  unsigned int diff = K - k;
40  beta_arg = eta + 0.5 * (diff - 1);
41  constant += diff * (2.0 * lgamma(beta_arg) - lgamma(2.0 * beta_arg));
42  sum += (2.0 * eta - 2.0 + diff) * diff;
43  }
44  constant += sum * LOG_TWO;
45  return constant;
46  }
47 
48  // LKJ_Corr(L|eta) [ L Cholesky factor of correlation matrix
49  // eta > 0; eta == 1 <-> uniform]
50  template <bool propto,
51  typename T_covar, typename T_shape,
52  class Policy>
53  typename boost::math::tools::promote_args<T_covar, T_shape>::type
54  lkj_corr_cholesky_log(
55  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L,
56  const T_shape& eta,
57  const Policy&) {
58  static const char* function
59  = "stan::prob::lkj_corr_cholesky_log(%1%)";
60 
61  using boost::math::tools::promote_args;
62  using stan::math::check_positive;
63 
64  typename promote_args<T_covar,T_shape>::type lp(0.0);
65  if (!check_positive(function, eta, "Shape parameter", &lp, Policy()))
66  return lp;
67 
68  const unsigned int K = L.rows();
69  if (K == 0)
70  return 0.0;
71 
72  if (include_summand<propto,T_shape>::value)
73  lp += do_lkj_constant(eta, K);
74  if (include_summand<propto,T_covar,T_shape>::value) {
75  if ( (eta == 1.0) &&
76  stan::is_constant<typename stan::scalar_type<T_shape> >::value)
77  return lp;
78  lp += (eta - 1.0) * 2.0 * L.diagonal().array().log().sum();
79  }
80 
81  return lp;
82  }
83 
84  template <bool propto,
85  typename T_covar, typename T_shape>
86  inline
87  typename boost::math::tools::promote_args<T_covar, T_shape>::type
88  lkj_corr_cholesky_log(
89  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L,
90  const T_shape& eta) {
91  return lkj_corr_cholesky_log<propto>(L,eta,stan::math::default_policy());
92  }
93 
94 
95  template <typename T_covar, typename T_shape,
96  class Policy>
97  inline
98  typename boost::math::tools::promote_args<T_covar, T_shape>::type
99  lkj_corr_cholesky_log(
100  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L,
101  const T_shape& eta,
102  const Policy&) {
103  return lkj_corr_cholesky_log<false>(L,eta,Policy());
104  }
105 
106  template <typename T_covar, typename T_shape>
107  inline
108  typename boost::math::tools::promote_args<T_covar, T_shape>::type
109  lkj_corr_cholesky_log(
110  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L,
111  const T_shape& eta) {
112  return lkj_corr_cholesky_log<false>(L,eta,stan::math::default_policy());
113  }
114 
115 
116 
117  // LKJ_Corr(y|eta) [ y correlation matrix (not covariance matrix)
118  // eta > 0; eta == 1 <-> uniform]
119  template <bool propto,
120  typename T_y, typename T_shape,
121  class Policy>
122  typename boost::math::tools::promote_args<T_y, T_shape>::type
123  lkj_corr_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
124  const T_shape& eta,
125  const Policy&) {
126  static const char* function = "stan::prob::lkj_corr_log(%1%)";
127 
128  using stan::math::check_size_match;
129  using stan::math::check_not_nan;
130  using stan::math::check_positive;
131  using stan::math::check_corr_matrix;
132  using boost::math::tools::promote_args;
133 
134  typename promote_args<T_y,T_shape>::type lp;
135  if (!check_positive(function, eta, "Shape parameter", &lp, Policy()))
136  return lp;
137  if (!check_size_match(function, y.rows(), y.cols(), &lp, Policy()))
138  return lp;
139  if (!check_not_nan(function, y, "Correlation matrix", &lp, Policy()))
140  return lp;
141  if (!check_corr_matrix(function, y, "Correlation matrix", &lp, Policy())) {
142  return lp;
143  }
144 
145  const unsigned int K = y.rows();
146  if (K == 0)
147  return 0.0;
148 
149  LLT< Matrix<T_y, Dynamic, Dynamic> > Cholesky = y.llt();
150  // FIXME: check_numerical_issue function?
151  if (Cholesky.info() == Eigen::NumericalIssue)
152  return lp;
153 
154  Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic> L = Cholesky.matrixL();
155  return lkj_corr_cholesky_log<propto>(L, eta, Policy());
156  }
157 
158 
159 
160 
161  template <bool propto,
162  typename T_y, typename T_shape>
163  inline
164  typename boost::math::tools::promote_args<T_y, T_shape>::type
165  lkj_corr_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
166  const T_shape& eta) {
167  return lkj_corr_log<propto>(y,eta,stan::math::default_policy());
168  }
169 
170 
171  template <typename T_y, typename T_shape,
172  class Policy>
173  inline
174  typename boost::math::tools::promote_args<T_y, T_shape>::type
175  lkj_corr_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
176  const T_shape& eta,
177  const Policy&) {
178  return lkj_corr_log<false>(y,eta,Policy());
179  }
180 
181 
182  template <typename T_y, typename T_shape>
183  inline
184  typename boost::math::tools::promote_args<T_y, T_shape>::type
185  lkj_corr_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
186  const T_shape& eta) {
187  return lkj_corr_log<false>(y,eta,stan::math::default_policy());
188  }
189 
190 
191  }
192 }
193 #endif
     [ Stan Home Page ] © 2011–2012, Stan Development Team.