Stan  1.3
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 <boost/random/normal_distribution.hpp>
5 #include <boost/random/variate_generator.hpp>
6 
7 #include <stan/math/matrix_error_handling.hpp>
8 #include <stan/math/error_handling.hpp>
9 #include <stan/prob/constants.hpp>
10 #include <stan/prob/traits.hpp>
11 #include <stan/agrad/agrad.hpp>
12 #include <stan/meta/traits.hpp>
13 #include <stan/agrad/matrix.hpp>
14 #include <stan/math/matrix/columns_dot_product.hpp>
15 #include <stan/math/matrix/columns_dot_self.hpp>
16 #include <stan/math/matrix/dot_product.hpp>
17 #include <stan/math/matrix/dot_self.hpp>
18 #include <stan/math/matrix/log.hpp>
19 #include <stan/math/matrix/log_determinant.hpp>
20 #include <stan/math/matrix/mdivide_left_tri_low.hpp>
21 #include <stan/math/matrix/mdivide_left_spd.hpp>
22 #include <stan/math/matrix/multiply.hpp>
23 #include <stan/math/matrix/subtract.hpp>
24 #include <stan/math/matrix/sum.hpp>
25 
26 namespace stan {
27  namespace prob {
45  template <bool propto,
46  typename T_y, typename T_loc, typename T_covar,
47  class Policy>
48  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
49  multi_normal_cholesky_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
50  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
51  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L,
52  const Policy&) {
53  static const char* function = "stan::prob::multi_normal_cholesky_log(%1%)";
54 
55  using stan::math::mdivide_left_tri_low;
56  using stan::math::dot_self;
57  using stan::math::multiply;
58  using stan::math::subtract;
59  using stan::math::sum;
60 
61  using stan::math::check_size_match;
62  using stan::math::check_finite;
63  using stan::math::check_not_nan;
64  using stan::math::check_cov_matrix;
65  using boost::math::tools::promote_args;
66 
67  typename promote_args<T_y,T_loc,T_covar>::type lp(0.0);
68 
69  if (!check_size_match(function,
70  y.size(), "Size of random variable",
71  mu.size(), "size of location parameter",
72  &lp, Policy()))
73  return lp;
74  if (!check_size_match(function,
75  y.size(), "Size of random variable",
76  L.rows(), "rows of covariance parameter",
77  &lp, Policy()))
78  return lp;
79  if (!check_size_match(function,
80  y.size(), "Size of random variable",
81  L.cols(), "columns of covariance parameter",
82  &lp, Policy()))
83  return lp;
84  if (!check_finite(function, mu, "Location parameter", &lp, Policy()))
85  return lp;
86  if (!check_not_nan(function, y, "Random variable", &lp, Policy()))
87  return lp;
88 
89  if (y.rows() == 0)
90  return lp;
91 
92  if (include_summand<propto>::value)
93  lp += NEG_LOG_SQRT_TWO_PI * y.rows();
94 
95  if (include_summand<propto,T_covar>::value) {
96  Eigen::Matrix<T_covar,Eigen::Dynamic,1> L_log_diag = L.diagonal().array().log().matrix();
97  lp -= sum(L_log_diag);
98  }
99 
100  if (include_summand<propto,T_y,T_loc,T_covar>::value) {
101  Eigen::Matrix<typename
102  boost::math::tools::promote_args<T_y,T_loc>::type,
103  Eigen::Dynamic, 1> y_minus_mu(y.size());
104  for (int i = 0; i < y.size(); i++)
105  y_minus_mu(i) = y(i)-mu(i);
106  Eigen::Matrix<typename
107  boost::math::tools::promote_args<T_covar,T_loc,T_y>::type,
108  Eigen::Dynamic, 1>
109  half(mdivide_left_tri_low(L,y_minus_mu));
110  // FIXME: this code does not compile. revert after fixing subtract()
111  // Eigen::Matrix<typename
112  // boost::math::tools::promote_args<T_covar,T_loc,T_y>::type,
113  // Eigen::Dynamic, 1>
114  // half(mdivide_left_tri_low(L,subtract(y,mu)));
115  lp -= 0.5 * dot_self(half);
116  }
117  return lp;
118  }
119 
120  template <bool propto,
121  typename T_y, typename T_loc, typename T_covar>
122  inline
123  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
124  multi_normal_cholesky_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
125  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
126  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L) {
127  return multi_normal_cholesky_log<propto>(y,mu,L,
128  stan::math::default_policy());
129  }
130 
131  template <typename T_y, typename T_loc, typename T_covar,
132  class Policy>
133  inline
134  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
135  multi_normal_cholesky_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
136  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
137  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L,
138  const Policy&) {
139  return multi_normal_cholesky_log<false>(y,mu,L,Policy());
140  }
141 
142  template <typename T_y, typename T_loc, typename T_covar>
143  inline
144  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
145  multi_normal_cholesky_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
146  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
147  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L) {
148  return multi_normal_cholesky_log<false>(y,mu,L,
149  stan::math::default_policy());
150  }
151 
154  template <bool propto,
155  typename T_y, typename T_loc, typename T_covar,
156  class Policy>
157  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
158  multi_normal_cholesky_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
159  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
160  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L,
161  const Policy&) {
162  static const char* function = "stan::prob::multi_normal_cholesky_log(%1%)";
163 
164  using stan::math::mdivide_left_tri_low;
165  using stan::math::columns_dot_self;
166  using stan::math::multiply;
167  using stan::math::subtract;
168  using stan::math::sum;
169  using stan::math::log;
170 
171  using stan::math::check_size_match;
172  using stan::math::check_finite;
173  using stan::math::check_not_nan;
174  using stan::math::check_cov_matrix;
175  using boost::math::tools::promote_args;
176 
177  typename promote_args<T_y,T_loc,T_covar>::type lp(0.0);
178 
179  if (!check_size_match(function,
180  y.cols(), "Columns of random variable",
181  mu.rows(), "rows of location parameter",
182  &lp, Policy()))
183  return lp;
184  if (!check_size_match(function,
185  y.cols(), "Columns of random variable",
186  L.rows(), "rows of covariance parameter",
187  &lp, Policy()))
188  return lp;
189  if (!check_size_match(function,
190  y.cols(), "Columns of random variable",
191  L.cols(), "columns of covariance parameter",
192  &lp, Policy()))
193  return lp;
194  if (!check_finite(function, mu, "Location parameter", &lp, Policy()))
195  return lp;
196  if (!check_not_nan(function, y, "Random variable", &lp, Policy()))
197  return lp;
198 
199  if (y.cols() == 0)
200  return lp;
201 
202  if (include_summand<propto>::value)
203  lp += NEG_LOG_SQRT_TWO_PI * y.cols() * y.rows();
204 
205  if (include_summand<propto,T_covar>::value) {
206  Eigen::Matrix<T_covar,Eigen::Dynamic,1> L_log_diag = L.diagonal().array().log().matrix();
207  lp -= sum(L_log_diag) * y.rows();
208  }
209 
210  if (include_summand<propto,T_y,T_loc,T_covar>::value) {
211  Eigen::Matrix<T_loc, Eigen::Dynamic, Eigen::Dynamic> MU(y.rows(),y.cols());
212  for(typename Eigen::Matrix<T_loc, Eigen::Dynamic, Eigen::Dynamic>::size_type i = 0; i < y.rows(); i++)
213  MU.row(i) = mu;
214 
215  Eigen::Matrix<typename
216  boost::math::tools::promote_args<T_loc,T_y>::type,
217  Eigen::Dynamic,Eigen::Dynamic>
218  y_minus_MU(y.rows(), y.cols());
219  for (int i = 0; i < y.size(); i++)
220  y_minus_MU(i) = y(i)-MU(i);
221 
222  Eigen::Matrix<typename
223  boost::math::tools::promote_args<T_loc,T_y>::type,
224  Eigen::Dynamic,Eigen::Dynamic>
225  z(y_minus_MU.transpose()); // was =
226 
227  // FIXME: revert this code when subtract() is fixed.
228  // Eigen::Matrix<typename
229  // boost::math::tools::promote_args<T_loc,T_y>::type,
230  // Eigen::Dynamic,Eigen::Dynamic>
231  // z(subtract(y,MU).transpose()); // was =
232 
233  Eigen::Matrix<typename
234  boost::math::tools::promote_args<T_covar,T_loc,T_y>::type,
235  Eigen::Dynamic,Eigen::Dynamic>
236  half(mdivide_left_tri_low(L,z));
237 
238  lp -= 0.5 * sum(columns_dot_self(half));
239  }
240  return lp;
241  }
242 
243  template <bool propto,
244  typename T_y, typename T_loc, typename T_covar>
245  inline
246  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
247  multi_normal_cholesky_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
248  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
249  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L) {
250  return multi_normal_cholesky_log<propto>(y,mu,L,
251  stan::math::default_policy());
252  }
253 
254  template <typename T_y, typename T_loc, typename T_covar,
255  class Policy>
256  inline
257  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
258  multi_normal_cholesky_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
259  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
260  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L,
261  const Policy&) {
262  return multi_normal_cholesky_log<false>(y,mu,L,Policy());
263  }
264 
265  template <typename T_y, typename T_loc, typename T_covar>
266  inline
267  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
268  multi_normal_cholesky_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
269  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
270  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L) {
271  return multi_normal_cholesky_log<false>(y,mu,L,
272  stan::math::default_policy());
273  }
274 
275  // MultiNormal(y|mu,Sigma) [y.rows() = mu.rows() = Sigma.rows();
276  // y.cols() = mu.cols() = 0;
277  // Sigma symmetric, non-negative, definite]
294  template <bool propto,
295  typename T_y, typename T_loc, typename T_covar,
296  class Policy>
297  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
298  multi_normal_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
299  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
300  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma,
301  const Policy&) {
302  static const char* function = "stan::prob::multi_normal_log(%1%)";
303  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type lp(0.0);
304 
305  using stan::math::check_not_nan;
306  using stan::math::check_size_match;
307  using stan::math::check_positive;
308  using stan::math::check_pos_definite;
309  using stan::math::check_finite;
310  using stan::math::check_symmetric;
311  using stan::math::dot_product;
312  using stan::math::mdivide_left_spd;
313  using stan::math::log_determinant;
314 
315  if (!check_size_match(function,
316  Sigma.rows(), "Rows of covariance parameter",
317  Sigma.cols(), "columns of covariance parameter",
318  &lp, Policy()))
319  return lp;
320  if (!check_positive(function, Sigma.rows(), "Covariance matrix rows", &lp, Policy()))
321  return lp;
322  if (!check_symmetric(function, Sigma, "Covariance matrix", &lp, Policy()))
323  return lp;
324  if (!check_pos_definite(function, Sigma, "Covariance matrix", &lp, Policy()))
325  return lp;
326  if (!check_size_match(function,
327  y.size(), "Size of random variable",
328  mu.size(), "size of location parameter",
329  &lp, Policy()))
330  return lp;
331  if (!check_size_match(function,
332  y.size(), "Size of random variable",
333  Sigma.rows(), "rows of covariance parameter",
334  &lp, Policy()))
335  return lp;
336  if (!check_size_match(function,
337  y.size(), "Size of random variable",
338  Sigma.cols(), "columns of covariance parameter",
339  &lp, Policy()))
340  return lp;
341  if (!check_finite(function, mu, "Location parameter", &lp, Policy()))
342  return lp;
343  if (!check_not_nan(function, y, "Random variable", &lp, Policy()))
344  return lp;
345 
346  if (y.rows() == 0)
347  return lp;
348 
349  if (include_summand<propto>::value)
350  lp += NEG_LOG_SQRT_TWO_PI * y.rows();
351 
352  if (include_summand<propto,T_covar>::value) {
353  lp -= 0.5 * log_determinant(Sigma);
354  }
355 
356  if (include_summand<propto,T_y,T_loc,T_covar>::value) {
357  Eigen::Matrix<typename
358  boost::math::tools::promote_args<T_y,T_loc>::type,
359  Eigen::Dynamic, 1> y_minus_mu(y.size());
360  for (int i = 0; i < y.size(); i++)
361  y_minus_mu(i) = y(i)-mu(i);
362  Eigen::Matrix<typename
363  boost::math::tools::promote_args<T_covar,T_loc,T_y>::type,
364  Eigen::Dynamic, 1> Sinv_y_minus_mu(mdivide_left_spd(Sigma,y_minus_mu));
365  lp -= 0.5 * dot_product(y_minus_mu,Sinv_y_minus_mu);
366  }
367  return lp;
368  }
369 
370  template <bool propto,
371  typename T_y, typename T_loc, typename T_covar>
372  inline
373  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
374  multi_normal_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
375  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
376  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma) {
377  return multi_normal_log<propto>(y,mu,Sigma,stan::math::default_policy());
378  }
379 
380 
381  template <typename T_y, typename T_loc, typename T_covar,
382  class Policy>
383  inline
384  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
385  multi_normal_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
386  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
387  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma,
388  const Policy&){
389  return multi_normal_log<false>(y,mu,Sigma,Policy());
390  }
391 
392 
393  template <typename T_y, typename T_loc, typename T_covar>
394  inline
395  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
396  multi_normal_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
397  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
398  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma) {
399  return multi_normal_log<false>(y,mu,Sigma,stan::math::default_policy());
400  }
401 
405  template <bool propto,
406  typename T_y, typename T_loc, typename T_covar,
407  class Policy>
408  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
409  multi_normal_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
410  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
411  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma,
412  const Policy&) {
413  static const char* function = "stan::prob::multi_normal_log(%1%)";
414  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type lp(0.0);
415 
416  using stan::math::check_size_match;
417  using stan::math::check_positive;
418  using stan::math::check_pos_definite;
419  using stan::math::check_finite;
420  using stan::math::check_symmetric;
421  using stan::math::check_not_nan;
422  using stan::math::sum;
423  using stan::math::mdivide_left_spd;
424  using stan::math::log_determinant;
425  using stan::math::columns_dot_product;
426 
427  if (!check_size_match(function,
428  Sigma.rows(), "Rows of covariance matrix",
429  Sigma.cols(), "columns of covariance matrix",
430  &lp, Policy()))
431  return lp;
432  if (!check_positive(function, Sigma.rows(), "Covariance matrix rows", &lp, Policy()))
433  return lp;
434  if (!check_symmetric(function, Sigma, "Covariance matrix", &lp, Policy()))
435  return lp;
436  if (!check_pos_definite(function, Sigma, "Covariance matrix", &lp, Policy()))
437  return lp;
438  if (!check_size_match(function,
439  y.cols(), "Columns of random variable",
440  mu.rows(), "rows of location parameter",
441  &lp, Policy()))
442  return lp;
443  if (!check_size_match(function,
444  y.cols(), "Columns of random variable",
445  Sigma.rows(), "rows of covariance parameter",
446  &lp, Policy()))
447  return lp;
448  if (!check_size_match(function,
449  y.cols(), "Columns of random variable",
450  Sigma.cols(), "columns of covariance parameter",
451  &lp, Policy()))
452  return lp;
453  if (!check_finite(function, mu, "Location parameter", &lp, Policy()))
454  return lp;
455  if (!check_not_nan(function, y, "Random variable", &lp, Policy()))
456  return lp;
457 
458  if (y.cols() == 0)
459  return lp;
460 
461  if (include_summand<propto>::value)
462  lp += NEG_LOG_SQRT_TWO_PI * y.cols() * y.rows();
463 
464  if (include_summand<propto,T_covar>::value) {
465  lp -= 0.5 * log_determinant(Sigma) * y.rows();
466  }
467 
468  if (include_summand<propto,T_y,T_loc,T_covar>::value) {
469  Eigen::Matrix<T_loc, Eigen::Dynamic, Eigen::Dynamic> MU(y.rows(),y.cols());
470  for(typename Eigen::Matrix<T_loc, Eigen::Dynamic, Eigen::Dynamic>::size_type i = 0; i < y.rows(); i++)
471  MU.row(i) = mu;
472 
473  Eigen::Matrix<typename
474  boost::math::tools::promote_args<T_loc,T_y>::type,
475  Eigen::Dynamic,Eigen::Dynamic> y_minus_MU(y.rows(), y.cols());
476 
477  for (int i = 0; i < y.size(); i++)
478  y_minus_MU(i) = y(i)-MU(i);
479 
480  Eigen::Matrix<typename
481  boost::math::tools::promote_args<T_loc,T_y>::type,
482  Eigen::Dynamic,Eigen::Dynamic> z(y_minus_MU.transpose()); // was =
483 
484  // FIXME: revert this code when subtract() is fixed.
485  // Eigen::Matrix<typename
486  // boost::math::tools::promote_args<T_loc,T_y>::type,
487  // Eigen::Dynamic,Eigen::Dynamic>
488  // z(subtract(y,MU).transpose()); // was =
489 
490  Eigen::Matrix<typename
491  boost::math::tools::promote_args<T_covar,T_loc,T_y>::type,
492  Eigen::Dynamic,Eigen::Dynamic> Sinv_z(mdivide_left_spd(Sigma,z));
493 
494  lp -= 0.5 * sum(columns_dot_product(z,Sinv_z));
495  }
496  return lp;
497  }
498 
499  template <bool propto,
500  typename T_y, typename T_loc, typename T_covar>
501  inline
502  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
503  multi_normal_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
504  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
505  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma) {
506  return multi_normal_log<propto>(y,mu,Sigma,stan::math::default_policy());
507  }
508 
509 
510  template <typename T_y, typename T_loc, typename T_covar,
511  class Policy>
512  inline
513  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
514  multi_normal_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
515  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
516  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma,
517  const Policy&){
518  return multi_normal_log<false>(y,mu,Sigma,Policy());
519  }
520 
521  template <typename T_y, typename T_loc, typename T_covar>
522  inline
523  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
524  multi_normal_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
525  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
526  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma) {
527  return multi_normal_log<false>(y,mu,Sigma,stan::math::default_policy());
528  }
529 
530  template <bool propto,
531  typename T_y, typename T_loc, typename T_covar,
532  class Policy>
533  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
534  multi_normal_prec_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
535  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
536  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma,
537  const Policy&) {
538  static const char* function = "stan::prob::multi_normal_prec_log(%1%)";
539  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type lp(0.0);
540 
541  using stan::math::check_not_nan;
542  using stan::math::check_symmetric;
543  using stan::math::check_size_match;
544  using stan::math::check_positive;
545  using stan::math::check_pos_definite;
546  using stan::math::check_finite;
547  using stan::math::log_determinant;
548  using stan::math::sum;
549  using stan::math::dot_product;
550  using stan::math::multiply;
551 
552  if (!check_size_match(function,
553  Sigma.rows(), "Rows of covariance parameter",
554  Sigma.cols(), "columns of covariance parameter",
555  &lp, Policy()))
556  return lp;
557  if (!check_positive(function, Sigma.rows(), "Covariance matrix rows", &lp, Policy()))
558  return lp;
559  if (!check_symmetric(function, Sigma, "Covariance matrix", &lp, Policy()))
560  return lp;
561  if (!check_pos_definite(function, Sigma, "Covariance matrix", &lp, Policy()))
562  return lp;
563  if (!check_size_match(function,
564  y.size(), "Size of random variable",
565  mu.size(), "size of location parameter",
566  &lp, Policy()))
567  return lp;
568  if (!check_size_match(function,
569  y.size(), "Size of random variable",
570  Sigma.rows(), "rows of covariance parameter",
571  &lp, Policy()))
572  return lp;
573  if (!check_size_match(function,
574  y.size(), "Size of random variable",
575  Sigma.cols(), "columns of covariance parameter",
576  &lp, Policy()))
577  return lp;
578  if (!check_finite(function, mu, "Location parameter", &lp, Policy()))
579  return lp;
580  if (!check_not_nan(function, y, "Random variable", &lp, Policy()))
581  return lp;
582 
583  if (y.rows() == 0)
584  return lp;
585 
586  if (include_summand<propto>::value)
587  lp += NEG_LOG_SQRT_TWO_PI * y.rows();
588 
589  if (include_summand<propto,T_covar>::value)
590  lp += 0.5*log_determinant(Sigma);
591 
592  if (include_summand<propto,T_y,T_loc,T_covar>::value) {
593  Eigen::Matrix<typename
594  boost::math::tools::promote_args<T_y,T_loc>::type,
595  Eigen::Dynamic, 1> y_minus_mu(y.size());
596  for (int i = 0; i < y.size(); i++)
597  y_minus_mu(i) = y(i)-mu(i);
598  Eigen::Matrix<typename
599  boost::math::tools::promote_args<T_covar,T_loc,T_y>::type,
600  Eigen::Dynamic, 1> Sinv_y_minus_mu(multiply(Sigma,y_minus_mu));
601  lp -= 0.5 * dot_product(y_minus_mu,Sinv_y_minus_mu);
602  }
603  return lp;
604  }
605 
606  template <bool propto,
607  typename T_y, typename T_loc, typename T_covar>
608  inline
609  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
610  multi_normal_prec_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
611  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
612  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma) {
613  return multi_normal_prec_log<propto>(y,mu,Sigma,stan::math::default_policy());
614  }
615 
616 
617  template <typename T_y, typename T_loc, typename T_covar,
618  class Policy>
619  inline
620  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
621  multi_normal_prec_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
622  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
623  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma,
624  const Policy&){
625  return multi_normal_prec_log<false>(y,mu,Sigma,Policy());
626  }
627 
628 
629  template <typename T_y, typename T_loc, typename T_covar>
630  inline
631  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
632  multi_normal_prec_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
633  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
634  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma) {
635  return multi_normal_prec_log<false>(y,mu,Sigma,stan::math::default_policy());
636  }
637 
641  template <bool propto,
642  typename T_y, typename T_loc, typename T_covar,
643  class Policy>
644  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
645  multi_normal_prec_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
646  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
647  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma,
648  const Policy&) {
649  static const char* function = "stan::prob::multi_normal_prec_log(%1%)";
650  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type lp(0.0);
651 
652  using stan::math::check_not_nan;
653  using stan::math::check_symmetric;
654  using stan::math::check_size_match;
655  using stan::math::check_positive;
656  using stan::math::check_pos_definite;
657  using stan::math::check_finite;
658  using stan::math::log_determinant;
659  using stan::math::sum;
660  using stan::math::columns_dot_product;
661  using stan::math::multiply;
662 
663  if (!check_size_match(function,
664  Sigma.rows(), "Rows of covariance matrix",
665  Sigma.cols(), "columns of covariance matrix",
666  &lp, Policy()))
667  return lp;
668  if (!check_positive(function, Sigma.rows(), "Covariance matrix rows", &lp, Policy()))
669  return lp;
670  if (!check_symmetric(function, Sigma, "Covariance matrix", &lp, Policy()))
671  return lp;
672  if (!check_pos_definite(function, Sigma, "Covariance matrix", &lp, Policy()))
673  return lp;
674  if (!check_size_match(function,
675  y.cols(), "Columns of random variable",
676  mu.rows(), "rows of location parameter",
677  &lp, Policy()))
678  return lp;
679  if (!check_size_match(function,
680  y.cols(), "Columns of random variable",
681  Sigma.rows(), "rows of covariance parameter",
682  &lp, Policy()))
683  return lp;
684  if (!check_size_match(function,
685  y.cols(), "Columns of random variable",
686  Sigma.cols(), "columns of covariance parameter",
687  &lp, Policy()))
688  return lp;
689  if (!check_finite(function, mu, "Location parameter", &lp, Policy()))
690  return lp;
691  if (!check_not_nan(function, y, "Random variable", &lp, Policy()))
692  return lp;
693 
694  if (y.cols() == 0)
695  return lp;
696 
697  if (include_summand<propto>::value)
698  lp += NEG_LOG_SQRT_TWO_PI * y.cols() * y.rows();
699 
700  if (include_summand<propto,T_covar>::value) {
701  lp += log_determinant(Sigma) * (0.5 * y.rows());
702  }
703 
704  if (include_summand<propto,T_y,T_loc,T_covar>::value) {
705  Eigen::Matrix<T_loc, Eigen::Dynamic, Eigen::Dynamic> MU(y.rows(),y.cols());
706  for(typename Eigen::Matrix<T_loc, Eigen::Dynamic, Eigen::Dynamic>::size_type i = 0; i < y.rows(); i++)
707  MU.row(i) = mu;
708 
709  Eigen::Matrix<typename
710  boost::math::tools::promote_args<T_loc,T_y>::type,
711  Eigen::Dynamic,Eigen::Dynamic> y_minus_MU(y.rows(), y.cols());
712 
713  for (int i = 0; i < y.size(); i++)
714  y_minus_MU(i) = y(i)-MU(i);
715 
716  Eigen::Matrix<typename
717  boost::math::tools::promote_args<T_loc,T_y>::type,
718  Eigen::Dynamic,Eigen::Dynamic> z(y_minus_MU.transpose()); // was =
719 
720  // FIXME: revert this code when subtract() is fixed.
721  // Eigen::Matrix<typename
722  // boost::math::tools::promote_args<T_loc,T_y>::type,
723  // Eigen::Dynamic,Eigen::Dynamic>
724  // z(subtract(y,MU).transpose()); // was =
725 
726  Eigen::Matrix<typename
727  boost::math::tools::promote_args<T_covar,T_loc,T_y>::type,
728  Eigen::Dynamic,Eigen::Dynamic> Sinv_y_minus_mu(multiply(Sigma,z));
729 
730  lp -= 0.5 * sum(columns_dot_product(z,Sinv_y_minus_mu));
731  }
732  return lp;
733  }
734 
735  template <bool propto,
736  typename T_y, typename T_loc, typename T_covar>
737  inline
738  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
739  multi_normal_prec_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
740  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
741  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma) {
742  return multi_normal_prec_log<propto>(y,mu,Sigma,stan::math::default_policy());
743  }
744 
745 
746  template <typename T_y, typename T_loc, typename T_covar,
747  class Policy>
748  inline
749  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
750  multi_normal_prec_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
751  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
752  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma,
753  const Policy&){
754  return multi_normal_prec_log<false>(y,mu,Sigma,Policy());
755  }
756 
757  template <typename T_y, typename T_loc, typename T_covar>
758  inline
759  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
760  multi_normal_prec_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
761  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
762  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma) {
763  return multi_normal_prec_log<false>(y,mu,Sigma,stan::math::default_policy());
764  }
765 
766  template <class RNG>
767  inline Eigen::VectorXd
768  multi_normal_rng(const Eigen::Matrix<double,Eigen::Dynamic,1>& mu,
769  const Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic>& S,
770  RNG& rng) {
771  using boost::variate_generator;
772  using boost::normal_distribution;
773  variate_generator<RNG&, normal_distribution<> >
774  std_normal_rng(rng, normal_distribution<>(0,1));
775 
776  Eigen::VectorXd z(S.cols());
777  for(int i = 0; i < S.cols(); i++)
778  z(i) = std_normal_rng();
779 
780  return mu + S.llt().matrixL() * z;
781  }
782  }
783 }
784 
785 #endif
     [ Stan Home Page ] © 2011–2013, Stan Development Team.