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