Automatic Differentiation
 
Loading...
Searching...
No Matches
von_mises_cdf.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_PROB_VON_MISES_CDF_HPP
2#define STAN_MATH_PRIM_PROB_VON_MISES_CDF_HPP
3
6#include <cmath>
7
8namespace stan {
9namespace math {
10
11namespace internal {
12
22template <typename T_x, typename T_k>
23return_type_t<T_x, T_k> von_mises_cdf_series(const T_x& x, const T_k& k) {
24 const double pi = stan::math::pi();
25 int p = value_of_rec(28 + 0.5 * k - 100 / (k + 5) + 1);
26 auto s = sin(x);
27 auto c = cos(x);
28 auto sn = sin(p * x);
29 auto cn = cos(p * x);
30
31 using return_t = return_type_t<T_x, T_k>;
32 return_t R = 0.0;
33 return_t V = 0.0;
34
35 int n;
36 for (n = 1; n < p; n++) {
37 auto sn_tmp = sn * c - cn * s;
38 cn = cn * c + sn * s;
39 sn = sn_tmp;
40 R = 1 / (2.0 * (p - n) / k + R);
41 V = R * (sn / (p - n) + V);
42 }
43 return 0.5 + x / TWO_PI + V / pi;
44}
45
51template <typename T_x, typename T_k>
53 using std::exp;
54 using std::sqrt;
55
56 const auto b = sqrt(2 / pi()) * exp(k) / modified_bessel_first_kind(0, k);
57 const auto z = b * sin(x / 2);
58 const double mu = 0;
59 const double sigma = 1;
60
61 return normal_cdf(z, mu, sigma);
62}
63
70template <typename T_x, typename T_k>
72 using return_t = return_type_t<T_x, T_k>;
73 return_t f;
74 if (k < 49) {
75 f = von_mises_cdf_series(x, k);
76 if (f < 0) {
77 f = 0.0;
78 return f;
79 }
80 if (f > 1) {
81 f = 1.0;
82 return f;
83 }
84 return f;
85 } else if (k < 50) {
86 f = (50.0 - k) * von_mises_cdf_series(x, 49.0)
87 + (k - 49.0) * von_mises_cdf_normalapprox(x, 50.0);
88 return f;
89 } else {
91 return f;
92 }
93}
94
95} // namespace internal
96
117template <typename T_x, typename T_mu, typename T_k>
118inline return_type_t<T_x, T_mu, T_k> von_mises_cdf(const T_x& x, const T_mu& mu,
119 const T_k& k) {
121 const double pi = stan::math::pi();
122
123 using T_return = return_type_t<T_x, T_mu, T_k>;
124 using T_x_ref = ref_type_t<T_x>;
125 using T_mu_ref = ref_type_t<T_mu>;
126 using T_k_ref = ref_type_t<T_k>;
127 static char const* function = "von_mises_cdf";
128 check_consistent_sizes(function, "Random variable", x, "Location parameter",
129 mu, "Scale parameter", k);
130
131 T_x_ref x_ref = x;
132 T_mu_ref mu_ref = mu;
133 T_k_ref k_ref = k;
134
135 check_bounded(function, "Random variable", value_of(x_ref), -pi, pi);
136 check_finite(function, "Location parameter", value_of(mu_ref));
137 check_positive(function, "Scale parameter", value_of(k_ref));
138
139 if (size_zero(x, mu, k)) {
140 return 1.0;
141 }
142
143 T_return res(1.0);
144
145 scalar_seq_view<T_x_ref> x_vec(x_ref);
146 scalar_seq_view<T_mu_ref> mu_vec(mu_ref);
147 scalar_seq_view<T_k_ref> k_vec(k_ref);
148 size_t N = max_size(x, mu, k);
149
150 for (size_t n = 0; n < N; ++n) {
151 auto x_n = x_vec[n];
152 auto mu_n = mu_vec[n];
153 auto k_n = k_vec[n];
154
155 if (x_n == -pi) {
156 res *= 0.0;
157 } else if (x_n == pi) {
158 res *= 1.0;
159 } else {
160 // shift x so that mean is 0
161 T_return x2 = x_n - mu_n;
162
163 // x2 is on an interval (2*n*pi, (2*n + 1)*pi), move it to (-pi, pi)
164 x2 += pi;
165 const auto x_floor = floor(x2 / TWO_PI);
166 const auto x_modded = x2 - x_floor * TWO_PI;
167 x2 = x_modded - pi;
168
169 // mu is on an interval (2*n*pi, (2*n + 1)*pi), move it to (-pi, pi)
170 T_return mu2 = mu_n + pi;
171 const auto mu_floor = floor(mu2 / TWO_PI);
172 const auto mu_modded = mu2 - mu_floor * TWO_PI;
173 mu2 = mu_modded - pi;
174
175 // find cut
176 T_return cut, bound_val;
177 if (mu2 < 0) {
178 cut = mu2 + pi;
179 bound_val = -pi - mu2;
180 }
181 if (mu2 >= 0) {
182 cut = mu2 - pi;
183 bound_val = pi - mu2;
184 }
185
186 T_return f_bound_val = von_mises_cdf_centered(bound_val, k_n);
187 T_return cdf_n;
188 if (x_n <= cut) {
189 cdf_n = von_mises_cdf_centered(x2, k_n) - f_bound_val;
190 } else {
191 cdf_n = von_mises_cdf_centered(x2, k_n) + 1 - f_bound_val;
192 }
193
194 if (cdf_n < 0.0)
195 cdf_n = 0.0;
196
197 res *= cdf_n;
198 }
199 }
200
201 return res;
202}
203
204} // namespace math
205} // namespace stan
206
207#endif
scalar_seq_view provides a uniform sequence-like wrapper around either a scalar or a sequence of scal...
return_type_t< T_y_cl, T_loc_cl, T_scale_cl > normal_cdf(const T_y_cl &y, const T_loc_cl &mu, const T_scale_cl &sigma)
Returns the normal cumulative distribution function for the given location, and scale.
return_type_t< T_x, T_mu, T_k > von_mises_cdf(const T_x &x, const T_mu &mu, const T_k &k)
Calculates the cumulative distribution function of the von Mises distribution:
typename return_type< Ts... >::type return_type_t
Convenience type for the return type of the specified template parameters.
return_type_t< T_x, T_k > von_mises_cdf_normalapprox(const T_x &x, const T_k &k)
conv_mises_cdf_normapprox(x, k) is used to approximate the von Mises cdf for k > 50.
return_type_t< T_x, T_k > von_mises_cdf_series(const T_x &x, const T_k &k)
This implementation of the von Mises cdf is a copy of scipy's.
return_type_t< T_x, T_k > von_mises_cdf_centered(const T_x &x, const T_k &k)
This function calculates the cdf of the von Mises distribution in the case where the distribution has...
double value_of_rec(const fvar< T > &v)
Return the value of the specified variable.
fvar< T > sin(const fvar< T > &x)
Definition sin.hpp:14
bool size_zero(const T &x)
Returns 1 if input is of length 0, returns 0 otherwise.
Definition size_zero.hpp:19
void check_bounded(const char *function, const char *name, const T_y &y, const T_low &low, const T_high &high)
Check if the value is between the low and high values, inclusively.
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition value_of.hpp:18
fvar< T > modified_bessel_first_kind(int v, const fvar< T > &z)
void check_consistent_sizes(const char *)
Trivial no input case, this function is a no-op.
fvar< T > sqrt(const fvar< T > &x)
Definition sqrt.hpp:17
void check_finite(const char *function, const char *name, const T_y &y)
Return true if all values in y are finite.
fvar< T > floor(const fvar< T > &x)
Definition floor.hpp:12
fvar< T > cos(const fvar< T > &x)
Definition cos.hpp:14
static constexpr double TWO_PI
Twice the value of , .
Definition constants.hpp:62
void check_positive(const char *function, const char *name, const T_y &y)
Check if y is positive.
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
fvar< T > exp(const fvar< T > &x)
Definition exp.hpp:13
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 ...