Stan Math Library  2.20.0
reverse mode automatic differentiation
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 
11 #include <boost/random/uniform_real_distribution.hpp>
12 #include <boost/random/variate_generator.hpp>
13 
14 namespace stan {
15 namespace math {
16 
45 template <typename T_loc, typename T_conc, class RNG>
47  const T_loc& mu, const T_conc& kappa, RNG& rng) {
48  using boost::random::uniform_real_distribution;
49  using boost::variate_generator;
50  static const char* function = "von_mises_rng";
51 
52  check_finite(function, "mean", mu);
53  check_positive_finite(function, "inverse of variance", kappa);
54  check_consistent_sizes(function, "Location parameter", mu,
55  "Concentration Parameter", kappa);
56 
57  scalar_seq_view<T_loc> mu_vec(mu);
58  scalar_seq_view<T_conc> kappa_vec(kappa);
59  size_t N = max_size(mu, kappa);
61 
62  variate_generator<RNG&, uniform_real_distribution<> > uniform_rng(
63  rng, uniform_real_distribution<>(0.0, 1.0));
64 
65  for (size_t n = 0; n < N; ++n) {
66  double r = 1 + std::pow((1 + 4 * kappa_vec[n] * kappa_vec[n]), 0.5);
67  double rho = 0.5 * (r - std::pow(2 * r, 0.5)) / kappa_vec[n];
68  double s = 0.5 * (1 + rho * rho) / rho;
69 
70  bool done = false;
71  double W;
72  while (!done) {
73  double Z = std::cos(pi() * uniform_rng());
74  W = (1 + s * Z) / (s + Z);
75  double Y = kappa_vec[n] * (s - W);
76  double U2 = uniform_rng();
77  done = Y * (2 - Y) - U2 > 0;
78 
79  if (!done)
80  done = std::log(Y / U2) + 1 - Y >= 0;
81  }
82 
83  double U3 = uniform_rng() - 0.5;
84  double sign = ((U3 >= 0) - (U3 <= 0));
85 
86  // it's really an fmod() with a positivity constraint
87  output[n]
88  = sign * std::acos(W)
89  + std::fmod(std::fmod(mu_vec[n], 2 * pi()) + 2 * stan::math::pi(),
90  2 * pi());
91  }
92 
93  return output.data();
94 }
95 
96 } // namespace math
97 } // namespace stan
98 #endif
fvar< T > cos(const fvar< T > &x)
Definition: cos.hpp:12
void check_finite(const char *function, const char *name, const T_y &y)
Check if y is finite.
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...
int sign(const T &z)
Definition: sign.hpp:10
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...
fvar< T > fmod(const fvar< T > &x1, const fvar< T > &x2)
Definition: fmod.hpp:14
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
void check_positive_finite(const char *function, const char *name, const T_y &y)
Check if y is positive and finite.
size_t max_size(const T1 &x1, const T2 &x2)
Definition: max_size.hpp:9
fvar< T > acos(const fvar< T > &x)
Definition: acos.hpp:13
VectorBuilder allocates type T1 values to be used as intermediate values.
double pi()
Return the value of pi.
Definition: constants.hpp:80
fvar< T > pow(const fvar< T > &x1, const fvar< T > &x2)
Definition: pow.hpp:16
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.