von_mises_rng.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_SCAL_PROB_VON_MISES_RNG_HPP
2 #define STAN_MATH_PRIM_SCAL_PROB_VON_MISES_RNG_HPP
3 
13 #include <boost/random/uniform_real_distribution.hpp>
14 #include <boost/random/variate_generator.hpp>
15 
16 namespace stan {
17 namespace math {
18 
19 /**
20  * Return a von Mises random variate for the given location and concentration
21  * using the specified random number generator.
22  *
23  * mu and kappa can each be a scalar or a vector. Any non-scalar
24  * inputs must be the same length.
25  *
26  * The algorithm used in von_mises_rng is a modified version of the
27  * algorithm in:
28  *
29  * Efficient Simulation of the von Mises Distribution
30  * D. J. Best and N. I. Fisher
31  * Journal of the Royal Statistical Society. Series C (Applied Statistics),
32  * Vol. 28, No. 2 (1979), pp. 152-157
33  *
34  * See licenses/stan-license.txt for Stan license.
35  *
36  * @tparam T_loc Type of location parameter
37  * @tparam T_conc Type of concentration parameter
38  * @tparam RNG type of random number generator
39  * @param mu (Sequence of) location parameter(s)
40  * @param kappa (Sequence of) positive concentration parameter(s)
41  * @param rng random number generator
42  * @return (Sequence of) von Mises random variate(s)
43  * @throw std::domain_error if mu is infinite or kappa is nonpositive
44  * @throw std::invalid_argument if non-scalar arguments are of different
45  * sizes
46  */
47 template <typename T_loc, typename T_conc, class RNG>
49  const T_loc& mu, const T_conc& kappa, RNG& rng) {
50  using boost::random::uniform_real_distribution;
51  using boost::variate_generator;
52  static const char* function = "von_mises_rng";
53 
54  check_finite(function, "mean", mu);
55  check_positive_finite(function, "inverse of variance", kappa);
56  check_consistent_sizes(function, "Location parameter", mu,
57  "Concentration Parameter", kappa);
58 
59  scalar_seq_view<T_loc> mu_vec(mu);
60  scalar_seq_view<T_conc> kappa_vec(kappa);
61  size_t N = max_size(mu, kappa);
63 
64  variate_generator<RNG&, uniform_real_distribution<> > uniform_rng(
65  rng, uniform_real_distribution<>(0.0, 1.0));
66 
67  for (size_t n = 0; n < N; ++n) {
68  double r = 1 + std::pow((1 + 4 * kappa_vec[n] * kappa_vec[n]), 0.5);
69  double rho = 0.5 * (r - std::pow(2 * r, 0.5)) / kappa_vec[n];
70  double s = 0.5 * (1 + rho * rho) / rho;
71 
72  bool done = false;
73  double W;
74  while (!done) {
75  double Z = std::cos(pi() * uniform_rng());
76  W = (1 + s * Z) / (s + Z);
77  double Y = kappa_vec[n] * (s - W);
78  double U2 = uniform_rng();
79  done = Y * (2 - Y) - U2 > 0;
80 
81  if (!done)
82  done = std::log(Y / U2) + 1 - Y >= 0;
83  }
84 
85  double U3 = uniform_rng() - 0.5;
86  double sign = ((U3 >= 0) - (U3 <= 0));
87 
88  // it's really an fmod() with a positivity constraint
89  output[n]
90  = sign * std::acos(W)
91  + std::fmod(std::fmod(mu_vec[n], 2 * pi()) + 2 * stan::math::pi(),
92  2 * pi());
93  }
94 
95  return output.data();
96 }
97 
98 } // namespace math
99 } // namespace stan
100 #endif
ofstream output
T1 fmod(T1 number1, T2 number2)
Definition: d0nt_math.hpp:102
void check_finite(const char *function, const char *name, const T_y &y)
VectorBuilder< true, double, T_loc, T_conc >::type von_mises_rng(const T_loc &mu, const T_conc &kappa, RNG &rng)
constexpr T pow(T x)
Definition: pow.h:75
int sign(const T &z)
Definition: sign.hpp:9
T acos(T number)
Definition: d0nt_math.hpp:54
Float_t Y
Definition: plot.C:38
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
Float_t Z
Definition: plot.C:38
const XML_Char * s
Definition: expat.h:262
void check_positive_finite(const char *function, const char *name, const T_y &y)
size_t max_size(const T1 &x1, const T2 &x2)
Definition: max_size.hpp:9
T cos(T number)
Definition: d0nt_math.hpp:78
TRandom3 r(0)
double pi()
Definition: constants.hpp:82
#define W(x)
void check_consistent_sizes(const char *function, const char *name1, const T1 &x1, const char *name2, const T2 &x2)