Automatic Differentiation
 
Loading...
Searching...
No Matches
von_mises_lcdf.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_PRIM_PROB_VON_MISES_LCDF_HPP
2#define STAN_MATH_PRIM_PROB_VON_MISES_LCDF_HPP
3
6#include <cmath>
7
8namespace stan {
9namespace math {
10
31template <typename T_x, typename T_mu, typename T_k>
33 const T_mu& mu,
34 const T_k& k) {
36 using std::log;
37 const double pi = stan::math::pi();
38
39 using T_return = return_type_t<T_x, T_mu, T_k>;
40 using T_x_ref = ref_type_t<T_x>;
41 using T_mu_ref = ref_type_t<T_mu>;
42 using T_k_ref = ref_type_t<T_k>;
43 static char const* function = "von_mises_lcdf";
44 check_consistent_sizes(function, "Random variable", x, "Location parameter",
45 mu, "Scale parameter", k);
46
47 T_x_ref x_ref = x;
48 T_mu_ref mu_ref = mu;
49 T_k_ref k_ref = k;
50
51 check_bounded(function, "Random variable", value_of(x_ref), -pi, pi);
52 check_finite(function, "Location parameter", value_of(mu_ref));
53 check_positive(function, "Scale parameter", value_of(k_ref));
54
55 if (size_zero(x, mu, k)) {
56 return 0.0;
57 }
58
59 T_return res(0.0);
60
61 scalar_seq_view<T_x_ref> x_vec(x_ref);
62 scalar_seq_view<T_mu_ref> mu_vec(mu_ref);
63 scalar_seq_view<T_k_ref> k_vec(k_ref);
64 size_t N = max_size(x, mu, k);
65
66 for (size_t n = 0; n < N; ++n) {
67 auto x_n = x_vec[n];
68 auto mu_n = mu_vec[n];
69 auto k_n = k_vec[n];
70
71 if (x_n == -pi) {
72 res += NEGATIVE_INFTY;
73 } else if (x_n == pi) {
74 res += 0.0;
75 } else {
76 // shift x so that mean is 0
77 T_return x2 = x_n - mu_n;
78
79 // x2 is on an interval (2*n*pi, (2*n + 1)*pi), move it to (-pi, pi)
80 x2 += pi;
81 const auto x_floor = floor(x2 / TWO_PI);
82 const auto x_modded = x2 - x_floor * TWO_PI;
83 x2 = x_modded - pi;
84
85 // mu is on an interval (2*n*pi, (2*n + 1)*pi), move it to (-pi, pi)
86 T_return mu2 = mu_n + pi;
87 const auto mu_floor = floor(mu2 / TWO_PI);
88 const auto mu_modded = mu2 - mu_floor * TWO_PI;
89 mu2 = mu_modded - pi;
90
91 // find cut
92 T_return cut, bound_val;
93 if (mu2 < 0) {
94 cut = mu2 + pi;
95 bound_val = -pi - mu2;
96 }
97 if (mu2 >= 0) {
98 cut = mu2 - pi;
99 bound_val = pi - mu2;
100 }
101
102 T_return f_bound_val = von_mises_cdf_centered(bound_val, k_n);
103 T_return cdf_n;
104 if (x_n <= cut) {
105 cdf_n = von_mises_cdf_centered(x2, k_n) - f_bound_val;
106 } else {
107 cdf_n = von_mises_cdf_centered(x2, k_n) + 1 - f_bound_val;
108 }
109
110 if (cdf_n < 0.0)
111 cdf_n = 0.0;
112
113 res += log(cdf_n);
114 }
115 }
116
117 return res;
118}
119
120} // namespace math
121} // namespace stan
122
123#endif
scalar_seq_view provides a uniform sequence-like wrapper around either a scalar or a sequence of scal...
return_type_t< T_x, T_mu, T_k > von_mises_lcdf(const T_x &x, const T_mu &mu, const T_k &k)
Calculates the log of 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_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...
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 > log(const fvar< T > &x)
Definition log.hpp:15
static constexpr double NEGATIVE_INFTY
Negative infinity.
Definition constants.hpp:51
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.
fvar< T > floor(const fvar< T > &x)
Definition floor.hpp:12
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
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 ...