diff --git a/doc/graphs/von_mises/von_mises_ulps_cdf_4d.svg b/doc/graphs/von_mises/von_mises_ulps_cdf_4d.svg
new file mode 100644
index 0000000000..4b708fc043
--- /dev/null
+++ b/doc/graphs/von_mises/von_mises_ulps_cdf_4d.svg
@@ -0,0 +1,1048 @@
+
+
diff --git a/doc/graphs/von_mises/von_mises_ulps_cdf_4f.svg b/doc/graphs/von_mises/von_mises_ulps_cdf_4f.svg
new file mode 100644
index 0000000000..95ee0bfa45
--- /dev/null
+++ b/doc/graphs/von_mises/von_mises_ulps_cdf_4f.svg
@@ -0,0 +1,1049 @@
+
+
diff --git a/doc/graphs/von_mises/von_mises_ulps_cdf_64d.svg b/doc/graphs/von_mises/von_mises_ulps_cdf_64d.svg
new file mode 100644
index 0000000000..ceebe4fdcf
--- /dev/null
+++ b/doc/graphs/von_mises/von_mises_ulps_cdf_64d.svg
@@ -0,0 +1,1070 @@
+
+
diff --git a/doc/graphs/von_mises/von_mises_ulps_cdf_64f.svg b/doc/graphs/von_mises/von_mises_ulps_cdf_64f.svg
new file mode 100644
index 0000000000..910dfda78f
--- /dev/null
+++ b/doc/graphs/von_mises/von_mises_ulps_cdf_64f.svg
@@ -0,0 +1,1066 @@
+
+
diff --git a/doc/graphs/von_mises/von_mises_ulps_pdf_4d.svg b/doc/graphs/von_mises/von_mises_ulps_pdf_4d.svg
new file mode 100644
index 0000000000..70c7809e1b
--- /dev/null
+++ b/doc/graphs/von_mises/von_mises_ulps_pdf_4d.svg
@@ -0,0 +1,1050 @@
+
+
diff --git a/doc/graphs/von_mises/von_mises_ulps_pdf_4f.svg b/doc/graphs/von_mises/von_mises_ulps_pdf_4f.svg
new file mode 100644
index 0000000000..faec0d8b55
--- /dev/null
+++ b/doc/graphs/von_mises/von_mises_ulps_pdf_4f.svg
@@ -0,0 +1,1050 @@
+
+
diff --git a/doc/graphs/von_mises/von_mises_ulps_pdf_64d.svg b/doc/graphs/von_mises/von_mises_ulps_pdf_64d.svg
new file mode 100644
index 0000000000..6fd47e2174
--- /dev/null
+++ b/doc/graphs/von_mises/von_mises_ulps_pdf_64d.svg
@@ -0,0 +1,1050 @@
+
+
diff --git a/doc/graphs/von_mises/von_mises_ulps_pdf_64f.svg b/doc/graphs/von_mises/von_mises_ulps_pdf_64f.svg
new file mode 100644
index 0000000000..38e600c580
--- /dev/null
+++ b/doc/graphs/von_mises/von_mises_ulps_pdf_64f.svg
@@ -0,0 +1,1050 @@
+
+
diff --git a/include/boost/math/distributions/detail/common_error_handling.hpp b/include/boost/math/distributions/detail/common_error_handling.hpp
index 486fb0b5c8..8d866d50d5 100644
--- a/include/boost/math/distributions/detail/common_error_handling.hpp
+++ b/include/boost/math/distributions/detail/common_error_handling.hpp
@@ -11,6 +11,7 @@
#include
#include
+#include
// using boost::math::isfinite;
// using boost::math::isnan;
@@ -96,6 +97,25 @@ inline bool check_location(
return true;
}
+template
+inline bool check_angle(
+ const char* function,
+ RealType angle,
+ RealType* result,
+ const Policy& pol)
+{
+ if(!(boost::math::isfinite)(angle)
+ || angle < -boost::math::constants::pi()
+ || angle > +boost::math::constants::pi())
+ {
+ *result = policies::raise_domain_error(
+ function,
+ "Angle parameter is %1%, but must be between -pi and +pi!", angle, pol);
+ return false;
+ }
+ return true;
+}
+
template
inline bool check_x(
const char* function,
diff --git a/include/boost/math/distributions/fwd.hpp b/include/boost/math/distributions/fwd.hpp
index 5b212c8ec6..1d88b347ea 100644
--- a/include/boost/math/distributions/fwd.hpp
+++ b/include/boost/math/distributions/fwd.hpp
@@ -111,6 +111,9 @@ class triangular_distribution;
template
class uniform_distribution;
+template
+class von_mises_distribution;
+
template
class weibull_distribution;
@@ -148,6 +151,7 @@ class weibull_distribution;
typedef boost::math::students_t_distribution students_t;\
typedef boost::math::triangular_distribution triangular;\
typedef boost::math::uniform_distribution uniform;\
+ typedef boost::math::von_mises_distribution von_mises;\
typedef boost::math::weibull_distribution weibull;
#endif // BOOST_MATH_DISTRIBUTIONS_FWD_HPP
diff --git a/include/boost/math/distributions/von_mises.hpp b/include/boost/math/distributions/von_mises.hpp
new file mode 100644
index 0000000000..00c928543f
--- /dev/null
+++ b/include/boost/math/distributions/von_mises.hpp
@@ -0,0 +1,796 @@
+// Copyright John Maddock 2006, 2007.
+// Copyright Paul A. Bristow 2006, 2007.
+// Copyright Philipp C. J. Münster, 2020.
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifndef BOOST_MATH_DISTRIBUTIONS_VON_MISES_HPP
+#define BOOST_MATH_DISTRIBUTIONS_VON_MISES_HPP
+
+// https://en.wikipedia.org/wiki/Von_Mises_distribution
+// From MathWorld--A Wolfram Web Resource.
+// http://mathworld.wolfram.com/VonMisesDistribution.html
+
+#include
+#include
+#include
+#include // for besseli0 and besseli1
+#include // for erf
+#include
+
+#include
+
+namespace boost{ namespace math{
+
+template >
+class von_mises_distribution
+{
+public:
+ typedef RealType value_type;
+ typedef Policy policy_type;
+
+ von_mises_distribution(RealType l_mean = 0, RealType concentration = 1)
+ : m_mean(l_mean), m_concentration(concentration)
+ { // Default is a 'standard' von Mises distribution vM01.
+ static const char* function = "boost::math::von_mises_distribution<%1%>::von_mises_distribution";
+
+ RealType result;
+ detail::check_positive_x(function, concentration, &result, Policy());
+ detail::check_angle(function, l_mean, &result, Policy());
+ }
+
+ RealType mean() const
+ { // alias for location.
+ return m_mean;
+ }
+
+ RealType concentration() const
+ { // alias for scale.
+ return m_concentration;
+ }
+
+ // Synonyms, provided to allow generic use of find_location and find_scale.
+ RealType location() const
+ { // location.
+ return m_mean;
+ }
+ RealType scale() const
+ { // scale.
+ return m_concentration;
+ }
+
+private:
+ //
+ // Data members:
+ //
+ RealType m_mean; // distribution mean or location.
+ RealType m_concentration; // distribution standard deviation or scale.
+}; // class von_mises_distribution
+
+typedef von_mises_distribution von_mises;
+
+#ifdef BOOST_MSVC
+#pragma warning(push)
+#pragma warning(disable:4127)
+#endif
+
+template
+inline const std::pair range(const von_mises_distribution& /*dist*/)
+{ // Range of permissible values for random variable x.
+ if (std::numeric_limits::has_infinity)
+ {
+ return std::pair(-std::numeric_limits::infinity(),
+ +std::numeric_limits::infinity()); // - to + infinity.
+ }
+ else
+ { // Can only use max_value.
+ using boost::math::tools::max_value;
+ return std::pair(-max_value(),
+ +max_value()); // - to + max value.
+ }
+}
+
+template
+inline const std::pair support(const von_mises_distribution& dist)
+{ // This is range values for random variable x where cdf rises from 0 to 1, and outside it, the pdf is zero.
+ return std::pair(dist.mean() - constants::pi(),
+ dist.mean() + constants::pi()); // [µ-π, µ+π)
+}
+
+#ifdef BOOST_MSVC
+#pragma warning(pop)
+#endif
+
+namespace detail {
+// float version of pdf_impl
+template
+inline RealType pdf_impl(const von_mises_distribution& dist, RealType x,
+ const boost::integral_constant&)
+{
+ RealType mean = dist.mean();
+ RealType conc = dist.concentration();
+
+ BOOST_MATH_STD_USING
+ if(conc < 87)
+ {
+ RealType bessel_i0 = cyl_bessel_i(0, conc, Policy());
+ return exp(conc * cos(x - mean)) / bessel_i0 / boost::math::constants::two_pi();
+ }
+ else // exp(88) > MAX_FLOAT
+ {
+ // we make use of I0(conc) = exp(conc) * P(conc), with polynomial P
+ // polynomial coefficients from boost/math/special_functions/detail/bessel_i0.hpp
+ static const float P[] = {
+ 3.98942280401432677e-01f,
+ 4.98677850501790847e-02f,
+ 2.80506290907257351e-02f,
+ 2.92194053028393074e-02f,
+ 4.47422143699726895e-02f,
+ };
+ RealType result = exp(conc * (cos(x - mean) - 1.f));
+ result /= boost::math::tools::evaluate_polynomial(P, RealType(1.f / conc)) / sqrt(conc)
+ * boost::math::constants::two_pi();
+ return result;
+ }
+}
+
+// double version of pdf_impl
+template
+inline RealType pdf_impl(const von_mises_distribution& dist, RealType x,
+ const boost::integral_constant&)
+{
+ RealType mean = dist.mean();
+ RealType conc = dist.concentration();
+
+ BOOST_MATH_STD_USING
+ if(conc < 709)
+ {
+ RealType bessel_i0 = cyl_bessel_i(0, conc, Policy());
+ return exp(conc * cos(x - mean)) / bessel_i0 / boost::math::constants::two_pi();
+ }
+ else // exp(709) > MAX_DOUBLE
+ {
+ // we make use of I0(conc) = exp(conc) * P(conc), with polynomial P
+ // polynomial coefficients from boost/math/special_functions/detail/bessel_i0.hpp
+ static const double P[] = {
+ 3.98942280401432905e-01,
+ 4.98677850491434560e-02,
+ 2.80506308916506102e-02,
+ 2.92179096853915176e-02,
+ 4.53371208762579442e-02
+ };
+ RealType result = exp(conc * (cos(x - mean) - 1.0));
+ result /= boost::math::tools::evaluate_polynomial(P, RealType(1.0 / conc)) / sqrt(conc)
+ * boost::math::constants::two_pi();
+ return result;
+ }
+}
+
+// long double version of pdf_impl
+template
+inline RealType pdf_impl(const von_mises_distribution& dist, RealType x,
+ const boost::integral_constant&)
+{
+ RealType mean = dist.mean();
+ RealType conc = dist.concentration();
+
+ BOOST_MATH_STD_USING
+ if (conc < 1000)
+ {
+ RealType bessel_i0 = cyl_bessel_i(0, conc, Policy());
+ return exp(conc * cos(x - mean)) / bessel_i0 / boost::math::constants::two_pi();
+ }
+ else
+ {
+ // Bessel I0 over[50, INF]
+ // Max error in interpolated form : 5.587e-20
+ // Max Error found at float80 precision = Poly : 8.776852e-20
+ static const RealType P[] = {
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +3.98942280401432677955074061e-01),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +4.98677850501789875615574058e-02),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +2.80506290908675604202206833e-02),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +2.92194052159035901631494784e-02),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +4.47422430732256364094681137e-02),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +9.05971614435738691235525172e-02),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +2.29180522595459823234266708e-01),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +6.15122547776140254569073131e-01),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +7.48491812136365376477357324e+00),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, -2.45569740166506688169730713e+02),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +9.66857566379480730407063170e+03),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, -2.71924083955641197750323901e+05),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +5.74276685704579268845870586e+06),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, -8.89753803265734681907148778e+07),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +9.82590905134996782086242180e+08),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, -7.30623197145529889358596301e+09),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +3.27310000726207055200805893e+10),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, -6.64365417189215599168817064e+10)
+ };
+ RealType result = exp(conc * (cos(x - mean) - 1.0));
+ result /= boost::math::tools::evaluate_polynomial(P, RealType(1.0 / conc)) / sqrt(conc)
+ * boost::math::constants::two_pi();
+ return result;
+ }
+}
+template
+inline RealType pdf_impl(const von_mises_distribution& dist, RealType x,
+ const boost::integral_constant&)
+{
+ BOOST_MATH_STD_USING // for ADL of std functions
+ RealType mean = dist.mean();
+ RealType conc = dist.concentration();
+ RealType bessel_i0 = cyl_bessel_i(0, conc, Policy());
+ return exp(conc * cos(x - mean)) / bessel_i0 / boost::math::constants::two_pi();
+}
+} // namespace detail
+
+template
+inline RealType pdf(const von_mises_distribution& dist, const RealType& x)
+{
+ BOOST_MATH_STD_USING // for ADL of std functions
+
+ RealType conc = dist.concentration();
+ RealType mean = dist.mean();
+
+ static const char* function = "boost::math::pdf(const von_mises_distribution<%1%>&, %1%)";
+
+ RealType result = 0;
+ if (false == detail::check_positive_x(function, conc, &result, Policy()))
+ return result;
+ if (false == detail::check_angle(function, mean, &result, Policy()))
+ return result;
+ if (false == detail::check_angle(function, x - mean, &result, Policy()))
+ return result;
+ // Below produces MSVC 4127 warnings, so the above used instead.
+ //if(std::numeric_limits::has_infinity && abs(x) == std::numeric_limits::infinity())
+ //{ // pdf + and - infinity is zero.
+ // return 0;
+ //}
+ typedef boost::integral_constant::digits == 0) || (std::numeric_limits::radix != 2)) ?
+ 0 :
+ std::numeric_limits::digits <= 24 ?
+ 24 :
+ std::numeric_limits::digits <= 53 ?
+ 53 :
+ std::numeric_limits::digits <= 64 ?
+ 64 :
+ std::numeric_limits::digits <= 113 ?
+ 113 : -1
+ > tag_type;
+
+ return detail::pdf_impl(dist, x, tag_type());
+} // pdf
+
+namespace detail {
+
+template
+inline RealType cdf_impl(const von_mises_distribution& dist, const RealType& x)
+{
+ BOOST_MATH_STD_USING // for ADL of std functions
+
+ RealType conc = dist.concentration();
+ RealType mean = dist.mean();
+
+ RealType u = x - mean;
+ if (u <= -boost::math::constants::pi())
+ return 0;
+ if (u >= +boost::math::constants::pi())
+ return 1;
+
+ // We use the Fortran algorithm designed by Geoffrey W. Hill in
+ // "Algorithm 518: Incomplete Bessel Function I0. The Von Mises Distribution", 1977, ACM
+ // DOI: 10.1145/355744.355753
+ RealType result = 0;
+
+ int digits = std::numeric_limits::max_digits10 - 1;
+ RealType ck = ((0.1611*digits - 2.8778)*digits + 18.45)*digits - 35.4;
+ if (conc > ck) {
+ RealType c = 24.0 * conc;
+ RealType v = c - 56;
+ RealType r = sqrt((54.0 / (347.0 / v + 26.0 - c) - 6.0 + c) / 12.0);
+ RealType z = sin(u / 2.0) * r;
+ RealType s = z * z * 2;
+ v = v - s + 3;
+ RealType y = (c - s - s - 16.0) / 3.0;
+ y = ((s + 1.75) * s + 83.5) / v - y;
+ result = boost::math::erf(z - s / (y * y) * z) / 2 + 0.5;
+ }
+ else {
+ RealType v = 0;
+ if(conc > 0) {
+ // extrapolation of the tables given in the paper
+ RealType a1 = (0.33 * digits - 2.6666) * digits + 12;
+ RealType a2 = std::max(0.5, std::min(1.5 - digits / 12, 1.0));
+ RealType a3 = 8;//digits <= 6 ? 3 : (1 << (digits - 5));
+ RealType a4 = digits <= 6 ? 1 : std::pow(1.5, digits - 8);
+
+ int iterations = static_cast(ceil(a1 + conc * a2 - a3 / (conc + a4)));
+ RealType r = 0;
+ RealType z = 2 / conc;
+ for (int j = iterations - 1; j > 0; --j) {
+ RealType sj = sin(j * u);
+ r = 1 / (j * z + r);
+ v = (sj / j + v) * r;
+ }
+ }
+ result = (x - mean + boost::math::constants::pi()) / 2;
+ result = (result + v) / boost::math::constants::pi();
+ }
+
+ return result <= 0 ? 0 : (1 <= result ? 1 : result);
+}
+} // namespace detail
+
+template
+inline RealType cdf(const von_mises_distribution& dist, const RealType& x)
+{
+ RealType conc = dist.concentration();
+ RealType mean = dist.mean();
+ static const char* function = "boost::math::cdf(const von_mises_distribution<%1%>&, %1%)";
+ RealType result = 0;
+ if (false == detail::check_positive_x(function, conc, &result, Policy()))
+ return result;
+ if (false == detail::check_angle(function, mean, &result, Policy()))
+ return result;
+ if (false == detail::check_angle(function, x - mean, &result, Policy()))
+ return result;
+
+ return detail::cdf_impl(dist, x);
+} // cdf
+
+template
+inline RealType quantile(const von_mises_distribution& dist, const RealType& p)
+{
+ BOOST_MATH_STD_USING // for ADL of std functions
+
+ RealType conc = dist.concentration();
+ RealType mean = dist.mean();
+ static const char* function
+ = "boost::math::quantile(const von_mises_distribution<%1%>&, %1%)";
+
+ RealType result = 0;
+ if (false == detail::check_positive_x(function, conc, &result, Policy()))
+ return result;
+ if (false == detail::check_angle(function, mean, &result, Policy()))
+ return result;
+ if (false == detail::check_probability(function, p, &result, Policy()))
+ return result;
+
+ if (p <= 0)
+ return -boost::math::constants::pi();
+ if (p >= 1)
+ return +boost::math::constants::pi();
+
+ typedef boost::integral_constant::digits == 0)
+ || (std::numeric_limits::radix != 2)) ? 0 :
+ std::numeric_limits::digits <= 24 ? 24 :
+ std::numeric_limits::digits <= 53 ? 53 :
+ std::numeric_limits::digits <= 64 ? 64 :
+ std::numeric_limits::digits <= 113 ? 113 :
+ -1
+ > tag_type;
+
+ struct step_func {
+ const von_mises_distribution& dist;
+ const RealType p;
+ std::pair operator()(RealType x) {
+ return std::make_pair(detail::cdf_impl(dist, x) - p, // f(x)
+ detail::pdf_impl(dist, x, tag_type())); // f'(x)
+ }
+ };
+ RealType lower = mean - boost::math::constants::pi();
+ RealType upper = mean + boost::math::constants::pi();
+ RealType zero = boost::math::tools::newton_raphson_iterate(
+ step_func{dist, p}, mean, lower, upper, 15 /* digits */);
+
+ return zero;
+} // quantile
+
+template
+inline RealType cdf(const complemented2_type, RealType>& c)
+{
+ RealType conc = c.dist.concentration();
+ RealType mean = c.dist.mean();
+ RealType x = c.param;
+ static const char* function
+ = "boost::math::cdf(const complement(von_mises_distribution<%1%>&), %1%)";
+
+ RealType result = 0;
+ if (false == detail::check_positive_x(function, conc, &result, Policy()))
+ return result;
+ if (false == detail::check_angle(function, mean, &result, Policy()))
+ return result;
+ if (false == detail::check_angle(function, x - mean, &result, Policy()))
+ return result;
+
+ return detail::cdf_impl(c.dist, 2 * mean - x);
+} // cdf complement
+
+template
+inline RealType quantile(const complemented2_type, RealType>& c)
+{
+ BOOST_MATH_STD_USING // for ADL of std functions
+
+ RealType conc = c.dist.concentration();
+ RealType mean = c.dist.mean();
+ static const char* function
+ = "boost::math::quantile(const complement(von_mises_distribution<%1%>&), %1%)";
+
+ RealType result = 0;
+ if (false == detail::check_positive_x(function, conc, &result, Policy()))
+ return result;
+ if (false == detail::check_angle(function, mean, &result, Policy()))
+ return result;
+ RealType q = c.param;
+ if (false == detail::check_probability(function, q, &result, Policy()))
+ return result;
+
+ if (q <= 0)
+ return +boost::math::constants::pi();
+ if (q >= 1)
+ return -boost::math::constants::pi();
+
+ typedef boost::integral_constant::digits == 0)
+ || (std::numeric_limits::radix != 2)) ? 0 :
+ std::numeric_limits::digits <= 24 ? 24 :
+ std::numeric_limits::digits <= 53 ? 53 :
+ std::numeric_limits::digits <= 64 ? 64 :
+ std::numeric_limits::digits <= 113 ? 113 :
+ -1
+ > tag_type;
+
+ struct step_func {
+ const complemented2_type, RealType>& c;
+ std::pair operator()(RealType x) {
+ RealType xc = 2 * c.dist.mean() - x;
+ return std::make_pair(detail::cdf_impl(c.dist, xc) - c.param, // f(x)
+ -detail::pdf_impl(c.dist, xc, tag_type())); // f'(x)
+ }
+ };
+ RealType lower = mean - boost::math::constants::pi();
+ RealType upper = mean + boost::math::constants::pi();
+ RealType zero = boost::math::tools::newton_raphson_iterate(
+ step_func{c}, mean, lower, upper, 15 /* digits */);
+
+ return zero;
+} // quantile
+
+template
+inline RealType mean(const von_mises_distribution& dist)
+{
+ return dist.mean();
+}
+
+template
+inline RealType standard_deviation(const von_mises_distribution& dist)
+{
+ BOOST_MATH_STD_USING
+ RealType bessel_quot = cyl_bessel_i(1, dist.concentration(), Policy())
+ / cyl_bessel_i(0, dist.concentration(), Policy());
+ return sqrt(-2 * log(bessel_quot));
+}
+
+template
+inline RealType mode(const von_mises_distribution& dist)
+{
+ return dist.mean();
+}
+
+template
+inline RealType median(const von_mises_distribution& dist)
+{
+ return dist.mean();
+}
+
+namespace detail {
+// float version of variance_impl
+template
+inline RealType variance_impl(const von_mises_distribution& dist,
+ const boost::integral_constant&)
+{
+ RealType conc = dist.concentration();
+ BOOST_MATH_STD_USING
+ if(conc < 7.75)
+ {
+ RealType bessel_i0 = cyl_bessel_i(0, conc, Policy());
+ RealType bessel_i1 = cyl_bessel_i(1, conc, Policy());
+ return 1 - bessel_i1 / bessel_i0;
+ }
+ else if(conc < 50)
+ {
+ // Polynomial coefficients from
+ // boost/math/special_functions/detail/bessel_i0.hpp and
+ // boost/math/special_functions/detail/bessel_i1.hpp
+ // compute numerator as I0(conc) - I1(conc) from single polynomial
+ static const float P_numer[] = {
+ +5.356107887570000000e-07f,
+ +1.994139882543095464e-01f,
+ +7.683426463016022940e-02f,
+ +4.007722563185265850e-02f,
+ +2.785578524715388070e-01f,
+ };
+ static const float P_denom[] = {
+ +3.98942651588301770e-01f,
+ +4.98327234176892844e-02f,
+ +2.91866904423115499e-02f,
+ +1.35614940793742178e-02f,
+ +1.31409251787866793e-01f,
+ };
+
+ RealType x = 1 / conc;
+ RealType numer = boost::math::tools::evaluate_polynomial(P_numer, x);
+ RealType denom = boost::math::tools::evaluate_polynomial(P_denom, x);
+
+ return numer / denom;
+ }
+ else
+ {
+ static const float P_numer[] = {
+ 1.644239196640000000e-07f,
+ 1.994490498867993467e-01f,
+ 7.569820327857441460e-02f,
+ 5.573513685531774810e-02f,
+ 1.918908150536447035e-01f
+ };
+ static const float P_denom[] = {
+ 3.98942280401432677e-01f,
+ 4.98677850501790847e-02f,
+ 2.80506290907257351e-02f,
+ 2.92194053028393074e-02f,
+ 4.47422143699726895e-02f,
+ };
+ RealType x = 1 / conc;
+ RealType numer = boost::math::tools::evaluate_polynomial(P_numer, x);
+ RealType denom = boost::math::tools::evaluate_polynomial(P_denom, x);
+
+ return numer / denom;
+ }
+}
+
+// double version of variance_impl
+template
+inline RealType variance_impl(const von_mises_distribution& dist,
+ const boost::integral_constant&)
+{
+ RealType conc = dist.concentration();
+ BOOST_MATH_STD_USING
+ if(conc < 7.75)
+ {
+ RealType bessel_i0 = cyl_bessel_i(0, conc, Policy());
+ RealType bessel_i1 = cyl_bessel_i(1, conc, Policy());
+ return 1 - bessel_i1 / bessel_i0;
+ }
+ else if(conc < 500)
+ {
+ // Polynomial coefficients from
+ // boost/math/special_functions/detail/bessel_i0.hpp and
+ // boost/math/special_functions/detail/bessel_i1.hpp
+ // compute numerator as I0(conc) - I1(conc) from single polynomial
+ static const double P_numer[] = {
+ -1.55174e-14,
+ +1.994711402218073518e-01,
+ +7.480166592881663552e-02,
+ +7.013008203242116521e-02,
+ +1.016110940936680100e-02,
+ +2.837895300433059925e-01,
+ -6.808807273294442296e+00
+ +4.756438487430168291e+02
+ -2.312449372965164219e+04,
+ +8.645232325622160644e+05,
+ -2.509248065298433807e+07,
+ +5.705709779789331679e+08,
+ -1.021122689535998576e+10,
+ +1.439407686911892518e+11,
+ -1.593043530624297061e+12,
+ +1.373585989454675265e+13,
+ -9.104389306577963414e+13,
+ +4.540402039126762439e+14,
+ -1.645981875199121119e+15,
+ +4.091196019695783875e+15,
+ -6.233158807315033853e+15,
+ +4.389193640817412685e+15,
+ };
+ static const double P_denom[] = {
+ 3.98942280401425088e-01,
+ 4.98677850604961985e-02,
+ 2.80506233928312623e-02,
+ 2.92211225166047873e-02,
+ 4.44207299493659561e-02,
+ 1.30970574605856719e-01,
+ -3.35052280231727022e+00,
+ 2.33025711583514727e+02,
+ -1.13366350697172355e+04,
+ 4.24057674317867331e+05,
+ -1.23157028595698731e+07,
+ 2.80231938155267516e+08,
+ -5.01883999713777929e+09,
+ 7.08029243015109113e+10,
+ -7.84261082124811106e+11,
+ 6.76825737854096565e+12,
+ -4.49034849696138065e+13,
+ 2.24155239966958995e+14,
+ -8.13426467865659318e+14,
+ 2.02391097391687777e+15,
+ -3.08675715295370878e+15,
+ 2.17587543863819074e+15
+ };
+ RealType x = 1 / conc;
+ RealType numer = boost::math::tools::evaluate_polynomial(P_numer, x);
+ RealType denom = boost::math::tools::evaluate_polynomial(P_denom, x);
+
+ return numer / denom;
+ }
+ else
+ {
+ static const double P_numer[] = {
+ +1.423e-15,
+ +1.994711401959018717e-01,
+ +7.480168411736836931e-02,
+ +7.012212565916144652e-02,
+ +1.037734243240472200e-01,
+ };
+ static const double P_denom[] = {
+ 3.98942280401432905e-01,
+ 4.98677850491434560e-02,
+ 2.80506308916506102e-02,
+ 2.92179096853915176e-02,
+ 4.53371208762579442e-02,
+ };
+ RealType x = 1 / conc;
+ RealType numer = boost::math::tools::evaluate_polynomial(P_numer, x);
+ RealType denom = boost::math::tools::evaluate_polynomial(P_denom, x);
+
+ return numer / denom;
+ }
+}
+
+// long double version of variance_impl
+template
+inline RealType variance_impl(const von_mises_distribution& dist,
+ const boost::integral_constant&)
+{
+ BOOST_MATH_STD_USING
+ RealType conc = dist.concentration();
+ if (conc < 50)
+ {
+ RealType bessel_i0 = cyl_bessel_i(0, conc, Policy());
+ RealType bessel_i1 = cyl_bessel_i(1, conc, Policy());
+ return 1 - bessel_i1 / bessel_i0;
+ }
+ else
+ {
+ static const RealType P[] = {
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +3.98942280401432677955074061e-01),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +4.98677850501789875615574058e-02),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +2.80506290908675604202206833e-02),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +2.92194052159035901631494784e-02),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +4.47422430732256364094681137e-02),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +9.05971614435738691235525172e-02),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +2.29180522595459823234266708e-01),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +6.15122547776140254569073131e-01),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +7.48491812136365376477357324e+00),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, -2.45569740166506688169730713e+02),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +9.66857566379480730407063170e+03),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, -2.71924083955641197750323901e+05),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +5.74276685704579268845870586e+06),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, -8.89753803265734681907148778e+07),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +9.82590905134996782086242180e+08),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, -7.30623197145529889358596301e+09),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, +3.27310000726207055200805893e+10),
+ BOOST_MATH_BIG_CONSTANT(RealType, 64, -6.64365417189215599168817064e+10)
+ };
+ return 0;
+ }
+}
+
+// quad version of variance_impl
+template
+inline RealType variance_impl(const von_mises_distribution& dist,
+ const boost::integral_constant&)
+{
+ BOOST_MATH_STD_USING
+ RealType conc = dist.concentration();
+ if (conc < 100)
+ {
+ RealType bessel_i0 = cyl_bessel_i(0, conc, Policy());
+ RealType bessel_i1 = cyl_bessel_i(1, conc, Policy());
+ return 1 - bessel_i1 / bessel_i0;
+ }
+ else
+ {
+ // Polynomial for Bessel I0 / exp(conc) / sqrt(conc)
+ static const RealType P_denom[] = {
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, 3.9894228040143267793994605993438166526772e-01),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, 4.9867785050179084742493257495245185241487e-02),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, 2.8050629090725735167652437695397756897920e-02),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, 2.9219405302839307466358297347675795965363e-02),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, 4.4742214369972689474366968442268908028204e-02),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, 9.0602984099194778006610058410222616383078e-02),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, 2.2839502241666629677015839125593079416327e-01),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, 6.8926354981801627920292655818232972385750e-01),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, 2.4231921590621824187100989532173995000655e+00),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, 9.7264260959693775207585700654645245723497e+00),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, 4.3890136225398811195878046856373030127018e+01),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, 2.1999720924619285464910452647408431234369e+02),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, 1.2076909538525038580501368530598517194748e+03),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, 7.5684635141332367730007149159063086133399e+03),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, 3.5178192543258299267923025833141286569141e+04),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, 6.2966297919851965784482163987240461837728e+05)
+ };
+
+ // Polynomial for (Bessel I0 - Bessel I1) / exp(conc) / sqrt(x)
+ static const RealType P_numer[] = {
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, -3.368165e-34),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, +0.1994711402007163389699730299733589146073),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, +0.07480167757526862711373984952175456888896),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, +0.07012657272681433791923412617430580227103),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, +0.10226791855993757596915805134093796837369),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, +0.2013399646648772644282448264813045181469),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, +0.4983164125454637874163933874417571110056),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, +1.4845676457582822590839158984567308297366),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, +5.169476606935535670414552625141409318434),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, +20.59713743665130419014001937211862931161),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, +92.40031163861578044111910646492625263225),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, +460.9440321063085921197536056693132597443),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, +2.520575547528944544570101030955801999019e+03),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, +1.5734053646329490893326466550032992462076e+04),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, +7.319778356894459435808347175389511056370e+04),
+ BOOST_MATH_BIG_CONSTANT(RealType, 113, +1.2997438696903014448182029282439919466883e+06),
+ };
+ RealType x = 1 / conc;
+ RealType numer = boost::math::tools::evaluate_polynomial(P_numer, x);
+ RealType denom = boost::math::tools::evaluate_polynomial(P_denom, x);
+
+ return numer / denom;
+ }
+}
+} // namespace detail
+
+template
+inline RealType variance(const von_mises_distribution& dist)
+{
+
+ typedef boost::integral_constant::digits == 0)
+ || (std::numeric_limits::radix != 2)) ? 0 :
+ std::numeric_limits::digits <= 24 ? 24 :
+ std::numeric_limits::digits <= 53 ? 53 :
+ std::numeric_limits::digits <= 64 ? 64 :
+ std::numeric_limits::digits <= 113 ? 113 :
+ -1
+ > tag_type;
+
+ return detail::variance_impl(dist, tag_type());
+}
+
+template
+inline RealType skewness(const von_mises_distribution& /*dist*/)
+{
+ return 0;
+}
+
+template
+inline RealType entropy(const von_mises_distribution & dist)
+{
+ BOOST_MATH_STD_USING
+ RealType arg = constants::two_pi() * cyl_bessel_i(0, dist.concentration(), Policy());
+ RealType bessel_quot = cyl_bessel_i(1, dist.concentration(), Policy())
+ / cyl_bessel_i(0, dist.concentration(), Policy());
+ return log(arg) - dist.concentration() * bessel_quot;
+}
+
+} // namespace math
+} // namespace boost
+
+// This include must be at the end, *after* the accessors
+// for this distribution have been defined, in order to
+// keep compilers that support two-phase lookup happy.
+#include
+
+#endif // BOOST_MATH_DISTRIBUTIONS_VON_MISES_HPP
+
+
diff --git a/reporting/accuracy/von_mises_ulps.cpp b/reporting/accuracy/von_mises_ulps.cpp
new file mode 100644
index 0000000000..a63b494bbf
--- /dev/null
+++ b/reporting/accuracy/von_mises_ulps.cpp
@@ -0,0 +1,81 @@
+#include
+
+#include
+#include
+
+using precise_real = long double;
+
+template
+void generate_ulps_plot_pdf(CoarseReal concentration)
+{
+ auto hi_conc = static_cast(concentration);
+ auto high_precision = [hi_conc](precise_real x)
+ {
+ boost::math::von_mises_distribution dist(0, hi_conc);
+ return boost::math::pdf(dist, x);
+ };
+
+ auto low_precision = [concentration](CoarseReal x)
+ {
+ boost::math::von_mises_distribution dist(0, concentration);
+ return boost::math::pdf(dist, x);
+ };
+
+ using hp_func = decltype (high_precision);
+
+ boost::math::tools::ulps_plot plot(
+ high_precision, 0, boost::math::constants::pi());
+ plot.add_fn(low_precision);
+ std::string filename = "von_mises_ulps_pdf_"
+ + std::to_string(static_cast(concentration))
+ + typeid(CoarseReal).name()
+ + ".svg";
+ plot.write(filename);
+}
+
+template
+void generate_ulps_plot_cdf(CoarseReal concentration)
+{
+ auto hi_conc = static_cast(concentration);
+ auto high_precision = [hi_conc](precise_real x)
+ {
+ boost::math::von_mises_distribution dist(0, hi_conc);
+ return boost::math::cdf(dist, x);
+ };
+
+ auto low_precision = [concentration](CoarseReal x)
+ {
+ boost::math::von_mises_distribution dist(0, concentration);
+ return boost::math::cdf(dist, x);
+ };
+
+ using hp_func = decltype (high_precision);
+
+ auto pi = boost::math::constants::pi();
+ boost::math::tools::ulps_plot plot(
+ high_precision, -pi, +pi);
+ plot.add_fn(low_precision);
+ std::string filename = "von_mises_ulps_cdf_"
+ + std::to_string(static_cast(concentration))
+ + typeid(CoarseReal).name()
+ + ".svg";
+ plot.write(filename);
+}
+
+void generate_ulps_plots(double concentration)
+{
+ float fconc = static_cast(concentration);
+ generate_ulps_plot_pdf(fconc);
+ generate_ulps_plot_cdf(fconc);
+
+ generate_ulps_plot_pdf(concentration);
+ generate_ulps_plot_cdf(concentration);
+}
+
+int main()
+{
+ generate_ulps_plots(4);
+ generate_ulps_plots(64);
+
+ return 0;
+}
diff --git a/test/Jamfile.v2 b/test/Jamfile.v2
index a98f8dc310..24350f0c3e 100644
--- a/test/Jamfile.v2
+++ b/test/Jamfile.v2
@@ -793,6 +793,9 @@ test-suite distribution_tests :
[ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : "-Bstatic -lquadmath -Bdynamic" ] ]
[ run test_triangular.cpp pch ../../test/build//boost_unit_test_framework ]
[ run test_uniform.cpp pch ../../test/build//boost_unit_test_framework ]
+ [ run test_von_mises.cpp pch ../../test/build//boost_unit_test_framework : : :
+ [ check-target-builds ../config//has_float128 "GCC libquadmath and __float128 support" : "-Bstatic -lquadmath -Bdynamic" ]
+ ]
[ run test_weibull.cpp ../../test/build//boost_unit_test_framework ]
[ run compile_test/dist_bernoulli_incl_test.cpp compile_test_main ]
diff --git a/test/test_von_mises.cpp b/test/test_von_mises.cpp
new file mode 100644
index 0000000000..3a89d17650
--- /dev/null
+++ b/test/test_von_mises.cpp
@@ -0,0 +1,397 @@
+// Copyright Philipp C. J. Münster 2020.
+// Copyright Paul A. Bristow 2010.
+// Copyright John Maddock 2007.
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// test_von_mises.cpp
+
+// https://en.wikipedia.org/wiki/Von_Mises_distribution
+// From MathWorld--A Wolfram Web Resource.
+// http://mathworld.wolfram.com/VonMisesDistribution.html
+
+#ifdef _MSC_VER
+# pragma warning (disable: 4127) // conditional expression is constant
+// caused by using if(std::numeric_limits::has_infinity)
+// and if (std::numeric_limits::has_quiet_NaN)
+#endif
+
+#include
+#include
+ using boost::multiprecision::float128;
+#include // for real_concept
+
+#define BOOST_TEST_MAIN
+#include // Boost.Test
+#include
+
+#include
+ using boost::math::von_mises_distribution;
+#include
+
+#include "math_unit_test.hpp"
+#include "pch.hpp"
+#include "test_out_of_range.hpp"
+
+#include
+#include
+ using std::cout;
+ using std::endl;
+ using std::setprecision;
+#include
+ using std::numeric_limits;
+
+template
+void check_von_mises(RealType mean, RealType conc, RealType x, RealType p, RealType q, RealType tol)
+{
+ BOOST_CHECK_CLOSE(
+ ::boost::math::cdf(
+ von_mises_distribution(mean, conc), // distribution.
+ x), // random variable.
+ p, // probability.
+ tol); // %tolerance.
+ BOOST_CHECK_CLOSE(
+ ::boost::math::cdf(
+ complement(
+ von_mises_distribution(mean, conc), // distribution.
+ x)), // random variable.
+ q, // probability complement.
+ tol); // %tolerance.
+ BOOST_CHECK_CLOSE(
+ ::boost::math::quantile(
+ von_mises_distribution(mean, conc), // distribution.
+ p), // probability.
+ x, // random variable.
+ tol); // %tolerance.
+ BOOST_CHECK_CLOSE(
+ ::boost::math::quantile(
+ complement(
+ von_mises_distribution(mean, conc), // distribution.
+ q)), // probability complement.
+ x, // random variable.
+ tol); // %tolerance.
+}
+
+template
+void test_spots(RealType)
+{
+ // Basic sanity checks
+ // Check some bad parameters to the distribution,
+#ifndef BOOST_NO_EXCEPTIONS
+ BOOST_MATH_CHECK_THROW(boost::math::von_mises_distribution nbad1(0, -1), std::domain_error); // negative conc
+#else
+ BOOST_MATH_CHECK_THROW(boost::math::von_mises_distribution(0, -1), std::domain_error); // negative conc
+#endif
+
+ // Tests on extreme values of random variate x, if has std::numeric_limits infinity etc.
+ von_mises_distribution N01;
+ if(std::numeric_limits::has_infinity)
+ {
+ BOOST_MATH_CHECK_THROW(pdf(N01, +std::numeric_limits::infinity()), std::domain_error); // x = + infinity, pdf = 0
+ BOOST_MATH_CHECK_THROW(pdf(N01, -std::numeric_limits::infinity()), std::domain_error); // x = - infinity, pdf = 0
+ BOOST_MATH_CHECK_THROW(cdf(N01, +std::numeric_limits::infinity()), std::domain_error); // x = + infinity, cdf = 1
+ BOOST_MATH_CHECK_THROW(cdf(N01, -std::numeric_limits::infinity()), std::domain_error); // x = - infinity, cdf = 0
+ BOOST_MATH_CHECK_THROW(cdf(complement(N01, +std::numeric_limits::infinity())), std::domain_error); // x = + infinity, c cdf = 0
+ BOOST_MATH_CHECK_THROW(cdf(complement(N01, -std::numeric_limits::infinity())), std::domain_error); // x = - infinity, c cdf = 1
+
+#ifndef BOOST_NO_EXCEPTIONS
+ BOOST_MATH_CHECK_THROW(boost::math::von_mises_distribution nbad1(std::numeric_limits::infinity(), static_cast(1)), std::domain_error); // +infinite mean
+ BOOST_MATH_CHECK_THROW(boost::math::von_mises_distribution nbad1(-std::numeric_limits::infinity(), static_cast(1)), std::domain_error); // -infinite mean
+ BOOST_MATH_CHECK_THROW(boost::math::von_mises_distribution nbad1(static_cast(0), std::numeric_limits::infinity()), std::domain_error); // infinite conc
+#else
+ BOOST_MATH_CHECK_THROW(boost::math::von_mises_distribution(std::numeric_limits::infinity(), static_cast(1)), std::domain_error); // +infinite mean
+ BOOST_MATH_CHECK_THROW(boost::math::von_mises_distribution(-std::numeric_limits::infinity(), static_cast(1)), std::domain_error); // -infinite mean
+ BOOST_MATH_CHECK_THROW(boost::math::von_mises_distribution(static_cast(0), std::numeric_limits::infinity()), std::domain_error); // infinite conc
+#endif
+ }
+
+ if (std::numeric_limits::has_quiet_NaN)
+ {
+ // No longer allow x to be NaN, then these tests should throw.
+ BOOST_MATH_CHECK_THROW(pdf(N01, +std::numeric_limits::quiet_NaN()), std::domain_error); // x = NaN
+ BOOST_MATH_CHECK_THROW(cdf(N01, +std::numeric_limits::quiet_NaN()), std::domain_error); // x = NaN
+ BOOST_MATH_CHECK_THROW(cdf(complement(N01, +std::numeric_limits::quiet_NaN())), std::domain_error); // x = + infinity
+ BOOST_MATH_CHECK_THROW(quantile(N01, +std::numeric_limits::quiet_NaN()), std::domain_error); // p = + infinity
+ BOOST_MATH_CHECK_THROW(quantile(complement(N01, +std::numeric_limits::quiet_NaN())), std::domain_error); // p = + infinity
+ }
+
+ //
+ // Tests for PDF: we know that the peak value is at e/(2*pi*I0(1)
+ //
+ RealType tolerance = boost::math::tools::epsilon() * 5 * 100; // 5 eps as a percentage
+ BOOST_CHECK_CLOSE(
+ pdf(von_mises_distribution(), static_cast(0)),
+ static_cast(0.34171048862346315949457814754706159394027L), // e/(2*pi*I0(1))
+ tolerance);
+ BOOST_CHECK_CLOSE(
+ pdf(von_mises_distribution(3, 0), static_cast(3)),
+ static_cast(0.15915494309189533576888376337251L), // 1/(2*pi)
+ tolerance);
+ BOOST_CHECK_CLOSE(
+ pdf(von_mises_distribution(3), static_cast(3)),
+ static_cast(0.34171048862346315949457814754706159394027L), // e/(2*pi*I0(1))
+ tolerance);
+ BOOST_CHECK_CLOSE(
+ pdf(von_mises_distribution(3, 5), static_cast(3)),
+ static_cast(0.86713652854235200257351846969777045343907L), // e^5/(2*pi*I0(5))
+ tolerance);
+ BOOST_CHECK_CLOSE(
+ pdf(von_mises_distribution(3, 25), static_cast(3)),
+ static_cast(1.98455543847726689510475504795539869409664L), // e^25/(2*pi*I0(25))
+ tolerance);
+ // edge case for single point precision
+ BOOST_CHECK_CLOSE(
+ pdf(von_mises_distribution(3, 86), static_cast(3)),
+ static_cast(3.69423343123704539549725123346713237943413L),
+ tolerance);
+ BOOST_CHECK_CLOSE(
+ pdf(von_mises_distribution(3, 87), static_cast(3)),
+ static_cast(3.71571226458759536792289974309199255119626L),
+ tolerance);
+ // edge case for double point precision
+ BOOST_CHECK_CLOSE(
+ pdf(von_mises_distribution(3, 708), static_cast(3)),
+ static_cast(10.6132883625399035032551439553585260831760L),
+ tolerance);
+ BOOST_CHECK_CLOSE(
+ pdf(von_mises_distribution(3, 709), static_cast(3)),
+ static_cast(10.6207836264247647802343430545802569228891L),
+ tolerance);
+
+ tolerance = 2e-3f; // 2e-5 (as %)
+
+ cout << "Tolerance for type " << typeid(RealType).name()
+ << " is " << tolerance << " %" << endl;
+
+ // test CDF for mean and interval edges
+ check_von_mises(
+ static_cast(0),
+ static_cast(1),
+ static_cast(0),
+ static_cast(0.5L),
+ static_cast(0.5L),
+ tolerance);
+
+ check_von_mises(
+ static_cast(0),
+ static_cast