Stan Math Library  2.20.0
reverse mode automatic differentiation
multi_normal_cholesky_lpdf.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_MAT_PROB_MULTI_NORMAL_CHOLESKY_LPDF_HPP
2 #define STAN_MATH_PRIM_MAT_PROB_MULTI_NORMAL_CHOLESKY_LPDF_HPP
3 
14 
15 namespace stan {
16 namespace math {
40 template <bool propto, typename T_y, typename T_loc, typename T_covar>
42  const T_y& y, const T_loc& mu, const T_covar& L) {
43  static const char* function = "multi_normal_cholesky_lpdf";
44  typedef typename scalar_type<T_covar>::type T_covar_elem;
45  typedef typename return_type<T_y, T_loc, T_covar>::type T_return;
47  T_partials_return;
48  typedef Eigen::Matrix<T_partials_return, Eigen::Dynamic, Eigen::Dynamic>
49  matrix_partials_t;
50  typedef Eigen::Matrix<T_partials_return, Eigen::Dynamic, 1> vector_partials_t;
51  typedef Eigen::Matrix<T_partials_return, 1, Eigen::Dynamic>
52  row_vector_partials_t;
53 
54  check_consistent_sizes_mvt(function, "y", y, "mu", mu);
55  size_t number_of_y = length_mvt(y);
56  size_t number_of_mu = length_mvt(mu);
57  if (number_of_y == 0 || number_of_mu == 0)
58  return 0;
59  vector_seq_view<T_y> y_vec(y);
60  vector_seq_view<T_loc> mu_vec(mu);
61  const size_t size_vec = max_size_mvt(y, mu);
62 
63  const int size_y = y_vec[0].size();
64  const int size_mu = mu_vec[0].size();
65  if (likely(size_vec > 1)) {
66  // check size consistency of all random variables y
67  int size_y_old = size_y;
68  for (size_t i = 1, size_ = length_mvt(y); i < size_; i++) {
69  int size_y_new = y_vec[i].size();
70  check_size_match(function,
71  "Size of one of the vectors of "
72  "the random variable",
73  size_y_new,
74  "Size of another vector of the "
75  "random variable",
76  size_y_old);
77  size_y_old = size_y_new;
78  }
79  // check size consistency of all means mu
80  int size_mu_old = size_mu;
81  for (size_t i = 1, size_ = length_mvt(mu); i < size_; i++) {
82  int size_mu_new = mu_vec[i].size();
83  check_size_match(function,
84  "Size of one of the vectors of "
85  "the location variable",
86  size_mu_new,
87  "Size of another vector of the "
88  "location variable",
89  size_mu_old);
90  size_mu_old = size_mu_new;
91  }
92  }
93 
94  check_size_match(function, "Size of random variable", size_y,
95  "size of location parameter", size_mu);
96  check_size_match(function, "Size of random variable", size_y,
97  "rows of covariance parameter", L.rows());
98  check_size_match(function, "Size of random variable", size_y,
99  "columns of covariance parameter", L.cols());
100 
101  for (size_t i = 0; i < size_vec; i++) {
102  check_finite(function, "Location parameter", mu_vec[i]);
103  check_not_nan(function, "Random variable", y_vec[i]);
104  }
105 
106  if (unlikely(size_y == 0))
107  return T_return(0);
108 
109  T_partials_return logp(0);
110  operands_and_partials<T_y, T_loc, T_covar> ops_partials(y, mu, L);
111 
113  logp += NEG_LOG_SQRT_TWO_PI * size_y * size_vec;
114 
115  const matrix_partials_t inv_L_dbl
116  = mdivide_left_tri<Eigen::Lower>(value_of(L));
117 
119  for (size_t i = 0; i < size_vec; i++) {
120  vector_partials_t y_minus_mu_dbl(size_y);
121  for (int j = 0; j < size_y; j++)
122  y_minus_mu_dbl(j) = value_of(y_vec[i](j)) - value_of(mu_vec[i](j));
123 
124  const row_vector_partials_t half
125  = (inv_L_dbl.template triangularView<Eigen::Lower>() * y_minus_mu_dbl)
126  .transpose();
127  const vector_partials_t scaled_diff
128  = (half * inv_L_dbl.template triangularView<Eigen::Lower>())
129  .transpose();
130 
131  logp -= 0.5 * dot_self(half);
132 
134  for (int j = 0; j < size_y; j++)
135  ops_partials.edge1_.partials_vec_[i](j) -= scaled_diff(j);
136  }
138  for (int j = 0; j < size_y; j++)
139  ops_partials.edge2_.partials_vec_[i](j) += scaled_diff(j);
140  }
142  ops_partials.edge3_.partials_ += scaled_diff * half;
143  }
144  }
145  }
146 
148  logp += inv_L_dbl.diagonal().array().log().sum() * size_vec;
150  ops_partials.edge3_.partials_ -= size_vec * inv_L_dbl.transpose();
151  }
152  }
153 
154  return ops_partials.build(logp);
155 }
156 
157 template <typename T_y, typename T_loc, typename T_covar>
159 multi_normal_cholesky_lpdf(const T_y& y, const T_loc& mu, const T_covar& L) {
160  return multi_normal_cholesky_lpdf<false>(y, mu, L);
161 }
162 
163 } // namespace math
164 } // namespace stan
165 #endif
return_type< T_y, T_loc, T_covar >::type multi_normal_cholesky_lpdf(const T_y &y, const T_loc &mu, const T_covar &L)
The log of the multivariate normal density for the given y, mu, and a Cholesky factor L of the varian...
void check_finite(const char *function, const char *name, const T_y &y)
Check if y is finite.
boost::math::tools::promote_args< double, typename partials_type< typename scalar_type< T >::type >::type, typename partials_return_type< T_pack... >::type >::type type
size_t max_size_mvt(const T1 &x1, const T2 &x2)
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition: value_of.hpp:17
Extends std::true_type when instantiated with zero or more template parameters, all of which extend t...
Definition: conjunction.hpp:14
This template builds partial derivatives with respect to a set of operands.
void check_size_match(const char *function, const char *name_i, T_size1 i, const char *name_j, T_size2 j)
Check if the provided sizes match.
fvar< T > dot_self(const Eigen::Matrix< fvar< T >, R, C > &v)
Definition: dot_self.hpp:13
Template metaprogram to calculate whether a summand needs to be included in a proportional (log) prob...
#define unlikely(x)
Definition: likely.hpp:9
boost::math::tools::promote_args< double, typename scalar_type< T >::type, typename return_type< Types_pack... >::type >::type type
Definition: return_type.hpp:36
void check_not_nan(const char *function, const char *name, const T_y &y)
Check if y is not NaN.
T_return_type build(double value)
Build the node to be stored on the autodiff graph.
const double NEG_LOG_SQRT_TWO_PI
Definition: constants.hpp:156
internal::ops_partials_edge< double, Op2 > edge2_
This class provides a low-cost wrapper for situations where you either need an Eigen Vector or RowVec...
void check_consistent_sizes_mvt(const char *function, const char *name1, const T1 &x1, const char *name2, const T2 &x2)
Check if the dimension of x1 is consistent with x2.
#define likely(x)
Definition: likely.hpp:8
size_t length_mvt(const Eigen::Matrix< T, R, C > &)
Definition: length_mvt.hpp:12
internal::ops_partials_edge< double, Op3 > edge3_
internal::ops_partials_edge< double, Op1 > edge1_

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