Automatic Differentiation
 
Loading...
Searching...
No Matches
wishart_rng.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_PROB_WISHART_RNG_HPP
2#define STAN_MATH_PRIM_PROB_WISHART_RNG_HPP
3
9#include <cmath>
10
11namespace stan {
12namespace math {
13
14template <class RNG>
15inline Eigen::MatrixXd wishart_rng(double nu, const Eigen::MatrixXd& S,
16 RNG& rng) {
17 using Eigen::MatrixXd;
18 static constexpr const char* function = "wishart_rng";
19 index_type_t<MatrixXd> k = S.rows();
20 check_square(function, "scale parameter", S);
21 check_symmetric(function, "scale parameter", S);
22 check_greater(function, "degrees of freedom > dims - 1", nu, k - 1);
23
24 Eigen::LLT<Eigen::MatrixXd> llt_of_S = S.llt();
25 check_pos_definite("wishart_rng", "scale parameter", llt_of_S);
26
27 MatrixXd B = MatrixXd::Zero(k, k);
28 for (int j = 0; j < k; ++j) {
29 for (int i = 0; i < j; ++i) {
30 B(i, j) = normal_rng(0, 1, rng);
31 }
32 B(j, j) = std::sqrt(chi_square_rng(nu - j, rng));
33 }
34 return crossprod(B * llt_of_S.matrixU());
35}
36
37} // namespace math
38} // namespace stan
39#endif
void check_symmetric(const char *function, const char *name, const matrix_cl< T > &y)
Check if the matrix_cl is symmetric.
VectorBuilder< true, double, T_loc, T_scale >::type normal_rng(const T_loc &mu, const T_scale &sigma, RNG &rng)
Return a Normal random variate for the given location and scale using the specified random number gen...
VectorBuilder< true, double, T_deg >::type chi_square_rng(const T_deg &nu, RNG &rng)
Return a chi squared random variate with nu degrees of freedom using the specified random number gene...
void check_square(const char *function, const char *name, const T_y &y)
Check if the specified matrix is square.
typename index_type< T >::type index_type_t
void check_pos_definite(const char *function, const char *name, const EigMat &y)
Check if the specified square, symmetric matrix is positive definite.
matrix_cl< typename std::decay_t< T_A >::Scalar > crossprod(T_A &&A)
Returns the result of pre-multiplying a matrix by its own transpose.
Definition crossprod.hpp:21
Eigen::MatrixXd wishart_rng(double nu, const Eigen::MatrixXd &S, RNG &rng)
void check_greater(const char *function, const char *name, const T_y &y, const T_low &low, Idxs... idxs)
Throw an exception if y is not strictly greater than low.
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...