Automatic Differentiation
 
Loading...
Searching...
No Matches
von_mises_rng.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_PROB_VON_MISES_RNG_HPP
2#define STAN_MATH_PRIM_PROB_VON_MISES_RNG_HPP
3
9#include <boost/random/uniform_real_distribution.hpp>
10#include <boost/random/variate_generator.hpp>
11#include <cmath>
12
13namespace stan {
14namespace math {
15
45template <typename T_loc, typename T_conc, class RNG>
47 const T_loc& mu, const T_conc& kappa, RNG& rng) {
48 using boost::variate_generator;
49 using boost::random::uniform_real_distribution;
50 using T_mu_ref = ref_type_t<T_loc>;
51 using T_kappa_ref = ref_type_t<T_conc>;
52 static constexpr const char* function = "von_mises_rng";
53 check_consistent_sizes(function, "Location parameter", mu, "Scale parameter",
54 kappa);
55 T_mu_ref mu_ref = mu;
56 T_kappa_ref kappa_ref = kappa;
57
58 check_finite(function, "Location parameter", mu_ref);
59 check_nonnegative(function, "Scale parameter", kappa_ref);
60 check_finite(function, "Scale parameter", kappa_ref);
61
62 scalar_seq_view<T_mu_ref> mu_vec(mu_ref);
63 scalar_seq_view<T_kappa_ref> kappa_vec(kappa_ref);
64 size_t N = max_size(mu, kappa_ref);
66
67 variate_generator<RNG&, uniform_real_distribution<> > uniform_rng(
68 rng, uniform_real_distribution<>(0.0, 1.0));
69
70 for (size_t n = 0; n < N; ++n) {
71 // for kappa sufficiently close to zero, it reduces to a
72 // circular uniform distribution centered at mu
73 if (kappa_vec[n] < 1.4e-8) {
74 output[n] = (uniform_rng() - 0.5) * TWO_PI
75 + std::fmod(std::fmod(mu_vec[n], TWO_PI) + TWO_PI, TWO_PI);
76 continue;
77 }
78
79 double r = 1 + std::pow((1 + 4 * kappa_vec[n] * kappa_vec[n]), 0.5);
80 double rho = 0.5 * (r - std::pow(2 * r, 0.5)) / kappa_vec[n];
81 double s = 0.5 * (1 + rho * rho) / rho;
82
83 bool done = false;
84 double W;
85 while (!done) {
86 double Z = std::cos(pi() * uniform_rng());
87 W = (1 + s * Z) / (s + Z);
88 double Y = kappa_vec[n] * (s - W);
89 double U2 = uniform_rng();
90 done = Y * (2 - Y) - U2 > 0;
91
92 if (!done) {
93 done = std::log(Y / U2) + 1 - Y >= 0;
94 }
95 }
96
97 double U3 = uniform_rng() - 0.5;
98 double sign = ((U3 >= 0) - (U3 <= 0));
99
100 // it's really an fmod() with a positivity constraint
101 output[n] = sign * std::acos(W)
102 + std::fmod(std::fmod(mu_vec[n], TWO_PI) + TWO_PI, TWO_PI);
103 }
104
105 return output.data();
106}
107
108} // namespace math
109} // namespace stan
110#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_loc, T_conc >::type von_mises_rng(const T_loc &mu, const T_conc &kappa, RNG &rng)
Return a von Mises random variate for the given location and concentration using the specified random...
void check_nonnegative(const char *function, const char *name, const T_y &y)
Check if y is non-negative.
auto sign(const T &x)
Returns signs of the arguments.
Definition sign.hpp:18
void check_consistent_sizes(const char *)
Trivial no input case, this function is a no-op.
void check_finite(const char *function, const char *name, const T_y &y)
Return true if all values in y are finite.
static constexpr double TWO_PI
Twice the value of , .
Definition constants.hpp:62
static constexpr double pi()
Return the value of pi.
Definition constants.hpp:36
int64_t max_size(const T1 &x1, const Ts &... xs)
Calculate the size of the largest input.
Definition max_size.hpp:20
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 ...