Automatic Differentiation
 
Loading...
Searching...
No Matches
multi_normal_rng.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_PROB_MULTI_NORMAL_RNG_HPP
2#define STAN_MATH_PRIM_PROB_MULTI_NORMAL_RNG_HPP
3
10#include <boost/random/normal_distribution.hpp>
11#include <boost/random/variate_generator.hpp>
12
13namespace stan {
14namespace math {
15
32template <typename T_loc, class RNG>
34multi_normal_rng(const T_loc& mu,
35 const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& S,
36 RNG& rng) {
37 using boost::normal_distribution;
38 using boost::variate_generator;
39 static constexpr const char* function = "multi_normal_rng";
40 check_positive(function, "Covariance matrix rows", S.rows());
41
42 vector_seq_view<T_loc> mu_vec(mu);
43 size_t size_mu = mu_vec[0].size();
44
45 size_t N = size_mvt(mu);
46 for (size_t i = 1; i < N; i++) {
47 check_size_match(function,
48 "Size of one of the vectors of "
49 "the location variable",
50 mu_vec[i].size(),
51 "Size of the first vector of the "
52 "location variable",
53 size_mu);
54 }
55
56 for (size_t i = 0; i < N; i++) {
57 check_finite(function, "Location parameter", mu_vec[i]);
58 }
59 const auto& S_ref = to_ref(S);
60 check_not_nan(function, "Covariance matrix", S_ref);
61 check_symmetric(function, "Covariance matrix", S_ref);
62 Eigen::LLT<Eigen::MatrixXd> llt_of_S = S_ref.llt();
63 check_pos_definite("multi_normal_rng", "covariance matrix argument",
64 llt_of_S);
65
67
68 variate_generator<RNG&, normal_distribution<> > std_normal_rng(
69 rng, normal_distribution<>(0, 1));
70
71 for (size_t n = 0; n < N; ++n) {
72 Eigen::VectorXd z(S.cols());
73 for (int i = 0; i < S.cols(); i++) {
74 z(i) = std_normal_rng();
75 }
76
77 output[n] = as_column_vector_or_scalar(mu_vec[n]) + llt_of_S.matrixL() * z;
78 }
79
80 return output.data();
81}
82
83} // namespace math
84} // namespace stan
85#endif
typename helper::type type
StdVectorBuilder allocates type T1 values to be used as intermediate values.
This class provides a low-cost wrapper for situations where you either need an Eigen Vector or RowVec...
void check_symmetric(const char *function, const char *name, const matrix_cl< T > &y)
Check if the matrix_cl is symmetric.
StdVectorBuilder< true, Eigen::VectorXd, T_loc >::type multi_normal_rng(const T_loc &mu, const Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > &S, RNG &rng)
Return a multivariate normal random variate with the given location and covariance using the specifie...
auto as_column_vector_or_scalar(T &&a)
as_column_vector_or_scalar of a kernel generator expression.
double std_normal_rng(RNG &rng)
Return a standard Normal random variate using the specified random number generator.
int64_t size_mvt(const ScalarT &)
Provides the size of a multivariate argument.
Definition size_mvt.hpp:25
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
void check_pos_definite(const char *function, const char *name, const EigMat &y)
Check if the specified square, symmetric matrix is positive definite.
ref_type_t< T && > to_ref(T &&a)
This evaluates expensive Eigen expressions.
Definition to_ref.hpp:17
void check_finite(const char *function, const char *name, const T_y &y)
Return true if all values in y are finite.
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.
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...