Stan  1.0
probability, sampling & optimization
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 
47  using stan::math::check_size_match;
48  using stan::math::check_finite;
49  using stan::math::check_not_nan;
50  using stan::math::check_cov_matrix;
51  using boost::math::tools::promote_args;
52 
53  typename promote_args<T_y,T_loc,T_covar>::type lp(0.0);
54 
55  if (!check_size_match(function, y.size(), mu.size(), &lp, Policy()))
56  return lp;
57  if (!check_size_match(function, y.size(), L.rows(), &lp, Policy()))
58  return lp;
59  if (!check_size_match(function, y.size(), L.cols(), &lp, Policy()))
60  return lp;
61  if (!check_finite(function, mu, "Location parameter", &lp, Policy()))
62  return lp;
63  if (!check_not_nan(function, y, "Random variable", &lp, Policy()))
64  return lp;
65 
66  if (y.rows() == 0)
67  return lp;
68 
69  if (include_summand<propto>::value)
70  lp += NEG_LOG_SQRT_TWO_PI * y.rows();
71 
72  if (include_summand<propto,T_covar>::value)
73  lp -= L.diagonal().array().log().sum();
74 
75  if (include_summand<propto,T_y,T_loc,T_covar>::value) {
76  Eigen::Matrix<typename
77  boost::math::tools::promote_args<T_y,T_loc>::type,
78  Eigen::Dynamic, 1> y_minus_mu(y.size());
79  for (int i = 0; i < y.size(); i++)
80  y_minus_mu(i) = y(i)-mu(i);
81  Eigen::Matrix<typename
82  boost::math::tools::promote_args<T_covar,T_loc,T_y>::type,
83  Eigen::Dynamic, 1>
84  half(mdivide_left_tri_low(L,y_minus_mu));
85  // FIXME: this code does not compile. revert after fixing subtract()
86  //Eigen::Matrix<typename
87  // boost::math::tools::promote_args<T_covar,T_loc,T_y>::type,
88  // Eigen::Dynamic, 1>
89  // half(mdivide_left_tri_low(L,subtract(y,mu)));
90  lp -= 0.5 * dot_self(half);
91  }
92  return lp;
93  }
94 
95  template <bool propto,
96  typename T_y, typename T_loc, typename T_covar>
97  inline
98  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
99  multi_normal_cholesky_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
100  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
101  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L) {
102  return multi_normal_cholesky_log<propto>(y,mu,L,
103  stan::math::default_policy());
104  }
105 
106  template <typename T_y, typename T_loc, typename T_covar,
107  class Policy>
108  inline
109  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
110  multi_normal_cholesky_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
111  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
112  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L,
113  const Policy&) {
114  return multi_normal_cholesky_log<false>(y,mu,L,Policy());
115  }
116 
117  template <typename T_y, typename T_loc, typename T_covar>
118  inline
119  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
120  multi_normal_cholesky_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
121  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
122  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L) {
123  return multi_normal_cholesky_log<false>(y,mu,L,
124  stan::math::default_policy());
125  }
126 
129  template <bool propto,
130  typename T_y, typename T_loc, typename T_covar,
131  class Policy>
132  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
133  multi_normal_cholesky_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
134  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
135  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L,
136  const Policy&) {
137  static const char* function = "stan::prob::multi_normal_cholesky_log(%1%)";
138 
139  using stan::math::mdivide_left_tri_low;
140  using stan::math::columns_dot_self;
141  using stan::math::multiply;
142  using stan::math::subtract;
143 
144  using stan::math::check_size_match;
145  using stan::math::check_finite;
146  using stan::math::check_not_nan;
147  using stan::math::check_cov_matrix;
148  using boost::math::tools::promote_args;
149 
150  typename promote_args<T_y,T_loc,T_covar>::type lp(0.0);
151 
152  if (!check_size_match(function, y.cols(), mu.rows(), &lp, Policy()))
153  return lp;
154  if (!check_size_match(function, y.cols(), L.rows(), &lp, Policy()))
155  return lp;
156  if (!check_size_match(function, y.cols(), L.cols(), &lp, Policy()))
157  return lp;
158  if (!check_finite(function, mu, "Location parameter", &lp, Policy()))
159  return lp;
160  if (!check_not_nan(function, y, "Random variable", &lp, Policy()))
161  return lp;
162 
163  if (y.cols() == 0)
164  return lp;
165 
166  if (include_summand<propto>::value)
167  lp += NEG_LOG_SQRT_TWO_PI * y.cols() * y.rows();
168 
169  if (include_summand<propto,T_covar>::value)
170  lp -= L.diagonal().array().log().sum() * y.rows();
171 
172  if (include_summand<propto,T_y,T_loc,T_covar>::value) {
173  Eigen::Matrix<T_loc, Eigen::Dynamic, Eigen::Dynamic> MU(y.rows(),y.cols());
174  for(typename Eigen::Matrix<T_loc, Eigen::Dynamic, Eigen::Dynamic>::size_type i = 0; i < y.rows(); i++)
175  MU.row(i) = mu;
176 
177  Eigen::Matrix<typename boost::math::tools::promote_args<T_loc,T_y>::type,
178  Eigen::Dynamic,Eigen::Dynamic>
179  y_minus_MU(y.rows(), y.cols());
180  for (int i = 0; i < y.size(); i++)
181  y_minus_MU(i) = y(i)-MU(i);
182  Eigen::Matrix<typename
183  boost::math::tools::promote_args<T_loc,T_y>::type,
184  Eigen::Dynamic,Eigen::Dynamic>
185  z(y_minus_MU.transpose()); // was =
186 
187  // FIXME: revert this code when subtract() is fixed.
188  // Eigen::Matrix<typename
189  // boost::math::tools::promote_args<T_loc,T_y>::type,
190  // Eigen::Dynamic,Eigen::Dynamic>
191  // z(subtract(y,MU).transpose()); // was =
192 
193  Eigen::Matrix<typename
194  boost::math::tools::promote_args<T_covar,T_loc,T_y>::type,
195  Eigen::Dynamic,Eigen::Dynamic>
196  half(mdivide_left_tri_low(L,z));
197 
198  lp -= 0.5 * columns_dot_self(half).sum();
199  }
200  return lp;
201  }
202 
203  template <bool propto,
204  typename T_y, typename T_loc, typename T_covar>
205  inline
206  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
207  multi_normal_cholesky_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
208  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
209  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L) {
210  return multi_normal_cholesky_log<propto>(y,mu,L,
211  stan::math::default_policy());
212  }
213 
214  template <typename T_y, typename T_loc, typename T_covar,
215  class Policy>
216  inline
217  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
218  multi_normal_cholesky_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
219  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
220  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L,
221  const Policy&) {
222  return multi_normal_cholesky_log<false>(y,mu,L,Policy());
223  }
224 
225  template <typename T_y, typename T_loc, typename T_covar>
226  inline
227  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
228  multi_normal_cholesky_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
229  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
230  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& L) {
231  return multi_normal_cholesky_log<false>(y,mu,L,
232  stan::math::default_policy());
233  }
234 
235  // MultiNormal(y|mu,Sigma) [y.rows() = mu.rows() = Sigma.rows();
236  // y.cols() = mu.cols() = 0;
237  // Sigma symmetric, non-negative, definite]
254  template <bool propto,
255  typename T_y, typename T_loc, typename T_covar,
256  class Policy>
257  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
258  multi_normal_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
259  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
260  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma,
261  const Policy&) {
262  static const char* function = "stan::prob::multi_normal_log(%1%)";
263  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type lp(0.0);
264 
265  using stan::math::check_size_match;
266  using stan::math::check_positive;
267  using stan::math::check_symmetric;
268 
269  if (!check_size_match(function, Sigma.rows(), Sigma.cols(), &lp, Policy()))
270  return lp;
271  if (!check_positive(function, Sigma.rows(), "Covariance matrix rows", &lp, Policy()))
272  return lp;
273  if (!check_symmetric(function, Sigma, "Covariance matrix", &lp, Policy()))
274  return lp;
275 
276  Eigen::LLT< Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic> > LLT = Sigma.llt();
277  if (LLT.info() != Eigen::Success) {
278  lp = stan::math::policies::raise_domain_error<T_covar>(function,
279  "Covariance matrix is not positive definite (%1%)",
280  0,Policy());
281  return lp;
282  }
283  Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic> L(LLT.matrixL());
284  return multi_normal_cholesky_log<propto>(y,mu,L,Policy());
285  }
286 
287  template <bool propto,
288  typename T_y, typename T_loc, typename T_covar>
289  inline
290  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
291  multi_normal_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
292  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
293  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma) {
294  return multi_normal_log<propto>(y,mu,Sigma,stan::math::default_policy());
295  }
296 
297 
298  template <typename T_y, typename T_loc, typename T_covar,
299  class Policy>
300  inline
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,1>& 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  return multi_normal_log<false>(y,mu,Sigma,Policy());
307  }
308 
309 
310  template <typename T_y, typename T_loc, typename T_covar>
311  inline
312  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
313  multi_normal_log(const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
314  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
315  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma) {
316  return multi_normal_log<false>(y,mu,Sigma,stan::math::default_policy());
317  }
318 
322  template <bool propto,
323  typename T_y, typename T_loc, typename T_covar,
324  class Policy>
325  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
326  multi_normal_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
327  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
328  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma,
329  const Policy&) {
330  static const char* function = "stan::prob::multi_normal_log(%1%)";
331  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type lp(0.0);
332 
333  using stan::math::check_size_match;
334  using stan::math::check_positive;
335  using stan::math::check_symmetric;
336 
337  if (!check_size_match(function, Sigma.rows(), Sigma.cols(), &lp, Policy()))
338  return lp;
339  if (!check_positive(function, Sigma.rows(), "Covariance matrix rows", &lp, Policy()))
340  return lp;
341  if (!check_symmetric(function, Sigma, "Covariance matrix", &lp, Policy()))
342  return lp;
343  Eigen::LLT< Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic> > LLT = Sigma.llt();
344  if (LLT.info() != Eigen::Success) {
345  lp = stan::math::policies::raise_domain_error<T_covar>(function,
346  "Sigma is not positive definite (%1%)",
347  0,Policy());
348  return lp;
349  }
350  Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic> L(LLT.matrixL());
351  return multi_normal_cholesky_log<propto>(y,mu,L,Policy());
352  }
353 
354  template <bool propto,
355  typename T_y, typename T_loc, typename T_covar>
356  inline
357  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
358  multi_normal_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
359  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
360  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma) {
361  return multi_normal_log<propto>(y,mu,Sigma,stan::math::default_policy());
362  }
363 
364 
365  template <typename T_y, typename T_loc, typename T_covar,
366  class Policy>
367  inline
368  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
369  multi_normal_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
370  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
371  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma,
372  const Policy&){
373  return multi_normal_log<false>(y,mu,Sigma,Policy());
374  }
375 
376  template <typename T_y, typename T_loc, typename T_covar>
377  inline
378  typename boost::math::tools::promote_args<T_y,T_loc,T_covar>::type
379  multi_normal_log(const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
380  const Eigen::Matrix<T_loc,Eigen::Dynamic,1>& mu,
381  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma) {
382  return multi_normal_log<false>(y,mu,Sigma,stan::math::default_policy());
383  }
384 
385  }
386 }
387 
388 #endif
size_type
internal::traits< Derived >::Index size_type
Definition: EigenDenseBaseAddons.hpp:1
error_handling.hpp
matrix_error_handling.hpp
special_functions.hpp
stan::agrad::dot_self
var dot_self(const Eigen::Matrix< var, R, C > &v)
Returns the dot product of a vector with itself.
Definition: matrix.hpp:702
stan::math::check_symmetric
bool check_symmetric(const char *function, const Eigen::Matrix< T_y, Eigen::Dynamic, Eigen::Dynamic > &y, const char *name, T_result *result, const Policy &)
Return true if the specified matrix is symmetric.
Definition: matrix_error_handling.hpp:70
stan::math::mdivide_left_tri_low
Eigen::Matrix< typename boost::math::tools::promote_args< T1, T2 >::type, R1, C2 > mdivide_left_tri_low(const Eigen::Matrix< T1, R1, C1 > &A, const Eigen::Matrix< T2, R2, C2 > &b)
Definition: matrix.hpp:1708
stan::math::columns_dot_self
Eigen::Matrix< T, Eigen::Dynamic, 1 > columns_dot_self(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &x)
Returns the dot product of each column of a matrix with itself.
Definition: matrix.hpp:861
stan::math::subtract
Eigen::Matrix< typename boost::math::tools::promote_args< T1, T2 >::type, R, C > subtract(const Eigen::Matrix< T1, R, C > &m1, const Eigen::Matrix< T2, R, C > &m2)
Return the result of subtracting the second specified matrix from the first specified matrix.
Definition: matrix.hpp:1300
stan::math::check_not_nan
bool check_not_nan(const char *function, const T_y &y, const char *name, T_result *result, const Policy &)
Checks if the variable y is nan.
Definition: error_handling.hpp:144
stan::math::check_cov_matrix
bool check_cov_matrix(const char *function, const Eigen::Matrix< T_y, Eigen::Dynamic, Eigen::Dynamic > &y, const char *name, T_result *result, const Policy &)
Return true if the specified matrix is a valid covariance matrix.
Definition: matrix_error_handling.hpp:206
stan::math::dot_self
double dot_self(const Eigen::Matrix< double, R, C > &v)
Returns the dot product of the specified vector with itself.
Definition: matrix.hpp:846
stan::math::multiply
Eigen::Matrix< double, R, C > multiply(const Eigen::Matrix< double, R, C > &m, double c)
Return specified matrix multiplied by specified scalar.
Definition: matrix.hpp:1460
stan::math::check_positive
bool check_positive(const char *function, const T_y &y, const char *name, T_result *result, const Policy &)
Definition: error_handling.hpp:740
stan::math::check_finite
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.
Definition: error_handling.hpp:224
stan::math::default_policy
boost::math::policies::policy default_policy
Default error-handling policy from Boost.
Definition: error_handling.hpp:30
stan::math::check_size_match
bool check_size_match(const char *function, T_size1 i, T_size2 j, T_result *result, const Policy &)
Definition: matrix_error_handling.hpp:16
stan::prob::multi_normal_cholesky_log
boost::math::tools::promote_args< T_y, T_loc, T_covar >::type multi_normal_cholesky_log(const Eigen::Matrix< T_y, Eigen::Dynamic, 1 > &y, const Eigen::Matrix< T_loc, Eigen::Dynamic, 1 > &mu, const Eigen::Matrix< T_covar, Eigen::Dynamic, Eigen::Dynamic > &L, const Policy &)
The log of the multivariate normal density for the given y, mu, and a Cholesky factor L of the varian...
Definition: multi_normal.hpp:36
stan::prob::multi_normal_log
boost::math::tools::promote_args< T_y, T_loc, T_covar >::type multi_normal_log(const Eigen::Matrix< T_y, Eigen::Dynamic, 1 > &y, const Eigen::Matrix< T_loc, Eigen::Dynamic, 1 > &mu, const Eigen::Matrix< T_covar, Eigen::Dynamic, Eigen::Dynamic > &Sigma, const Policy &)
The log of the multivariate normal density for the given y, mu, and variance matrix.
Definition: multi_normal.hpp:258
stan
Probability, optimization and sampling library.
Definition: agrad.cpp:6
stan::prob::include_summand
Template metaprogram to calculate whether a summand needs to be included in a proportional (log) prob...
Definition: traits.hpp:33
     [ Stan Home Page ] © 2011–2012, Stan Development Team.