Stan  1.0
probability, sampling & optimization
lkj_cov.hpp
Go to the documentation of this file.
1 #ifndef __STAN__PROB__DISTRIBUTIONS__MULTIVARIATE__CONTINUOUS__LKJ_COV_HPP__
2 #define __STAN__PROB__DISTRIBUTIONS__MULTIVARIATE__CONTINUOUS__LKJ_COV_HPP__
3 
5 #include <stan/math/matrix.hpp>
9 #include <stan/meta/traits.hpp>
10 #include <stan/prob/traits.hpp>
11 
14 
15 namespace stan {
16 
17  namespace prob {
18 
19  // LKJ_cov(y|mu,sigma,eta) [ y covariance matrix (not correlation matrix)
20  // mu vector, sigma > 0 vector, eta > 0 ]
21  template <bool propto,
22  typename T_y, typename T_loc, typename T_scale, typename T_shape,
23  class Policy>
24  typename boost::math::tools::promote_args<T_y,T_loc,T_scale,T_shape>::type
25  lkj_cov_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
26  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
27  const Eigen::Matrix<T_scale,Eigen::Dynamic,1>& sigma,
28  const T_shape& eta,
29  const Policy&) {
30  static const char* function = "stan::prob::lkj_cov_log(%1%)";
31 
35  using boost::math::tools::promote_args;
36 
37  typename promote_args<T_y,T_loc,T_scale,T_shape>::type lp(0.0);
38  if (!check_size_match(function, mu.rows(), sigma.rows(), &lp, Policy()))
39  return lp;
40  if (!check_size_match(function, y.cols(), y.rows(), &lp, Policy()))
41  return lp;
42  if (!check_size_match(function, mu.rows(), y.rows(), &lp, Policy()))
43  return lp;
44  if (!check_positive(function, eta, "Shape parameter", &lp, Policy()))
45  return lp;
46  if (!check_finite(function, mu, "Location parameter", &lp, Policy()))
47  return lp;
48  if (!check_finite(function, sigma, "Scale parameter", &lp, Policy()))
49  return lp;
50  // FIXME: build vectorized versions
51  for (int m = 0; m < y.rows(); ++m)
52  for (int n = 0; n < y.cols(); ++n)
53  if (!check_finite(function, y(m,n), "Covariance matrix", &lp, Policy()))
54  return lp;
55 
56  const unsigned int K = y.rows();
57  const Eigen::Array<T_y,Eigen::Dynamic,1> sds
58  = y.diagonal().array().sqrt();
59  for (unsigned int k = 0; k < K; k++) {
60  lp += lognormal_log<propto>(sds(k), mu(k), sigma(k), Policy());
61  }
62  if (stan::is_constant<typename stan::scalar_type<T_shape> >::value
63  && eta == 1.0) {
64  // no need to rescale y into a correlation matrix
65  lp += lkj_corr_log<propto,T_y,T_shape,Policy>(y, eta, Policy());
66  return lp;
67  }
68  Eigen::DiagonalMatrix<T_y,Eigen::Dynamic> D(K);
69  D.diagonal() = sds.inverse();
70  lp += lkj_corr_log<propto,T_y,T_shape,Policy>(D * y * D, eta, Policy());
71  return lp;
72  }
73 
74  template <bool propto,
75  typename T_y, typename T_loc, typename T_scale, typename T_shape>
76  inline
77  typename boost::math::tools::promote_args<T_y,T_loc,T_scale,T_shape>::type
78  lkj_cov_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
79  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
80  const Eigen::Matrix<T_scale,Eigen::Dynamic,1>& sigma,
81  const T_shape& eta) {
82  return lkj_cov_log<propto>(y,mu,sigma,eta,stan::math::default_policy());
83  }
84 
85 
86  template <typename T_y, typename T_loc, typename T_scale, typename T_shape,
87  class Policy>
88  inline
89  typename boost::math::tools::promote_args<T_y,T_loc,T_scale,T_shape>::type
90  lkj_cov_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
91  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
92  const Eigen::Matrix<T_scale,Eigen::Dynamic,1>& sigma,
93  const T_shape& eta,
94  const Policy&) {
95  return lkj_cov_log<false>(y,mu,sigma,eta,Policy());
96  }
97 
98 
99  template <typename T_y, typename T_loc, typename T_scale, typename T_shape>
100  inline
101  typename boost::math::tools::promote_args<T_y,T_loc,T_scale,T_shape>::type
102  lkj_cov_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
103  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
104  const Eigen::Matrix<T_scale,Eigen::Dynamic,1>& sigma,
105  const T_shape& eta) {
106  return lkj_cov_log<false>(y,mu,sigma,eta,stan::math::default_policy());
107  }
108 
109 
110 
111  // LKJ_Cov(y|mu,sigma,eta) [ y covariance matrix (not correlation matrix)
112  // mu scalar, sigma > 0 scalar, eta > 0 ]
113  template <bool propto,
114  typename T_y, typename T_loc, typename T_scale, typename T_shape,
115  class Policy>
116  typename boost::math::tools::promote_args<T_y,T_loc,T_scale,T_shape>::type
117  lkj_cov_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
118  const T_loc& mu,
119  const T_scale& sigma,
120  const T_shape& eta,
121  const Policy&) {
122  static const char* function = "stan::prob::lkj_cov_log(%1%)";
123 
126  using boost::math::tools::promote_args;
127 
128  typename promote_args<T_y,T_loc,T_scale,T_shape>::type lp(0.0);
129  if (!check_positive(function, eta, "Shape parameter", &lp, Policy()))
130  return lp;
131  if (!check_finite(function, mu, "Location parameter", &lp, Policy()))
132  return lp;
133  if (!check_finite(function, sigma, "Scale parameter",
134  &lp, Policy()))
135  return lp;
136 
137  const unsigned int K = y.rows();
138  const Eigen::Array<T_y,Eigen::Dynamic,1> sds
139  = y.diagonal().array().sqrt();
140  for (unsigned int k = 0; k < K; k++) {
141  lp += lognormal_log<propto>(sds(k), mu, sigma, Policy());
142  }
143  if (stan::is_constant<typename stan::scalar_type<T_shape> >::value
144  && eta == 1.0) {
145  // no need to rescale y into a correlation matrix
146  lp += lkj_corr_log<propto>(y,eta,Policy());
147  return lp;
148  }
149  Eigen::DiagonalMatrix<T_y,Eigen::Dynamic> D(K);
150  D.diagonal() = sds.inverse();
151  lp += lkj_corr_log<propto,T_y,T_shape,Policy>(D * y * D, eta, Policy());
152  return lp;
153  }
154 
155  template <bool propto,
156  typename T_y, typename T_loc, typename T_scale, typename T_shape>
157  inline
158  typename boost::math::tools::promote_args<T_y,T_loc,T_scale,T_shape>::type
159  lkj_cov_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
160  const T_loc& mu,
161  const T_scale& sigma,
162  const T_shape& eta) {
163  return lkj_cov_log<propto>(y,mu,sigma,eta,stan::math::default_policy());
164  }
165 
166  template <typename T_y, typename T_loc, typename T_scale, typename T_shape,
167  class Policy>
168  inline
169  typename boost::math::tools::promote_args<T_y,T_loc,T_scale,T_shape>::type
170  lkj_cov_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
171  const T_loc& mu,
172  const T_scale& sigma,
173  const T_shape& eta,
174  const Policy&) {
175  return lkj_cov_log<false>(y,mu,sigma,eta,Policy());
176  }
177 
178  template <typename T_y, typename T_loc, typename T_scale, typename T_shape>
179  inline
180  typename boost::math::tools::promote_args<T_y,T_loc,T_scale,T_shape>::type
181  lkj_cov_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
182  const T_loc& mu,
183  const T_scale& sigma,
184  const T_shape& eta) {
185  return lkj_cov_log<false>(y,mu,sigma,eta,stan::math::default_policy());
186  }
187 
188 
189  }
190 }
191 #endif
bool check_positive(const char *function, const T_y &y, const char *name, T_result *result, const Policy &)
bool check_finite(const char *function, const T_y &y, const char *name, T_result *result, const Policy &)
Checks if the variable y is finite.
boost::math::policies::policy default_policy
Default error-handling policy from Boost.
bool check_size_match(const char *function, T_size1 i, T_size2 j, T_result *result, const Policy &)
boost::math::tools::promote_args< T_y, T_loc, T_scale, T_shape >::type lkj_cov_log(const Eigen::Matrix< T_y, Eigen::Dynamic, Eigen::Dynamic > &y, const Eigen::Matrix< T_loc, Eigen::Dynamic, 1 > &mu, const Eigen::Matrix< T_scale, Eigen::Dynamic, 1 > &sigma, const T_shape &eta, const Policy &)
Definition: lkj_cov.hpp:25
Probability, optimization and sampling library.
Definition: agrad.cpp:6
Metaprogramming struct to detect whether a given type is constant in the mathematical sense (not the ...
Definition: traits.hpp:25
Metaprogram structure to determine the base scalar type of a template argument.
Definition: traits.hpp:104

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