Stan  1.1
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,
138  y.rows(), "Rows of correlation matrix",
139  y.cols(), "columns of correlation matrix",
140  &lp, Policy()))
141  return lp;
142  if (!check_not_nan(function, y, "Correlation matrix", &lp, Policy()))
143  return lp;
144  if (!check_corr_matrix(function, y, "Correlation matrix", &lp, Policy())) {
145  return lp;
146  }
147 
148  const unsigned int K = y.rows();
149  if (K == 0)
150  return 0.0;
151 
152  LLT< Matrix<T_y, Dynamic, Dynamic> > Cholesky = y.llt();
153  // FIXME: check_numerical_issue function?
154  if (Cholesky.info() == Eigen::NumericalIssue)
155  return lp;
156 
157  Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic> L = Cholesky.matrixL();
158  return lkj_corr_cholesky_log<propto>(L, eta, Policy());
159  }
160 
161 
162 
163 
164  template <bool propto,
165  typename T_y, typename T_shape>
166  inline
167  typename boost::math::tools::promote_args<T_y, T_shape>::type
168  lkj_corr_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
169  const T_shape& eta) {
170  return lkj_corr_log<propto>(y,eta,stan::math::default_policy());
171  }
172 
173 
174  template <typename T_y, typename T_shape,
175  class Policy>
176  inline
177  typename boost::math::tools::promote_args<T_y, T_shape>::type
178  lkj_corr_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
179  const T_shape& eta,
180  const Policy&) {
181  return lkj_corr_log<false>(y,eta,Policy());
182  }
183 
184 
185  template <typename T_y, typename T_shape>
186  inline
187  typename boost::math::tools::promote_args<T_y, T_shape>::type
188  lkj_corr_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
189  const T_shape& eta) {
190  return lkj_corr_log<false>(y,eta,stan::math::default_policy());
191  }
192 
193 
194  }
195 }
196 #endif
     [ Stan Home Page ] © 2011–2013, Stan Development Team.