von_mises_cdf.hpp

Go to the documentation of this file.

#ifndef STAN_MATH_PRIM_PROB_VON_MISES_CDF_HPP
#define STAN_MATH_PRIM_PROB_VON_MISES_CDF_HPP
 
#include <stan/math/prim/fun/constants.hpp>
#include <stan/math/prim/fun/modified_bessel_first_kind.hpp>
#include <cmath>
 
namespace stan {
namespace math {
 
namespace internal {
 
template <typename T_x, typename T_k>
inline return_type_t<T_x, T_k> von_mises_cdf_series(const T_x& x,
                                                    const T_k& k) {
  const double pi = stan::math::pi();
  int p = value_of_rec(28 + 0.5 * k - 100 / (k + 5) + 1);
  auto s = sin(x);
  auto c = cos(x);
  auto sn = sin(p * x);
  auto cn = cos(p * x);
 
  using return_t = return_type_t<T_x, T_k>;
  return_t R = 0.0;
  return_t V = 0.0;
 
  int n;
  for (n = 1; n < p; n++) {
    auto sn_tmp = sn * c - cn * s;
    cn = cn * c + sn * s;
    sn = sn_tmp;
    R = 1 / (2.0 * (p - n) / k + R);
    V = R * (sn / (p - n) + V);
  }
  return 0.5 + x / TWO_PI + V / pi;
}
 
template <typename T_x, typename T_k>
inline return_type_t<T_x, T_k> von_mises_cdf_normalapprox(const T_x& x,
                                                          const T_k& k) {
  using std::exp;
  using std::sqrt;
 
  const auto b = sqrt(2 / pi()) * exp(k) / modified_bessel_first_kind(0, k);
  const auto z = b * sin(x / 2);
  const double mu = 0;
  const double sigma = 1;
 
  return normal_cdf(z, mu, sigma);
}
 
template <typename T_x, typename T_k>
inline return_type_t<T_x, T_k> von_mises_cdf_centered(const T_x& x,
                                                      const T_k& k) {
  using return_t = return_type_t<T_x, T_k>;
  return_t f;
  if (k < 49) {
    f = von_mises_cdf_series(x, k);
    if (f < 0) {
      f = 0.0;
      return f;
    }
    if (f > 1) {
      f = 1.0;
      return f;
    }
    return f;
  } else if (k < 50) {
    f = (50.0 - k) * von_mises_cdf_series(x, 49.0)
        + (k - 49.0) * von_mises_cdf_normalapprox(x, 50.0);
    return f;
  } else {
    f = von_mises_cdf_normalapprox(x, k);
    return f;
  }
}
 
}  // namespace internal
 
template <typename T_x, typename T_mu, typename T_k>
inline return_type_t<T_x, T_mu, T_k> von_mises_cdf(const T_x& x, const T_mu& mu,
                                                   const T_k& k) {
  using internal::von_mises_cdf_centered;
  const double pi = stan::math::pi();
 
  using T_return = return_type_t<T_x, T_mu, T_k>;
  using T_x_ref = ref_type_t<T_x>;
  using T_mu_ref = ref_type_t<T_mu>;
  using T_k_ref = ref_type_t<T_k>;
  static char const* function = "von_mises_cdf";
  check_consistent_sizes(function, "Random variable", x, "Location parameter",
                         mu, "Scale parameter", k);
 
  T_x_ref x_ref = x;
  T_mu_ref mu_ref = mu;
  T_k_ref k_ref = k;
 
  check_bounded(function, "Random variable", value_of(x_ref), -pi, pi);
  check_finite(function, "Location parameter", value_of(mu_ref));
  check_positive(function, "Scale parameter", value_of(k_ref));
 
  if (size_zero(x, mu, k)) {
    return 1.0;
  }
 
  T_return res(1.0);
 
  scalar_seq_view<T_x_ref> x_vec(x_ref);
  scalar_seq_view<T_mu_ref> mu_vec(mu_ref);
  scalar_seq_view<T_k_ref> k_vec(k_ref);
  size_t N = max_size(x, mu, k);
 
  for (size_t n = 0; n < N; ++n) {
    auto x_n = x_vec[n];
    auto mu_n = mu_vec[n];
    auto k_n = k_vec[n];
 
    if (x_n == -pi) {
      res *= 0.0;
    } else if (x_n == pi) {
      res *= 1.0;
    } else {
      // shift x so that mean is 0
      T_return x2 = x_n - mu_n;
 
      // x2 is on an interval (2*n*pi, (2*n + 1)*pi), move it to (-pi, pi)
      x2 += pi;
      const auto x_floor = floor(x2 / TWO_PI);
      const auto x_modded = x2 - x_floor * TWO_PI;
      x2 = x_modded - pi;
 
      // mu is on an interval (2*n*pi, (2*n + 1)*pi), move it to (-pi, pi)
      T_return mu2 = mu_n + pi;
      const auto mu_floor = floor(mu2 / TWO_PI);
      const auto mu_modded = mu2 - mu_floor * TWO_PI;
      mu2 = mu_modded - pi;
 
      // find cut
      T_return cut, bound_val;
      if (mu2 < 0) {
        cut = mu2 + pi;
        bound_val = -pi - mu2;
      }
      if (mu2 >= 0) {
        cut = mu2 - pi;
        bound_val = pi - mu2;
      }
 
      T_return f_bound_val = von_mises_cdf_centered(bound_val, k_n);
      T_return cdf_n;
      if (x_n <= cut) {
        cdf_n = von_mises_cdf_centered(x2, k_n) - f_bound_val;
      } else {
        cdf_n = von_mises_cdf_centered(x2, k_n) + 1 - f_bound_val;
      }
 
      if (cdf_n < 0.0)
        cdf_n = 0.0;
 
      res *= cdf_n;
    }
  }
 
  return res;
}
 
}  // namespace math
}  // namespace stan
 
#endif