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 
4 #include <boost/random/gamma_distribution.hpp>
5 #include <boost/random/uniform_real_distribution.hpp>
6 #include <boost/random/variate_generator.hpp>
13 
14 namespace stan {
15 namespace math {
16 
17 /**
18  * Return a Beta random variate with the supplied success and failure
19  * parameters using the given random number generator.
20  *
21  * alpha and beta can each be a scalar or a one-dimensional container. Any
22  * non-scalar inputs must be the same size.
23  *
24  * @tparam T_shape1 Type of success parameter
25  * @tparam T_shape2 Type of failure parameter
26  * @tparam RNG type of random number generator
27  * @param alpha (Sequence of) positive finite success parameter(s)
28  * @param beta (Sequence of) positive finite failure parameter(s)
29  * @param rng random number generator
30  * @return (Sequence of) beta random variate(s)
31  * @throw std::domain_error if alpha or beta are nonpositive
32  * @throw std::invalid_argument if non-scalar arguments are of different
33  * sizes
34  */
35 template <typename T_shape1, typename T_shape2, class RNG>
37  const T_shape1 &alpha, const T_shape2 &beta, RNG &rng) {
38  using boost::random::gamma_distribution;
39  using boost::random::uniform_real_distribution;
40  using boost::variate_generator;
41  static const char *function = "beta_rng";
42 
43  check_positive_finite(function, "First shape parameter", alpha);
44  check_positive_finite(function, "Second shape parameter", beta);
45  check_consistent_sizes(function, "First shape parameter", alpha,
46  "Second shape Parameter", beta);
47 
48  scalar_seq_view<T_shape1> alpha_vec(alpha);
49  scalar_seq_view<T_shape2> beta_vec(beta);
50  size_t N = max_size(alpha, beta);
52 
53  variate_generator<RNG &, uniform_real_distribution<>> uniform_rng(
54  rng, uniform_real_distribution<>(0.0, 1.0));
55  for (size_t n = 0; n < N; ++n) {
56  // If alpha and beta are large, trust the usual ratio of gammas
57  // method for generating beta random variables. If any parameter
58  // is small, work in log space and use Marsaglia and Tsang's trick
59  if (alpha_vec[n] > 1.0 && beta_vec[n] > 1.0) {
60  variate_generator<RNG &, gamma_distribution<>> rng_gamma_alpha(
61  rng, gamma_distribution<>(alpha_vec[n], 1.0));
62  variate_generator<RNG &, gamma_distribution<>> rng_gamma_beta(
63  rng, gamma_distribution<>(beta_vec[n], 1.0));
64  double a = rng_gamma_alpha();
65  double b = rng_gamma_beta();
66  output[n] = a / (a + b);
67  } else {
68  variate_generator<RNG &, gamma_distribution<>> rng_gamma_alpha(
69  rng, gamma_distribution<>(alpha_vec[n] + 1, 1.0));
70  variate_generator<RNG &, gamma_distribution<>> rng_gamma_beta(
71  rng, gamma_distribution<>(beta_vec[n] + 1, 1.0));
72  double log_a = std::log(uniform_rng()) / alpha_vec[n]
73  + std::log(rng_gamma_alpha());
74  double log_b
75  = std::log(uniform_rng()) / beta_vec[n] + std::log(rng_gamma_beta());
76  double log_sum = log_sum_exp(log_a, log_b);
77  output[n] = std::exp(log_a - log_sum);
78  }
79  }
80 
81  return output.data();
82 }
83 
84 } // namespace math
85 } // namespace stan
86 #endif
ofstream output
Double_t beta
VectorBuilder< true, double, T_shape1, T_shape2 >::type beta_rng(const T_shape1 &alpha, const T_shape2 &beta, RNG &rng)
Definition: beta_rng.hpp:36
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)
Definition: uniform_rng.hpp:37
void check_positive_finite(const char *function, const char *name, const T_y &y)
const double a
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
size_t max_size(const T1 &x1, const T2 &x2)
Definition: max_size.hpp:9
const hit & b
Definition: hits.cxx:21
void check_consistent_sizes(const char *function, const char *name1, const T1 &x1, const char *name2, const T2 &x2)