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