Automatic Differentiation
 
Loading...
Searching...
No Matches
multi_normal_prec_rng.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_PROB_MULTI_NORMAL_PREC_RNG_HPP
2#define STAN_MATH_PRIM_PROB_MULTI_NORMAL_PREC_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_prec_rng(const T_loc &mu, const Eigen::MatrixXd &S, RNG &rng) {
35 using boost::normal_distribution;
36 using boost::variate_generator;
37 static constexpr const char *function = "multi_normal_prec_rng";
38 check_positive(function, "Precision matrix rows", S.rows());
39
40 vector_seq_view<T_loc> mu_vec(mu);
41 check_positive(function, "number of location parameter vectors",
42 mu_vec.size());
43 size_t size_mu = mu_vec[0].size();
44
45 size_t N = mu_vec.size();
46
47 for (size_t i = 1; i < N; i++) {
48 check_size_match(function,
49 "Size of one of the vectors of "
50 "the location variable",
51 mu_vec[i].size(),
52 "Size of the first vector of the "
53 "location variable",
54 size_mu);
55 }
56 check_size_match(function, "Rows of location parameter", size_mu, "Rows of S",
57 S.rows());
58
59 for (size_t i = 0; i < N; i++) {
60 check_finite(function, "Location parameter", mu_vec[i]);
61 }
62 const auto &S_ref = to_ref(S);
63 check_finite(function, "Precision matrix", S_ref);
64 check_symmetric(function, "Precision matrix", S_ref);
65 Eigen::LLT<Eigen::MatrixXd> llt_of_S = S_ref.llt();
66 check_pos_definite(function, "precision matrix argument", llt_of_S);
67
69
70 variate_generator<RNG &, normal_distribution<>> std_normal_rng(
71 rng, normal_distribution<>(0, 1));
72
73 for (size_t n = 0; n < N; ++n) {
74 Eigen::VectorXd z(S.cols());
75 for (int i = 0; i < S.cols(); i++) {
76 z(i) = std_normal_rng();
77 }
78
79 output[n]
80 = as_column_vector_or_scalar(mu_vec[n]) + llt_of_S.matrixU().solve(z);
81 }
82
83 return output.data();
84}
85
86} // namespace math
87} // namespace stan
88#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_prec_rng(const T_loc &mu, const Eigen::MatrixXd &S, RNG &rng)
Return a multivariate normal random variate with the given location and precision using the specified...
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(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_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 ...