Automatic Differentiation
 
Loading...
Searching...
No Matches
multi_normal_lpdf.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_PROB_MULTI_NORMAL_LPDF_HPP
2#define STAN_MATH_PRIM_PROB_MULTI_NORMAL_LPDF_HPP
3
20
21namespace stan {
22namespace math {
23
24template <bool propto, typename T_y, typename T_loc, typename T_covar,
25 require_any_not_vector_vt<is_stan_scalar, T_y, T_loc>* = nullptr,
27 T_y, T_loc, T_covar>* = nullptr>
29 const T_loc& mu,
30 const T_covar& Sigma) {
31 using T_covar_elem = typename scalar_type<T_covar>::type;
33 using T_partials_return = partials_return_t<T_y, T_loc, T_covar>;
34 using matrix_partials_t
35 = Eigen::Matrix<T_partials_return, Eigen::Dynamic, Eigen::Dynamic>;
36 using vector_partials_t = Eigen::Matrix<T_partials_return, Eigen::Dynamic, 1>;
37 using T_y_ref = ref_type_t<T_y>;
38 using T_mu_ref = ref_type_t<T_loc>;
39 using T_Sigma_ref = ref_type_t<T_covar>;
40
41 static const char* function = "multi_normal_lpdf";
42 check_positive(function, "Covariance matrix rows", Sigma.rows());
43
44 check_consistent_sizes_mvt(function, "y", y, "mu", mu);
45 size_t number_of_y = size_mvt(y);
46 size_t number_of_mu = size_mvt(mu);
47 if (number_of_y == 0 || number_of_mu == 0) {
48 return 0.0;
49 }
50
51 T_y_ref y_ref = y;
52 T_mu_ref mu_ref = mu;
53 T_Sigma_ref Sigma_ref = Sigma;
54 vector_seq_view<T_y_ref> y_vec(y_ref);
55 vector_seq_view<T_mu_ref> mu_vec(mu_ref);
56 const size_t size_vec = max_size_mvt(y, mu);
57 const int K = Sigma.rows();
58
59 int size_y = y_vec[0].size();
60 int size_mu = mu_vec[0].size();
61
62 // check size consistency of all random variables y
63 for (size_t i = 1, size_mvt_y = size_mvt(y); i < size_mvt_y; i++) {
64 check_size_match(function,
65 "Size of one of the vectors of "
66 "the random variable",
67 y_vec[i].size(),
68 "Size of the first vector of the "
69 "random variable",
70 size_y);
71 }
72 // check size consistency of all means mu
73 for (size_t i = 1, size_mvt_mu = size_mvt(mu); i < size_mvt_mu; i++) {
74 check_size_match(function,
75 "Size of one of the vectors of "
76 "the location variable",
77 mu_vec[i].size(),
78 "Size of the first vector of the "
79 "location variable",
80 size_mu);
81 }
82
83 check_size_match(function, "Size of random variable", size_y,
84 "size of location parameter", size_mu);
85 check_size_match(function, "Size of random variable", size_y,
86 "rows of covariance parameter", Sigma.rows());
87 check_size_match(function, "Size of random variable", size_y,
88 "columns of covariance parameter", Sigma.cols());
89
90 for (size_t i = 0; i < size_vec; i++) {
91 check_finite(function, "Location parameter", mu_vec[i]);
92 check_not_nan(function, "Random variable", y_vec[i]);
93 }
94 check_symmetric(function, "Covariance matrix", Sigma_ref);
95
96 auto ldlt_Sigma = make_ldlt_factor(value_of(Sigma_ref));
97
98 check_ldlt_factor(function, "LDLT_Factor of covariance parameter",
99 ldlt_Sigma);
100
101 if (unlikely(size_y == 0)) {
102 return T_return(0);
103 }
104
105 auto ops_partials = make_partials_propagator(y_ref, mu_ref, Sigma_ref);
106
107 T_partials_return logp(0);
108
110 logp += NEG_LOG_SQRT_TWO_PI * size_y * size_vec;
111 }
112
114 vector_partials_t half(size_vec);
115 vector_partials_t y_val_minus_mu_val(size_vec);
116
117 T_partials_return sum_lp_vec(0.0);
118 for (size_t i = 0; i < size_vec; i++) {
119 const auto& y_val = as_value_column_vector_or_scalar(y_vec[i]);
120 const auto& mu_val = as_value_column_vector_or_scalar(mu_vec[i]);
121 y_val_minus_mu_val = eval(y_val - mu_val);
122 half = mdivide_left_ldlt(ldlt_Sigma, y_val_minus_mu_val);
123
124 sum_lp_vec += dot_product(y_val_minus_mu_val, half);
125
127 partials_vec<0>(ops_partials)[i] += -half;
128 }
130 partials_vec<1>(ops_partials)[i] += half;
131 }
133 partials_vec<2>(ops_partials)[i] += 0.5 * half * half.transpose();
134 }
135 }
136
137 logp += -0.5 * sum_lp_vec;
138
139 // If the covariance is not autodiff, we can avoid computing a matrix
140 // inverse
143 logp += -0.5 * log_determinant_ldlt(ldlt_Sigma) * size_vec;
144 }
145 } else {
146 matrix_partials_t inv_Sigma
147 = mdivide_left_ldlt(ldlt_Sigma, Eigen::MatrixXd::Identity(K, K));
148
149 logp += -0.5 * log_determinant_ldlt(ldlt_Sigma) * size_vec;
150
151 partials<2>(ops_partials) += -0.5 * size_vec * inv_Sigma;
152 }
153 }
154
155 return ops_partials.build(logp);
156}
157
158template <bool propto, typename T_y, typename T_loc, typename T_covar,
161 T_y, T_loc, T_covar>* = nullptr>
163 const T_loc& mu,
164 const T_covar& Sigma) {
165 using T_covar_elem = typename scalar_type<T_covar>::type;
166 using T_return = return_type_t<T_y, T_loc, T_covar>;
167 using T_partials_return = partials_return_t<T_y, T_loc, T_covar>;
168 using matrix_partials_t
169 = Eigen::Matrix<T_partials_return, Eigen::Dynamic, Eigen::Dynamic>;
170 using vector_partials_t = Eigen::Matrix<T_partials_return, Eigen::Dynamic, 1>;
171 using T_y_ref = ref_type_t<T_y>;
172 using T_mu_ref = ref_type_t<T_loc>;
173 using T_Sigma_ref = ref_type_t<T_covar>;
174
175 static const char* function = "multi_normal_lpdf";
176 check_positive(function, "Covariance matrix rows", Sigma.rows());
177
178 T_y_ref y_ref = y;
179 T_mu_ref mu_ref = mu;
180 T_Sigma_ref Sigma_ref = Sigma;
181
182 decltype(auto) y_val = as_value_column_vector_or_scalar(y_ref);
183 decltype(auto) mu_val = as_value_column_vector_or_scalar(mu_ref);
184
185 const int size_y = y_ref.size();
186 const int size_mu = mu_ref.size();
187 const unsigned int K = Sigma.rows();
188
189 check_finite(function, "Location parameter", mu_val);
190 check_not_nan(function, "Random variable", y_val);
191
192 check_size_match(function, "Size of random variable", size_y,
193 "size of location parameter", size_mu);
194 check_size_match(function, "Size of random variable", size_y,
195 "rows of covariance parameter", Sigma.rows());
196 check_size_match(function, "Size of random variable", size_y,
197 "columns of covariance parameter", Sigma.cols());
198
199 check_symmetric(function, "Covariance matrix", Sigma_ref);
200
201 auto ldlt_Sigma = make_ldlt_factor(value_of(Sigma_ref));
202 check_ldlt_factor(function, "LDLT_Factor of covariance parameter",
203 ldlt_Sigma);
204
205 if (unlikely(size_y == 0)) {
206 return T_return(0);
207 }
208
209 auto ops_partials = make_partials_propagator(y_ref, mu_ref, Sigma_ref);
210
211 T_partials_return logp(0);
212
213 if (include_summand<propto>::value) {
214 logp += NEG_LOG_SQRT_TWO_PI * size_y;
215 }
216
217 if (include_summand<propto, T_y, T_loc, T_covar_elem>::value) {
218 vector_partials_t half(size_y);
219 vector_partials_t y_val_minus_mu_val = eval(y_val - mu_val);
220
221 // If the covariance is not autodiff, we can avoid computing a matrix
222 // inverse
223 if (is_constant<T_covar_elem>::value) {
224 half = mdivide_left_ldlt(ldlt_Sigma, y_val_minus_mu_val);
225
226 if (include_summand<propto>::value) {
227 logp += -0.5 * log_determinant_ldlt(ldlt_Sigma);
228 }
229 } else {
230 matrix_partials_t inv_Sigma
231 = mdivide_left_ldlt(ldlt_Sigma, Eigen::MatrixXd::Identity(K, K));
232
233 half.noalias() = inv_Sigma * y_val_minus_mu_val;
234
235 logp += -0.5 * log_determinant_ldlt(ldlt_Sigma);
236
237 edge<2>(ops_partials).partials_
238 += 0.5 * (half * half.transpose() - inv_Sigma);
239 }
240
241 logp += -0.5 * dot_product(y_val_minus_mu_val, half);
242
243 if (!is_constant_all<T_y>::value) {
244 partials<0>(ops_partials) += -half;
245 }
246 if (!is_constant_all<T_loc>::value) {
247 partials<1>(ops_partials) += half;
248 }
249 }
250
251 return ops_partials.build(logp);
252}
253
254template <typename T_y, typename T_loc, typename T_covar>
256 const T_y& y, const T_loc& mu, const T_covar& Sigma) {
257 return multi_normal_lpdf<false>(y, mu, Sigma);
258}
259
260} // namespace math
261} // namespace stan
262#endif
This class provides a low-cost wrapper for situations where you either need an Eigen Vector or RowVec...
#define unlikely(x)
void check_symmetric(const char *function, const char *name, const matrix_cl< T > &y)
Check if the matrix_cl is symmetric.
require_all_not_t< is_nonscalar_prim_or_rev_kernel_expression< std::decay_t< Types > >... > require_all_not_nonscalar_prim_or_rev_kernel_expression_t
Require none of the types satisfy is_nonscalar_prim_or_rev_kernel_expression.
int64_t size_mvt(const ScalarT &)
Provides the size of a multivariate argument.
Definition size_mvt.hpp:25
typename return_type< Ts... >::type return_type_t
Convenience type for the return type of the specified template parameters.
int64_t size(const T &m)
Returns the size (number of the elements) of a matrix_cl or var_value<matrix_cl<T>>.
Definition size.hpp:19
require_all_t< container_type_check_base< is_vector, value_type_t, TypeCheck, Check >... > require_all_vector_vt
Require all of the types satisfy is_vector.
value_type_t< T > log_determinant_ldlt(LDLT_factor< T > &A)
Returns log(abs(det(A))) given a LDLT_factor of A.
auto make_ldlt_factor(const T &A)
Make an LDLT_factor with matrix type plain_type_t<T>
return_type_t< T_y, T_loc, T_covar > multi_normal_lpdf(const T_y &y, const T_loc &mu, const T_covar &Sigma)
void check_consistent_sizes_mvt(const char *)
Trivial no input case, this function is a no-op.
T eval(T &&arg)
Inputs which have a plain_type equal to the own time are forwarded unmodified (for Eigen expressions ...
Definition eval.hpp:20
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition value_of.hpp:18
const double NEG_LOG_SQRT_TWO_PI
The value of minus the natural logarithm of the square root of , .
int64_t max_size_mvt(const T1 &x1, const Ts &... xs)
Calculate the size of the largest multivariate input.
Eigen::Matrix< value_type_t< EigMat >, Eigen::Dynamic, EigMat::ColsAtCompileTime > mdivide_left_ldlt(LDLT_factor< T > &A, const EigMat &b)
Returns the solution of the system Ax=b given an LDLT_factor of A.
void check_finite(const char *function, const char *name, const T_y &y)
Return true if all values in y are finite.
void check_ldlt_factor(const char *function, const char *name, LDLT_factor< T > &A)
Raise domain error if the specified LDLT factor is invalid.
void check_not_nan(const char *function, const char *name, const T_y &y)
Check if y is not NaN.
void check_positive(const char *function, const char *name, const T_y &y)
Check if y is positive.
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.
auto dot_product(const T_a &a, const T_b &b)
Returns the dot product of the specified vectors.
auto as_value_column_vector_or_scalar(T &&a)
Extract values from input argument and transform to a column vector.
auto make_partials_propagator(Ops &&... ops)
Construct an partials_propagator.
typename ref_type_if< true, T >::type ref_type_t
Definition ref_type.hpp:55
typename partials_return_type< Args... >::type partials_return_t
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...
Metaprogramming struct to detect whether a given type is constant in the mathematical sense (not the ...
Extends std::true_type when instantiated with zero or more template parameters, all of which extend t...
Template metaprogram to calculate whether a summand needs to be included in a proportional (log) prob...
std::decay_t< T > type