Automatic Differentiation
 
Loading...
Searching...
No Matches
beta_rng.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_PROB_BETA_RNG_HPP
2#define STAN_MATH_PRIM_PROB_BETA_RNG_HPP
3
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
35template <typename T_shape1, typename T_shape2, class RNG>
37 const T_shape1 &alpha, const T_shape2 &beta, RNG &rng) {
38 using boost::variate_generator;
39 using boost::random::gamma_distribution;
40 using boost::random::uniform_real_distribution;
41 using T_alpha_ref = ref_type_t<T_shape1>;
42 using T_beta_ref = ref_type_t<T_shape2>;
43 static constexpr const char *function = "beta_rng";
44 check_consistent_sizes(function, "First shape parameter", alpha,
45 "Second shape Parameter", beta);
46 T_alpha_ref alpha_ref = alpha;
47 T_beta_ref beta_ref = beta;
48 check_positive_finite(function, "First shape parameter", alpha_ref);
49 check_positive_finite(function, "Second shape parameter", beta_ref);
50
51 scalar_seq_view<T_alpha_ref> alpha_vec(alpha_ref);
52 scalar_seq_view<T_beta_ref> beta_vec(beta_ref);
53 size_t N = max_size(alpha, beta);
55
56 variate_generator<RNG &, uniform_real_distribution<>> uniform_rng(
57 rng, uniform_real_distribution<>(0.0, 1.0));
58 for (size_t n = 0; n < N; ++n) {
59 // If alpha and beta are large, trust the usual ratio of gammas
60 // method for generating beta random variables. If any parameter
61 // is small, work in log space and use Marsaglia and Tsang's trick
62 if (alpha_vec[n] > 1.0 && beta_vec[n] > 1.0) {
63 variate_generator<RNG &, gamma_distribution<>> rng_gamma_alpha(
64 rng, gamma_distribution<>(alpha_vec[n], 1.0));
65 variate_generator<RNG &, gamma_distribution<>> rng_gamma_beta(
66 rng, gamma_distribution<>(beta_vec[n], 1.0));
67 double a = rng_gamma_alpha();
68 double b = rng_gamma_beta();
69 output[n] = a / (a + b);
70 } else {
71 variate_generator<RNG &, gamma_distribution<>> rng_gamma_alpha(
72 rng, gamma_distribution<>(alpha_vec[n] + 1, 1.0));
73 variate_generator<RNG &, gamma_distribution<>> rng_gamma_beta(
74 rng, gamma_distribution<>(beta_vec[n] + 1, 1.0));
75 double log_a = std::log(uniform_rng()) / alpha_vec[n]
76 + std::log(rng_gamma_alpha());
77 double log_b
78 = std::log(uniform_rng()) / beta_vec[n] + std::log(rng_gamma_beta());
79 double log_sum = log_sum_exp(log_a, log_b);
80 output[n] = std::exp(log_a - log_sum);
81 }
82 }
83
84 return output.data();
85}
86
87} // namespace math
88} // namespace stan
89#endif
typename helper::type type
VectorBuilder allocates type T1 values to be used as intermediate values.
scalar_seq_view provides a uniform sequence-like wrapper around either a scalar or a sequence of scal...
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_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:36
void check_consistent_sizes(const char *)
Trivial no input case, this function is a no-op.
int64_t max_size(const T1 &x1, const Ts &... xs)
Calculate the size of the largest input.
Definition max_size.hpp:20
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 > log_sum_exp(const fvar< T > &x1, const fvar< T > &x2)
typename ref_type_if< true, T >::type ref_type_t
Definition ref_type.hpp:55
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...