Automatic Differentiation
 
Loading...
Searching...
No Matches
dirichlet_rng.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_PROB_DIRICHLET_RNG_HPP
2#define STAN_MATH_PRIM_PROB_DIRICHLET_RNG_HPP
3
8#include <boost/math/special_functions/gamma.hpp>
9#include <boost/random/gamma_distribution.hpp>
10#include <boost/random/uniform_real_distribution.hpp>
11#include <boost/random/variate_generator.hpp>
12#include <cmath>
13
14namespace stan {
15namespace math {
16
38template <class RNG>
39inline Eigen::VectorXd dirichlet_rng(
40 const Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha, RNG& rng) {
41 using boost::gamma_distribution;
42 using boost::variate_generator;
43 using boost::random::uniform_real_distribution;
44 using Eigen::VectorXd;
45 using std::exp;
46 using std::log;
47
48 // separate algorithm if any parameter is less than 1
49 if (alpha.minCoeff() < 1) {
50 variate_generator<RNG&, uniform_real_distribution<> > uniform_rng(
51 rng, uniform_real_distribution<>(0.0, 1.0));
52 VectorXd log_y(alpha.size());
53 for (int i = 0; i < alpha.size(); ++i) {
54 variate_generator<RNG&, gamma_distribution<> > gamma_rng(
55 rng, gamma_distribution<>(alpha(i) + 1, 1));
56 double log_u = log(uniform_rng());
57 log_y(i) = log(gamma_rng()) + log_u / alpha(i);
58 }
59 double log_sum_y = log_sum_exp(log_y);
60 VectorXd theta(alpha.size());
61 for (int i = 0; i < alpha.size(); ++i) {
62 theta(i) = exp(log_y(i) - log_sum_y);
63 }
64 return theta;
65 }
66
67 // standard normalized gamma algorithm
68 Eigen::VectorXd y(alpha.rows());
69 for (int i = 0; i < alpha.rows(); i++) {
70 variate_generator<RNG&, gamma_distribution<> > gamma_rng(
71 rng, gamma_distribution<>(alpha(i, 0), 1e-7));
72 y(i) = gamma_rng();
73 }
74 return y / y.sum();
75}
76
77} // namespace math
78} // namespace stan
79#endif
Eigen::VectorXd dirichlet_rng(const Eigen::Matrix< double, Eigen::Dynamic, 1 > &alpha, RNG &rng)
Return a draw from a Dirichlet distribution with specified parameters and pseudo-random number genera...
VectorBuilder< true, double, T_alpha, T_beta >::type uniform_rng(const T_alpha &alpha, const T_beta &beta, RNG &rng)
Return a uniform random variate for the given upper and lower bounds using the specified random numbe...
VectorBuilder< true, double, T_shape, T_inv >::type gamma_rng(const T_shape &alpha, const T_inv &beta, RNG &rng)
Return a gamma random variate for the given shape and inverse scale parameters using the specified ra...
Definition gamma_rng.hpp:33
static constexpr double e()
Return the base of the natural logarithm.
Definition constants.hpp:20
fvar< T > log(const fvar< T > &x)
Definition log.hpp:15
fvar< T > log_sum_exp(const fvar< T > &x1, const fvar< T > &x2)
fvar< T > exp(const fvar< T > &x)
Definition exp.hpp:13
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...