Stan  1.0
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
multi_normal.hpp
Go to the documentation of this file.
1 #ifndef __STAN__PROB__DISTRIBUTIONS__MULTIVARIATE__CONTINUOUS__MULTI_NORMAL_HPP__
2 #define __STAN__PROB__DISTRIBUTIONS__MULTIVARIATE__CONTINUOUS__MULTI_NORMAL_HPP__
3 
4 #include <stan/math/matrix_error_handling.hpp>
5 #include <stan/math/error_handling.hpp>
6 #include <stan/math/special_functions.hpp>
7 #include <stan/prob/constants.hpp>
8 #include <stan/prob/traits.hpp>
9 
10 namespace stan {
11  namespace prob {
29  template <bool propto,
30  typename T_y, typename T_loc, typename T_covar,
31  class Policy>
32  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
33  multi_normal_cholesky_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
34  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
35  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L,
36  const Policy&) {
37  static const char* function = "stan::prob::multi_normal_cholesky_log(%1%)";
38 
39  using stan::math::mdivide_left_tri_low;
40  using stan::math::dot_self;
41  using stan::math::multiply;
42  using stan::math::subtract;
43 
44  using stan::math::check_size_match;
45  using stan::math::check_finite;
46  using stan::math::check_not_nan;
47  using stan::math::check_cov_matrix;
48  using boost::math::tools::promote_args;
49 
50  typename promote_args<T_y,T_loc,T_covar>::type lp(0.0);
51 
52  if (!check_size_match(function, y.size(), mu.size(), &lp, Policy()))
53  return lp;
54  if (!check_size_match(function, y.size(), L.rows(), &lp, Policy()))
55  return lp;
56  if (!check_size_match(function, y.size(), L.cols(), &lp, Policy()))
57  return lp;
58  if (!check_finite(function, mu, "Location parameter", &lp, Policy()))
59  return lp;
60  if (!check_not_nan(function, y, "Random variable", &lp, Policy()))
61  return lp;
62 
63  if (y.rows() == 0)
64  return lp;
65 
66  if (include_summand<propto>::value)
67  lp += NEG_LOG_SQRT_TWO_PI * y.rows();
68 
69  if (include_summand<propto,T_covar>::value)
70  lp -= L.diagonal().array().log().sum();
71 
72  if (include_summand<propto,T_y,T_loc,T_covar>::value) {
73  Eigen::Matrix<typename
74  boost::math::tools::promote_args<T_covar,T_loc,T_y>::type,
75  Eigen::Dynamic, 1>
76  half(mdivide_left_tri_low(L,subtract(y,mu)));
77 
78  lp -= 0.5 * dot_self(half);
79  }
80  return lp;
81  }
82 
83  template <bool propto,
84  typename T_y, typename T_loc, typename T_covar>
85  inline
86  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
87  multi_normal_cholesky_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
88  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
89  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L) {
90  return multi_normal_cholesky_log<propto>(y,mu,L,
91  stan::math::default_policy());
92  }
93 
94  template <typename T_y, typename T_loc, typename T_covar,
95  class Policy>
96  inline
97  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
98  multi_normal_cholesky_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
99  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
100  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L,
101  const Policy&) {
102  return multi_normal_cholesky_log<false>(y,mu,L,Policy());
103  }
104 
105  template <typename T_y, typename T_loc, typename T_covar>
106  inline
107  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
108  multi_normal_cholesky_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
109  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
110  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L) {
111  return multi_normal_cholesky_log<false>(y,mu,L,
112  stan::math::default_policy());
113  }
114 
117  template <bool propto,
118  typename T_y, typename T_loc, typename T_covar,
119  class Policy>
120  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
121  multi_normal_cholesky_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
122  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
123  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L,
124  const Policy&) {
125  static const char* function = "stan::prob::multi_normal_cholesky_log(%1%)";
126 
127  using stan::math::mdivide_left_tri_low;
128  using stan::math::columns_dot_self;
129  using stan::math::multiply;
130  using stan::math::subtract;
131 
132  using stan::math::check_size_match;
133  using stan::math::check_finite;
134  using stan::math::check_not_nan;
135  using stan::math::check_cov_matrix;
136  using boost::math::tools::promote_args;
137 
138  typename promote_args<T_y,T_loc,T_covar>::type lp(0.0);
139 
140  if (!check_size_match(function, y.cols(), mu.rows(), &lp, Policy()))
141  return lp;
142  if (!check_size_match(function, y.cols(), L.rows(), &lp, Policy()))
143  return lp;
144  if (!check_size_match(function, y.cols(), L.cols(), &lp, Policy()))
145  return lp;
146  if (!check_finite(function, mu, "Location parameter", &lp, Policy()))
147  return lp;
148  if (!check_not_nan(function, y, "Random variable", &lp, Policy()))
149  return lp;
150 
151  if (y.cols() == 0)
152  return lp;
153 
154  if (include_summand<propto>::value)
155  lp += NEG_LOG_SQRT_TWO_PI * y.cols() * y.rows();
156 
157  if (include_summand<propto,T_covar>::value)
158  lp -= L.diagonal().array().log().sum() * y.rows();
159 
160  if (include_summand<propto,T_y,T_loc,T_covar>::value) {
161  Eigen::Matrix<T_loc, Eigen::Dynamic, Eigen::Dynamic> MU(y.rows(),y.cols());
162  for(typename Eigen::Matrix<T_loc, Eigen::Dynamic, Eigen::Dynamic>::size_type i = 0; i < y.rows(); i++)
163  MU.row(i) = mu;
164 
165  Eigen::Matrix<typename
166  boost::math::tools::promote_args<T_loc,T_y>::type,
167  Eigen::Dynamic,Eigen::Dynamic>
168  z(subtract(y,MU).transpose()); // was =
169 
170  Eigen::Matrix<typename
171  boost::math::tools::promote_args<T_covar,T_loc,T_y>::type,
172  Eigen::Dynamic,Eigen::Dynamic>
173  half(mdivide_left_tri_low(L,z));
174 
175  lp -= 0.5 * columns_dot_self(half).sum();
176  }
177  return lp;
178  }
179 
180  template <bool propto,
181  typename T_y, typename T_loc, typename T_covar>
182  inline
183  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
184  multi_normal_cholesky_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
185  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
186  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L) {
187  return multi_normal_cholesky_log<propto>(y,mu,L,
188  stan::math::default_policy());
189  }
190 
191  template <typename T_y, typename T_loc, typename T_covar,
192  class Policy>
193  inline
194  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
195  multi_normal_cholesky_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
196  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
197  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L,
198  const Policy&) {
199  return multi_normal_cholesky_log<false>(y,mu,L,Policy());
200  }
201 
202  template <typename T_y, typename T_loc, typename T_covar>
203  inline
204  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
205  multi_normal_cholesky_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
206  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
207  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L) {
208  return multi_normal_cholesky_log<false>(y,mu,L,
209  stan::math::default_policy());
210  }
211 
212  // MultiNormal(y|mu,Sigma) [y.rows() = mu.rows() = Sigma.rows();
213  // y.cols() = mu.cols() = 0;
214  // Sigma symmetric, non-negative, definite]
231  template <bool propto,
232  typename T_y, typename T_loc, typename T_covar,
233  class Policy>
234  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
235  multi_normal_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
236  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
237  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma,
238  const Policy&) {
239  static const char* function = "stan::prob::multi_normal_log(%1%)";
240  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type lp(0.0);
241 
242  using stan::math::check_size_match;
243  using stan::math::check_positive;
244  using stan::math::check_symmetric;
245 
246  if (!check_size_match(function, Sigma.rows(), Sigma.cols(), &lp, Policy()))
247  return lp;
248  if (!check_positive(function, Sigma.rows(), "Covariance matrix rows", &lp, Policy()))
249  return lp;
250  if (!check_symmetric(function, Sigma, "Covariance matrix", &lp, Policy()))
251  return lp;
252  Eigen::LLT< Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic> > LLT = Sigma.llt();
253  if (LLT.info() != Eigen::Success) {
254  lp = stan::math::policies::raise_domain_error<T_covar>(function,
255  "Covariance matrix is not positive definite (%1%)",
256  0,Policy());
257  return lp;
258  }
259  Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic> L(LLT.matrixL());
260  return multi_normal_cholesky_log<propto>(y,mu,L,Policy());
261  }
262 
263  template <bool propto,
264  typename T_y, typename T_loc, typename T_covar>
265  inline
266  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
267  multi_normal_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
268  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
269  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma) {
270  return multi_normal_log<propto>(y,mu,Sigma,stan::math::default_policy());
271  }
272 
273 
274  template <typename T_y, typename T_loc, typename T_covar,
275  class Policy>
276  inline
277  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
278  multi_normal_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
279  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
280  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma,
281  const Policy&){
282  return multi_normal_log<false>(y,mu,Sigma,Policy());
283  }
284 
285 
286  template <typename T_y, typename T_loc, typename T_covar>
287  inline
288  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
289  multi_normal_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
290  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
291  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma) {
292  return multi_normal_log<false>(y,mu,Sigma,stan::math::default_policy());
293  }
294 
298  template <bool propto,
299  typename T_y, typename T_loc, typename T_covar,
300  class Policy>
301  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
302  multi_normal_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
303  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
304  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma,
305  const Policy&) {
306  static const char* function = "stan::prob::multi_normal_log(%1%)";
307  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type lp(0.0);
308 
309  using stan::math::check_size_match;
310  using stan::math::check_positive;
311  using stan::math::check_symmetric;
312 
313  if (!check_size_match(function, Sigma.rows(), Sigma.cols(), &lp, Policy()))
314  return lp;
315  if (!check_positive(function, Sigma.rows(), "Covariance matrix rows", &lp, Policy()))
316  return lp;
317  if (!check_symmetric(function, Sigma, "Covariance matrix", &lp, Policy()))
318  return lp;
319  Eigen::LLT< Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic> > LLT = Sigma.llt();
320  if (LLT.info() != Eigen::Success) {
321  lp = stan::math::policies::raise_domain_error<T_covar>(function,
322  "Sigma is not positive definite (%1%)",
323  0,Policy());
324  return lp;
325  }
326  Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic> L(LLT.matrixL());
327  return multi_normal_cholesky_log<propto>(y,mu,L,Policy());
328  }
329 
330  template <bool propto,
331  typename T_y, typename T_loc, typename T_covar>
332  inline
333  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
334  multi_normal_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
335  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
336  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma) {
337  return multi_normal_log<propto>(y,mu,Sigma,stan::math::default_policy());
338  }
339 
340 
341  template <typename T_y, typename T_loc, typename T_covar,
342  class Policy>
343  inline
344  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
345  multi_normal_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
346  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
347  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma,
348  const Policy&){
349  return multi_normal_log<false>(y,mu,Sigma,Policy());
350  }
351 
352  template <typename T_y, typename T_loc, typename T_covar>
353  inline
354  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
355  multi_normal_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
356  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
357  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma) {
358  return multi_normal_log<false>(y,mu,Sigma,stan::math::default_policy());
359  }
360 
361  }
362 }
363 
364 #endif
     [ Stan Home Page ] © 2011–2012, Stan Development Team.