Stan Math Library  2.20.0
reverse mode automatic differentiation
beta_rng.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_SCAL_PROB_BETA_RNG_HPP
2 #define STAN_MATH_PRIM_SCAL_PROB_BETA_RNG_HPP
3 
5 #include <boost/random/gamma_distribution.hpp>
6 #include <boost/random/uniform_real_distribution.hpp>
7 #include <boost/random/variate_generator.hpp>
11 
12 namespace stan {
13 namespace math {
14 
33 template <typename T_shape1, typename T_shape2, class RNG>
35  const T_shape1 &alpha, const T_shape2 &beta, RNG &rng) {
36  using boost::random::gamma_distribution;
37  using boost::random::uniform_real_distribution;
38  using boost::variate_generator;
39  static const char *function = "beta_rng";
40 
41  check_positive_finite(function, "First shape parameter", alpha);
42  check_positive_finite(function, "Second shape parameter", beta);
43  check_consistent_sizes(function, "First shape parameter", alpha,
44  "Second shape Parameter", beta);
45 
46  scalar_seq_view<T_shape1> alpha_vec(alpha);
47  scalar_seq_view<T_shape2> beta_vec(beta);
48  size_t N = max_size(alpha, beta);
50 
51  variate_generator<RNG &, uniform_real_distribution<>> uniform_rng(
52  rng, uniform_real_distribution<>(0.0, 1.0));
53  for (size_t n = 0; n < N; ++n) {
54  // If alpha and beta are large, trust the usual ratio of gammas
55  // method for generating beta random variables. If any parameter
56  // is small, work in log space and use Marsaglia and Tsang's trick
57  if (alpha_vec[n] > 1.0 && beta_vec[n] > 1.0) {
58  variate_generator<RNG &, gamma_distribution<>> rng_gamma_alpha(
59  rng, gamma_distribution<>(alpha_vec[n], 1.0));
60  variate_generator<RNG &, gamma_distribution<>> rng_gamma_beta(
61  rng, gamma_distribution<>(beta_vec[n], 1.0));
62  double a = rng_gamma_alpha();
63  double b = rng_gamma_beta();
64  output[n] = a / (a + b);
65  } else {
66  variate_generator<RNG &, gamma_distribution<>> rng_gamma_alpha(
67  rng, gamma_distribution<>(alpha_vec[n] + 1, 1.0));
68  variate_generator<RNG &, gamma_distribution<>> rng_gamma_beta(
69  rng, gamma_distribution<>(beta_vec[n] + 1, 1.0));
70  double log_a = std::log(uniform_rng()) / alpha_vec[n]
71  + std::log(rng_gamma_alpha());
72  double log_b
73  = std::log(uniform_rng()) / beta_vec[n] + std::log(rng_gamma_beta());
74  double log_sum = log_sum_exp(log_a, log_b);
75  output[n] = std::exp(log_a - log_sum);
76  }
77  }
78 
79  return output.data();
80 }
81 
82 } // namespace math
83 } // namespace stan
84 #endif
fvar< T > log(const fvar< T > &x)
Definition: log.hpp:12
scalar_seq_view provides a uniform sequence-like wrapper around either a scalar or a sequence of scal...
VectorBuilder< true, double, T_shape1, T_shape2 >::type beta_rng(const T_shape1 &alpha, const T_shape2 &beta, RNG &rng)
Return a Beta random variate with the supplied success and failure parameters using the given random ...
Definition: beta_rng.hpp:34
fvar< T > log_sum_exp(const std::vector< fvar< T > > &v)
Definition: log_sum_exp.hpp:12
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...
Definition: uniform_rng.hpp:36
fvar< T > beta(const fvar< T > &x1, const fvar< T > &x2)
Return fvar with the beta function applied to the specified arguments and its gradient.
Definition: beta.hpp:51
void check_positive_finite(const char *function, const char *name, const T_y &y)
Check if y is positive and finite.
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:11
size_t max_size(const T1 &x1, const T2 &x2)
Definition: max_size.hpp:9
VectorBuilder allocates type T1 values to be used as intermediate values.
void check_consistent_sizes(const char *function, const char *name1, const T1 &x1, const char *name2, const T2 &x2)
Check if the dimension of x1 is consistent with x2.

     [ Stan Home Page ] © 2011–2018, Stan Development Team.