From 18a871853a513ac236912b7c6d92d0c1fb4d981c Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Sun, 14 Aug 2022 20:39:38 -0700 Subject: [PATCH 1/5] Von Mises Distribution --- .../detail/common_error_handling.hpp | 22 + include/boost/math/distributions/fwd.hpp | 6 +- .../boost/math/distributions/von_mises.hpp | 907 ++++++++++++++++++ reporting/accuracy/von_mises_ulps.cpp | 88 ++ test/Jamfile.v2 | 1 + test/test_von_mises.cpp | 393 ++++++++ 6 files changed, 1416 insertions(+), 1 deletion(-) create mode 100644 include/boost/math/distributions/von_mises.hpp create mode 100644 reporting/accuracy/von_mises_ulps.cpp create mode 100644 test/test_von_mises.cpp diff --git a/include/boost/math/distributions/detail/common_error_handling.hpp b/include/boost/math/distributions/detail/common_error_handling.hpp index b91d2036f9..412209066f 100644 --- a/include/boost/math/distributions/detail/common_error_handling.hpp +++ b/include/boost/math/distributions/detail/common_error_handling.hpp @@ -11,6 +11,8 @@ #include #include +#include + // using boost::math::isfinite; // using boost::math::isnan; @@ -96,6 +98,26 @@ 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 a3c1a41df5..66dc54690d 100644 --- a/include/boost/math/distributions/fwd.hpp +++ b/include/boost/math/distributions/fwd.hpp @@ -117,6 +117,9 @@ class uniform_distribution; template class weibull_distribution; +template +class von_mises_distribution; + }} // namespaces #define BOOST_MATH_DECLARE_DISTRIBUTIONS(Type, Policy)\ @@ -152,6 +155,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::weibull_distribution weibull; + typedef boost::math::weibull_distribution weibull; \ + typedef boost::math::von_mises_distribution von_mises; #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..15334c51b0 --- /dev/null +++ b/include/boost/math/distributions/von_mises.hpp @@ -0,0 +1,907 @@ +// Copyright John Maddock 2006, 2007. +// Copyright Paul A. Bristow 2006, 2007. +// Copyright Philipp C. J. Muenster, 2020. +// Copyright Matt Borland, 2022. + +// 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 // for besseli0 and besseli1 +#include // for erf +#include +#include +#include + +#include +#include +#include +#include + +namespace boost { namespace math { + +template > +class von_mises_distribution +{ +public: + using value_type = RealType; + using policy_type = Policy; + + 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()); + } + + inline RealType mean() const + { // alias for location. + return m_mean; + } + + inline RealType concentration() const + { // alias for scale. + return m_concentration; + } + + // Synonyms, provided to allow generic use of find_location and find_scale. + inline RealType location() const + { // location. + return m_mean; + } + inline 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 + +using von_mises = von_mises_distribution; + +#ifdef __cpp_deduction_guides +template +von_mises_distribution(RealType)->von_mises_distribution>; +template +von_mises_distribution(RealType,RealType)->von_mises_distribution>; +#endif + +#ifdef BOOST_MSVC +#pragma warning(push) +#pragma warning(disable:4127) +#endif + +template +inline std::pair range(const von_mises_distribution& /*dist*/) +{ // Range of permissible values for random variable x. + BOOST_IF_CONSTEXPR (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 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. + constexpr RealType pi = boost::math::constants::pi(); + return std::pair(dist.mean() - pi, dist.mean() + 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 std::integral_constant&) +{ + const RealType mean = dist.mean(); + const 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 constexpr std::array 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 std::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 constexpr std::array 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 std::integral_constant&) +{ + const RealType mean = dist.mean(); + const 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 constexpr std::array 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 std::integral_constant&) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + const RealType mean = dist.mean(); + const RealType conc = dist.concentration(); + const 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 + + const RealType conc = dist.concentration(); + const RealType mean = dist.mean(); + + static const char* function = "boost::math::pdf(const von_mises_distribution<%1%>&, %1%)"; + + RealType result = 0; + if (!detail::check_positive_x(function, conc, &result, Policy())) + { + return result; + } + else if (!detail::check_angle(function, mean, &result, Policy())) + { + return result; + } + else if (!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; + //} + using tag_type = std::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 + >; + + 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); + + auto 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 (!detail::check_positive_x(function, conc, &result, Policy())) + { + return result; + } + else if (!detail::check_angle(function, mean, &result, Policy())) + { + return result; + } + if (!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 (!detail::check_positive_x(function, conc, &result, Policy())) + { + return result; + } + else if (!detail::check_angle(function, mean, &result, Policy())) + { + return result; + } + else if (!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(); + + using tag_type = std::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 + >; + + 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 (!detail::check_positive_x(function, conc, &result, Policy())) + { + return result; + } + if (!detail::check_angle(function, mean, &result, Policy())) + { + return result; + } + if (!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 (!detail::check_positive_x(function, conc, &result, Policy())) + { + return result; + } + else if (!detail::check_angle(function, mean, &result, Policy())) + { + return result; + } + + RealType q = c.param; + if (!detail::check_probability(function, q, &result, Policy())) + { + return result; + } + + if (q <= 0) + { + return +boost::math::constants::pi(); + } + else if (q >= 1) + { + return -boost::math::constants::pi(); + } + + using tag_type = std::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 + >; + + 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 std::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 constexpr std::array P_numer + { + +5.356107887570000000e-07f, + +1.994139882543095464e-01f, + +7.683426463016022940e-02f, + +4.007722563185265850e-02f, + +2.785578524715388070e-01f + }; + static constexpr std::array 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 constexpr std::array P_numer + { + 1.644239196640000000e-07f, + 1.994490498867993467e-01f, + 7.569820327857441460e-02f, + 5.573513685531774810e-02f, + 1.918908150536447035e-01f + }; + static constexpr std::array 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 std::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 constexpr std::array 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 constexpr std::array 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 constexpr std::array P_numer + { + +1.423e-15, + +1.994711401959018717e-01, + +7.480168411736836931e-02, + +7.012212565916144652e-02, + +1.037734243240472200e-01 + }; + static constexpr std::array 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 std::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 + { + // Polynomial for Bessel I0 / exp(conc) / sqrt(conc) + static constexpr std::array P_denom + { + BOOST_MATH_BIG_CONSTANT(RealType, 64, 3.9894228040143267793994605993438166526772e-01), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.9867785050179084742493257495245185241487e-02), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 2.8050629090725735167652437695397756897920e-02), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 2.9219405302839307466358297347675795965363e-02), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.4742214369972689474366968442268908028204e-02), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 9.0602984099194778006610058410222616383078e-02), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 2.2839502241666629677015839125593079416327e-01), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 6.8926354981801627920292655818232972385750e-01), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 2.4231921590621824187100989532173995000655e+00), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 9.7264260959693775207585700654645245723497e+00), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.3890136225398811195878046856373030127018e+01), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 2.1999720924619285464910452647408431234369e+02), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.2076909538525038580501368530598517194748e+03), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 7.5684635141332367730007149159063086133399e+03), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 3.5178192543258299267923025833141286569141e+04), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 6.2966297919851965784482163987240461837728e+05) + }; + + // Polynomial for (Bessel I0 - Bessel I1) / exp(conc) / sqrt(x) + static constexpr std::array P_numer + { + BOOST_MATH_BIG_CONSTANT(RealType, 64, -3.368165e-34), + BOOST_MATH_BIG_CONSTANT(RealType, 64, +0.1994711402007163389699730299733589146073), + BOOST_MATH_BIG_CONSTANT(RealType, 64, +0.07480167757526862711373984952175456888896), + BOOST_MATH_BIG_CONSTANT(RealType, 64, +0.07012657272681433791923412617430580227103), + BOOST_MATH_BIG_CONSTANT(RealType, 64, +0.10226791855993757596915805134093796837369), + BOOST_MATH_BIG_CONSTANT(RealType, 64, +0.2013399646648772644282448264813045181469), + BOOST_MATH_BIG_CONSTANT(RealType, 64, +0.4983164125454637874163933874417571110056), + BOOST_MATH_BIG_CONSTANT(RealType, 64, +1.4845676457582822590839158984567308297366), + BOOST_MATH_BIG_CONSTANT(RealType, 64, +5.169476606935535670414552625141409318434), + BOOST_MATH_BIG_CONSTANT(RealType, 64, +20.59713743665130419014001937211862931161), + BOOST_MATH_BIG_CONSTANT(RealType, 64, +92.40031163861578044111910646492625263225), + BOOST_MATH_BIG_CONSTANT(RealType, 64, +460.9440321063085921197536056693132597443), + BOOST_MATH_BIG_CONSTANT(RealType, 64, +2.520575547528944544570101030955801999019e+03), + BOOST_MATH_BIG_CONSTANT(RealType, 64, +1.5734053646329490893326466550032992462076e+04), + BOOST_MATH_BIG_CONSTANT(RealType, 64, +7.319778356894459435808347175389511056370e+04), + BOOST_MATH_BIG_CONSTANT(RealType, 64, +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; + + return 0; + } +} + +// quad version of variance_impl +template +inline RealType variance_impl(const von_mises_distribution& dist, + const std::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 std::array 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 std::array 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) +{ + + using tag_type = std::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 + >; + + 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 \ No newline at end of file diff --git a/reporting/accuracy/von_mises_ulps.cpp b/reporting/accuracy/von_mises_ulps.cpp new file mode 100644 index 0000000000..a50953c17d --- /dev/null +++ b/reporting/accuracy/von_mises_ulps.cpp @@ -0,0 +1,88 @@ +// Copyright Philipp C. J. Muenster, 2020. +// Copyright Matt Borland, 2022. + +// 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) + +#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) +{ + auto 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; +} \ No newline at end of file diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 17aa8484a0..3dfc131388 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -845,6 +845,7 @@ test-suite distribution_tests : [ run test_triangular.cpp pch ../../test/build//boost_unit_test_framework ] [ run test_uniform.cpp pch ../../test/build//boost_unit_test_framework ] [ run test_weibull.cpp ../../test/build//boost_unit_test_framework ] + [ run test_von_mises.cpp ../../test/build//boost_unit_test_framework ] [ run compile_test/dist_bernoulli_incl_test.cpp compile_test_main : : : [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] [ run compile_test/dist_beta_incl_test.cpp compile_test_main : : : [ check-target-builds ../config//is_ci_sanitizer_run "Sanitizer CI run" : no ] ] diff --git a/test/test_von_mises.cpp b/test/test_von_mises.cpp new file mode 100644 index 0000000000..1002bd084c --- /dev/null +++ b/test/test_von_mises.cpp @@ -0,0 +1,393 @@ + +// Copyright Philipp C. J. Muenster 2020. +// Copyright Paul A. Bristow 2010. +// Copyright John Maddock 2007. +// Copyright Matt Borland 2022. + +// 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 // 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(1), + static_cast(-boost::math::constants::pi()), + static_cast(0.0L), + static_cast(1.0L), + tolerance); + + check_von_mises( + static_cast(0), + static_cast(1), + static_cast(boost::math::constants::pi()), + static_cast(1.0L), + static_cast(0.0L), + tolerance); + + // test CDF for low concentrations + check_von_mises( + static_cast(0), + static_cast(1), + static_cast(1), + static_cast(0.794355307434683479987678129735260058645499262629455722769L), + static_cast(0.205644692565316520012321870264739941354500737370544277230L), + tolerance); + + check_von_mises( + static_cast(0), + static_cast(1), + static_cast(-1), + static_cast(0.205644692565316520012321870264739941354500737370544277230L), + static_cast(0.794355307434683479987678129735260058645499262629455722769L), + tolerance); + + check_von_mises( + static_cast(0), + static_cast(5), + static_cast(1), + static_cast(0.98096204546814689054581384796251763480020360394758184271L), + static_cast(0.01903795453185310945418615203748236519979639605241815729L), + tolerance); + + check_von_mises( + static_cast(0), + static_cast(5), + static_cast(-1), + static_cast(0.01903795453185310945418615203748236519979639605241815729L), + static_cast(0.98096204546814689054581384796251763480020360394758184271L), + tolerance); + + // test CDF for high concentrations + //~ check_von_mises( + //~ static_cast(0), + //~ static_cast(25), + //~ static_cast(1), + //~ static_cast(0.999999062404464440452233504489299782776166264467210572765L), + //~ static_cast(9.37595535559547766495510700217223833735532789427234e-7L), + //~ tolerance); + + //~ check_von_mises( + //~ static_cast(0), + //~ static_cast(25), + //~ static_cast(-1), + //~ static_cast(9.37595535559547766495510700217223833735532789427234e-7L), + //~ static_cast(0.999999062404464440452233504489299782776166264467210572765L), + //~ tolerance); + + //~ check_von_mises( + //~ static_cast(0), + //~ static_cast(125), + //~ static_cast(1), + //~ static_cast(0.99999999999999999996645363431349332910951002333081490389L), + //~ static_cast(3.3546365686506670890489976669185096109e-20L), + //~ tolerance); + + //~ check_von_mises( + //~ static_cast(0), + //~ static_cast(125), + //~ static_cast(-1), + //~ static_cast(3.3546365686506670890489976669185096109e-20L), + //~ static_cast(0.99999999999999999996645363431349332910951002333081490389L), + //~ tolerance); + + + RealType tol2 = boost::math::tools::epsilon() * 500; + von_mises_distribution dist(2, 3); + RealType x = static_cast(0.125); + + BOOST_MATH_STD_USING // ADL of std math lib names + + // mean: + BOOST_CHECK_CLOSE( + mean(dist) + , static_cast(2), tol2); + // variance: + BOOST_CHECK_CLOSE( + variance(von_mises_distribution(2, 0)) + , static_cast(1), tol2); + BOOST_CHECK_CLOSE( + variance(von_mises_distribution(2, 1)) + , static_cast(0.553610034103465492952318204807357330223746852599612177138L), tol2); + BOOST_CHECK_CLOSE( + variance(von_mises_distribution(2, 5)) + , static_cast(0.106616862955914778412994992775050827245562823701700738786L), tol2); + //~ BOOST_CHECK_CLOSE( + //~ variance(von_mises_distribution(2, 25)) + //~ , static_cast(0.020208546509484068780303296944566388571293465011951716357L), tol2); + //~ BOOST_CHECK_CLOSE( + //~ variance(von_mises_distribution(2, 125)) + //~ , static_cast(0.004008064813593637057963834031301840442156903531539609019L), tol2); + + // std deviation: + BOOST_CHECK_CLOSE( + standard_deviation(dist) + , static_cast(0.649213658343262740252572725779410261163523458085389547123L), tol2); + // hazard: + BOOST_CHECK_CLOSE( + hazard(dist, x) + , pdf(dist, x) / cdf(complement(dist, x)), tol2); + // cumulative hazard: + BOOST_CHECK_CLOSE( + chf(dist, x) + , -log(cdf(complement(dist, x))), tol2); + // coefficient_of_variation: + BOOST_CHECK_CLOSE( + coefficient_of_variation(dist) + , standard_deviation(dist) / mean(dist), tol2); + // mode: + BOOST_CHECK_CLOSE( + mode(dist) + , static_cast(2), tol2); + + BOOST_CHECK_CLOSE( + median(dist) + , static_cast(2), tol2); + + // skewness: + BOOST_CHECK_CLOSE( + skewness(dist) + , static_cast(0), tol2); + + BOOST_CHECK_CLOSE( + entropy(dist) + , 0.993228806353252817873771736349142645476889275244864748502L, tol2); + + von_mises_distribution norm01(0, 1); // Test default (0, 1) + BOOST_CHECK_CLOSE( + mean(norm01), + static_cast(0), 0); // Mean == zero + + von_mises_distribution defsd_norm01(0); // Test default (0, sd = 1) + BOOST_CHECK_CLOSE( + mean(defsd_norm01), + static_cast(0), 0); // Mean == zero + + von_mises_distribution def_norm01; // Test default (0, sd = 1) + BOOST_CHECK_CLOSE( + mean(def_norm01), + static_cast(0), 0); // Mean == zero + + // Error tests: + check_out_of_range >(0, 1); // (All) valid constructor parameter values. + + BOOST_MATH_CHECK_THROW(quantile(von_mises_distribution(0, 1), -1), std::domain_error); + BOOST_MATH_CHECK_THROW(quantile(von_mises_distribution(0, 1), 2), std::domain_error); +} // template void test_spots(RealType) + +template +void test_symmetry(RealType) +{ + RealType const pi = boost::math::constants::pi(); + RealType delta = 1.0 / (1 << 4); + for (RealType mean = 0; mean < pi; mean += delta) { + for (RealType conc = 0; conc < 100; conc = (conc + 1) * 1.5 - 1) { + von_mises_distribution dist(mean, conc); + for (RealType x = 0; x < pi; x += delta) { + CHECK_ULP_CLOSE(pdf(dist, mean + x), + pdf(dist, mean - x), 2); + CHECK_ULP_CLOSE(cdf(dist, mean + x) - static_cast(0.5), + static_cast(0.5) - cdf(dist, mean - x), 32); + } + } + } +} + +BOOST_AUTO_TEST_CASE( test_main ) +{ + // Check that we can generate von_mises distribution using the two convenience methods: + boost::math::von_mises myf1(1., 2); // Using typedef + von_mises_distribution<> myf2(1., 2); // Using default RealType double. + boost::math::von_mises myn01; // Use default values. + // Note NOT myn01() as the compiler will interpret as a function! + + // Check the synonyms, provided to allow generic use of find_location and find_scale. + BOOST_CHECK_EQUAL(myn01.mean(), myn01.location()); + BOOST_CHECK_EQUAL(myn01.concentration(), myn01.scale()); + + // Basic sanity-check spot values. + // (Parameter value, arbitrarily zero, only communicates the floating point type). + test_spots(0.0F); // Test float. OK at decdigits = 0 tolerance = 0.0001 % + test_spots(0.0); // Test double. OK at decdigits 7, tolerance = 1e07 % + test_spots(0.0L); // Test long double. + + // Check symmetry of PDF and CDF + test_symmetry(0.0F); + test_symmetry(0.0); + test_symmetry(0.0L); +} // BOOST_AUTO_TEST_CASE( test_main ) + +/* +./test_von_mises.exe +Output: +Running 1 test case... +Tolerance for type f is 0.002 % +Tolerance for type d is 0.002 % +Tolerance for type e is 0.002 % +*** No errors detected */ From 9f3626a153f873e720451ef25326fd65514637ac Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Tue, 16 Aug 2022 08:52:01 -0700 Subject: [PATCH 2/5] Fix formatting issues [ci skip] --- include/boost/math/distributions/von_mises.hpp | 4 ++-- reporting/accuracy/von_mises_ulps.cpp | 2 +- test/test_von_mises.cpp | 12 ++++++------ 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/include/boost/math/distributions/von_mises.hpp b/include/boost/math/distributions/von_mises.hpp index 15334c51b0..c07fb59114 100644 --- a/include/boost/math/distributions/von_mises.hpp +++ b/include/boost/math/distributions/von_mises.hpp @@ -107,7 +107,7 @@ template inline 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. constexpr RealType pi = boost::math::constants::pi(); - return std::pair(dist.mean() - pi, dist.mean() + pi); // [µ-π, µ+π) + return std::pair(dist.mean() - pi, dist.mean() + pi); // u-pi, u+pi } #ifdef BOOST_MSVC @@ -904,4 +904,4 @@ inline RealType entropy(const von_mises_distribution & dist) // keep compilers that support two-phase lookup happy. #include -#endif // BOOST_MATH_DISTRIBUTIONS_VON_MISES_HPP \ No newline at end of file +#endif // BOOST_MATH_DISTRIBUTIONS_VON_MISES_HPP diff --git a/reporting/accuracy/von_mises_ulps.cpp b/reporting/accuracy/von_mises_ulps.cpp index a50953c17d..9a4e79c633 100644 --- a/reporting/accuracy/von_mises_ulps.cpp +++ b/reporting/accuracy/von_mises_ulps.cpp @@ -85,4 +85,4 @@ int main() generate_ulps_plots(64); return 0; -} \ No newline at end of file +} diff --git a/test/test_von_mises.cpp b/test/test_von_mises.cpp index 1002bd084c..f3f28fd433 100644 --- a/test/test_von_mises.cpp +++ b/test/test_von_mises.cpp @@ -142,8 +142,8 @@ void test_spots(RealType) 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( + // edge case for single point precision + BOOST_CHECK_CLOSE( pdf(von_mises_distribution(3, 86), static_cast(3)), static_cast(3.69423343123704539549725123346713237943413L), tolerance); @@ -151,12 +151,12 @@ void test_spots(RealType) pdf(von_mises_distribution(3, 87), static_cast(3)), static_cast(3.71571226458759536792289974309199255119626L), tolerance); - // edge case for double point precision + // 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( + BOOST_CHECK_CLOSE( pdf(von_mises_distribution(3, 709), static_cast(3)), static_cast(10.6207836264247647802343430545802569228891L), tolerance); @@ -344,8 +344,8 @@ void test_spots(RealType) template void test_symmetry(RealType) { - RealType const pi = boost::math::constants::pi(); - RealType delta = 1.0 / (1 << 4); + RealType const pi = boost::math::constants::pi(); + RealType delta = 1.0 / (1 << 4); for (RealType mean = 0; mean < pi; mean += delta) { for (RealType conc = 0; conc < 100; conc = (conc + 1) * 1.5 - 1) { von_mises_distribution dist(mean, conc); From 51171c1731f43c621607405595c4bea18e91c542 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Wed, 17 Aug 2022 20:26:33 -0700 Subject: [PATCH 3/5] Begin re-writing cdf_impl [ci skip] --- .../boost/math/distributions/von_mises.hpp | 100 +++++++++++++++++- 1 file changed, 98 insertions(+), 2 deletions(-) diff --git a/include/boost/math/distributions/von_mises.hpp b/include/boost/math/distributions/von_mises.hpp index c07fb59114..95b0530a73 100644 --- a/include/boost/math/distributions/von_mises.hpp +++ b/include/boost/math/distributions/von_mises.hpp @@ -106,7 +106,7 @@ inline std::pair range(const von_mises_distribution inline 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. - constexpr RealType pi = boost::math::constants::pi(); + const RealType pi = boost::math::constants::pi(); return std::pair(dist.mean() - pi, dist.mean() + pi); // u-pi, u+pi } @@ -292,6 +292,101 @@ inline RealType pdf(const von_mises_distribution& dist, const namespace detail { + +// 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 +template +RealType cdf_impl(const von_mises_distribution& dist, const RealType& x) +{ + BOOST_MATH_STD_USING + + const RealType conc = dist.concentration(); + const RealType mean = dist.mean(); + const RealType pi = boost::math::constants::pi(); + const RealType two_pi = boost::math::constants::two_pi(); + + RealType u = x - mean; + + if (u <= -pi) + { + return 0; + } + else if (u >= pi) + { + return 1; + } + + constexpr auto digits = std::numeric_limits::max_digits10 - 1; + constexpr auto digits_cubed = digits*digits*digits; + constexpr auto digits_sq = digits*digits; + + // Interplation of critical kappa values described in the above paper + // https://www.wolframalpha.com/input?i=interpolate%5B%286%2C+6.5%29%2C+%287%2C+8.0%29%2C+%288%2C+10.5%29%2C+%289%2C+15%29%2C+%2810%2C+21%29%2C+%2811%2C+32%29%2C+%2812%2C+50%29%5D + constexpr auto critical_kappa = -0.0152778*digits_cubed*digits_cubed + 0.825*digits_cubed*digits_sq - + 18.3194*digits_sq*digits_sq + 214.25*digits_cubed - 1391.67*digits_sq + + 4760.43*digits - 6694.5; + + // We can use the pre-calculated values from the paper if the precision is less than twelve digits + // otherwise we are going to have to calulate ourselves. + RealType result; + BOOST_IF_CONSTEXPR (digits <= 12) + { + if (conc > critical_kappa) + { + 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); + + auto 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(); + } + } + else + { + + } + + if (result < 0) + { + result = 0; + } + else if (result > 1) + { + result = 1; + } + + return result; +} + +/* template inline RealType cdf_impl(const von_mises_distribution& dist, const RealType& x) { @@ -316,7 +411,7 @@ inline RealType cdf_impl(const von_mises_distribution& dist, c // DOI: 10.1145/355744.355753 RealType result = 0; - int digits = std::numeric_limits::max_digits10 - 1; + constexpr int digits = std::numeric_limits::max_digits10 - 1; RealType ck = ((0.1611*digits - 2.8778)*digits + 18.45)*digits - 35.4; if (conc > ck) { @@ -355,6 +450,7 @@ inline RealType cdf_impl(const von_mises_distribution& dist, c return result <= 0 ? 0 : (1 <= result ? 1 : result); } +*/ } // namespace detail template From 2149eeba25f02c2f2ba4252ab86c2661e43841a2 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Tue, 23 Aug 2022 19:56:39 -0700 Subject: [PATCH 4/5] Remove double precision constants from CDF impl --- .../boost/math/distributions/von_mises.hpp | 194 +++++------------- 1 file changed, 52 insertions(+), 142 deletions(-) diff --git a/include/boost/math/distributions/von_mises.hpp b/include/boost/math/distributions/von_mises.hpp index 95b0530a73..6b3fe8e1e7 100644 --- a/include/boost/math/distributions/von_mises.hpp +++ b/include/boost/math/distributions/von_mises.hpp @@ -33,8 +33,8 @@ template > class von_mises_distribution { public: - using value_type = RealType; - using policy_type = Policy; + using value_type = RealType; + using policy_type = Policy; von_mises_distribution(RealType l_mean = 0, RealType concentration = 1) : m_mean {l_mean}, m_concentration {concentration} @@ -275,16 +275,11 @@ inline RealType pdf(const von_mises_distribution& dist, const // return 0; //} using tag_type = std::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 + ((std::numeric_limits::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 >; return detail::pdf_impl(dist, x, tag_type()); @@ -292,165 +287,86 @@ inline RealType pdf(const von_mises_distribution& dist, const namespace detail { - // 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 template -RealType cdf_impl(const von_mises_distribution& dist, const RealType& x) -{ - BOOST_MATH_STD_USING - - const RealType conc = dist.concentration(); - const RealType mean = dist.mean(); - const RealType pi = boost::math::constants::pi(); - const RealType two_pi = boost::math::constants::two_pi(); - - RealType u = x - mean; - - if (u <= -pi) - { - return 0; - } - else if (u >= pi) - { - return 1; - } - - constexpr auto digits = std::numeric_limits::max_digits10 - 1; - constexpr auto digits_cubed = digits*digits*digits; - constexpr auto digits_sq = digits*digits; - - // Interplation of critical kappa values described in the above paper - // https://www.wolframalpha.com/input?i=interpolate%5B%286%2C+6.5%29%2C+%287%2C+8.0%29%2C+%288%2C+10.5%29%2C+%289%2C+15%29%2C+%2810%2C+21%29%2C+%2811%2C+32%29%2C+%2812%2C+50%29%5D - constexpr auto critical_kappa = -0.0152778*digits_cubed*digits_cubed + 0.825*digits_cubed*digits_sq - - 18.3194*digits_sq*digits_sq + 214.25*digits_cubed - 1391.67*digits_sq - + 4760.43*digits - 6694.5; - - // We can use the pre-calculated values from the paper if the precision is less than twelve digits - // otherwise we are going to have to calulate ourselves. - RealType result; - BOOST_IF_CONSTEXPR (digits <= 12) - { - if (conc > critical_kappa) - { - 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); - - auto 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(); - } - } - else - { - - } - - if (result < 0) - { - result = 0; - } - else if (result > 1) - { - result = 1; - } - - return result; -} - -/* -template inline RealType cdf_impl(const von_mises_distribution& dist, const RealType& x) { BOOST_MATH_STD_USING // for ADL of std functions + constexpr RealType pi = boost::math::constants::pi(); + RealType conc = dist.concentration(); RealType mean = dist.mean(); RealType u = x - mean; - if (u <= -boost::math::constants::pi()) + if (u <= -pi) { return 0; } - if (u >= +boost::math::constants::pi()) + if (u >= 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; + // Extrapolation of critical kappa value constexpr int digits = std::numeric_limits::max_digits10 - 1; - RealType ck = ((0.1611*digits - 2.8778)*digits + 18.45)*digits - 35.4; + constexpr RealType ck = ((0.1611 * digits - 2.8778) * digits + 18.45) * digits - 35.4; + if (conc > ck) { - RealType c = 24.0 * conc; + RealType c = 24 * 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 r = sqrt((54 / (347 / v + 26 - c) - 6 + c) / 12); + RealType z = sin(u / 2) * 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; + RealType y = (c - s - s - 16) / 3; + y = ((s + static_cast(7)/4) * s + static_cast(167)/2) / v - y; + result = (1 + boost::math::erf(z - s / (y * y) * z)) / 2; } - else + 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); + if(conc > 0) + { + // Extrapolation of the tables given in the paper used to estimate the number of iterations needed + constexpr RealType a1 = 0.357143 * digits * digits - 3.13286 * digits + 13.9286; + constexpr RealType a2 = 1.53571 - 0.0857143 * digits; + constexpr RealType a3 = 1.01389 * digits * digits * digits - 23.1488 * digits * digits + 177.349 * digits - 448.19; + constexpr RealType a4 = 0.178571 * digits * digits - 2.55 * digits + 9.87143; auto 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) { + for (auto 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(); + result = (x - mean + pi) / 2; + result = (result + v) / pi; + } + + if (result < 0) + { + result = 0; + } + else if (result > 1) + { + result = 1; } - return result <= 0 ? 0 : (1 <= result ? 1 : result); + return result; } -*/ + } // namespace detail template @@ -509,13 +425,11 @@ inline RealType quantile(const von_mises_distribution& dist, c return +boost::math::constants::pi(); using tag_type = std::integral_constant::digits == 0) - || (std::numeric_limits::radix != 2)) ? 0 : + ((std::numeric_limits::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 + std::numeric_limits::digits <= 113 ? 113 : -1 >; struct step_func @@ -523,8 +437,8 @@ inline RealType quantile(const von_mises_distribution& dist, c 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) + return std::make_pair(detail::cdf_impl(dist, x) - p, // f(x) + detail::pdf_impl(dist, x, tag_type())); // f'(x) } }; @@ -600,13 +514,11 @@ inline RealType quantile(const complemented2_type::digits == 0) - || (std::numeric_limits::radix != 2)) ? 0 : + ((std::numeric_limits::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 + std::numeric_limits::digits <= 113 ? 113 : -1 >; struct step_func @@ -964,13 +876,11 @@ inline RealType variance(const von_mises_distribution& dist) { using tag_type = std::integral_constant::digits == 0) - || (std::numeric_limits::radix != 2)) ? 0 : + ((std::numeric_limits::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 + std::numeric_limits::digits <= 113 ? 113 : -1 >; return detail::variance_impl(dist, tag_type()); From aad56c9c3edd4240f27bddfa7a4b4606962fa714 Mon Sep 17 00:00:00 2001 From: Matt Borland Date: Fri, 26 Aug 2022 15:12:47 -0700 Subject: [PATCH 5/5] Add documentation [ci skip] --- doc/distributions/dist_reference.qbk | 1 + doc/distributions/von_mises.qbk | 102 +++ doc/images/von_mises_ulps_cdf_4d.svg | 1036 ++++++++++++++++++++++++ doc/images/von_mises_ulps_cdf_64d.svg | 1038 +++++++++++++++++++++++++ doc/images/von_mises_ulps_pdf_4d.svg | 1038 +++++++++++++++++++++++++ doc/images/von_mises_ulps_pdf_64d.svg | 1038 +++++++++++++++++++++++++ 6 files changed, 4253 insertions(+) create mode 100644 doc/distributions/von_mises.qbk create mode 100644 doc/images/von_mises_ulps_cdf_4d.svg create mode 100644 doc/images/von_mises_ulps_cdf_64d.svg create mode 100644 doc/images/von_mises_ulps_pdf_4d.svg create mode 100644 doc/images/von_mises_ulps_pdf_64d.svg diff --git a/doc/distributions/dist_reference.qbk b/doc/distributions/dist_reference.qbk index c225d1953e..93d520db4f 100644 --- a/doc/distributions/dist_reference.qbk +++ b/doc/distributions/dist_reference.qbk @@ -38,6 +38,7 @@ [include students_t.qbk] [include triangular.qbk] [include uniform.qbk] +[include von_mises.qbk] [include weibull.qbk] [endsect] [/section:dists Distributions] diff --git a/doc/distributions/von_mises.qbk b/doc/distributions/von_mises.qbk new file mode 100644 index 0000000000..7789ef449d --- /dev/null +++ b/doc/distributions/von_mises.qbk @@ -0,0 +1,102 @@ +[/ + Copyright 2022 Matt Borland. + Distributed under 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). +] + +[section:von_mises_dist Von Mises Distribution] + + +``#include `` + + namespace boost{ namespace math{ + template + class von_mises_distribution; + + using von_mises = von_mises_distribution<>; + + template + class von_mises_distribution + { + public: + using value_type = RealType; + using policy_type = Policy; + + von_mises_distribution(RealType mean = 0, RealType concentration = 1); // Constructor. + : m_mean {maen}, m_conc {concentration} // Default is standard von_mises distribution. + // Accessor functions. + RealType mean() const; + RealType concentration() const; + RealType location() const; + RealType scale() const; + }; // class von_mises_distribution + + }} // namespaces + +The [@https://en.wikipedia.org/wiki/Von_Mises_distribution von Mises distribution], also known as a the circular normal distribution, +is a continuous probability distribution on a circle. + +[h4 Member Functions] + + von_mises_distribution(RealType mean = 0, RealType concentration = 1); ; + +Constructs a [@https://en.wikipedia.org/wiki/Von_Mises_distribution +von Mises Distribution] with mean /mean/, and concentration /concentration/. + +Requires that the /mean/ and /concentration/ parameters are both finite; +otherwise if infinity or NaN then calls __domain_error. + + RealType mean() const; + +Returns the /mean/ parameter of this distribution. + + RealType concentration() const; + +Returns the /concentration/ parameter of this distribution. + + RealType location() const; + +Returns the /mean/ parameter as the location is synonymous but used in some definitions. + + RealType scale() const; + +Returns the /concentration/ parameter as the location is synonymous but used in some definitions. + +[h4 Non-member Accessors] + +All the [link math_toolkit.dist_ref.nmp usual non-member accessor functions] +that are generic to all distributions are supported: __usual_accessors. + +The domain of the random variable is O <= x <= 2[pi]. + +[h4 Accuracy] + +The CDF of this distribution is not analytic so an approximation is made. +For the accuracy on your specific platform please run reporting/accuracy/von_mises_ulps.cpp + +An example on Apple M1 Pro is below using double precsision calculations: + +*PDF with concentration of 4* + +[$../images/von_mises_ulps_pdf_4d.svg] + +*PDF with concentration of 64* + +[$../images/von_mises_ulps_pdf_64d.svg] + +*CDF with concentration of 4* + +[$../images/von_mises_ulps_cdf_4d.svg] + +*CDF with concentration of 64* + +[$../images/von_mises_ulps_cdf_64d.svg] + +[h4 References] +* [@https://en.wikipedia.org/wiki/Von_Mises_distribution Wikipedia von Mises distribution] +* [@https://mathworld.wolfram.com/vonMisesDistribution.html Weisstein, Eric W. "von Mises Distribution." From MathWorld--A Wolfram Web Resource.] +* [@https://dl.acm.org/doi/pdf/10.1145/355744.355753] + +[endsect] [/section:von_mises_dist Von Mises] diff --git a/doc/images/von_mises_ulps_cdf_4d.svg b/doc/images/von_mises_ulps_cdf_4d.svg new file mode 100644 index 0000000000..5a043aca19 --- /dev/null +++ b/doc/images/von_mises_ulps_cdf_4d.svg @@ -0,0 +1,1036 @@ + + + + + + + +-0.5 + +0.5 + +-2.513 + +-1.885 + +-1.257 + +-0.6283 + +0 + +0.6283 + +1.257 + +1.885 + +2.513 + +3.142 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc/images/von_mises_ulps_cdf_64d.svg b/doc/images/von_mises_ulps_cdf_64d.svg new file mode 100644 index 0000000000..843fcc2be2 --- /dev/null +++ b/doc/images/von_mises_ulps_cdf_64d.svg @@ -0,0 +1,1038 @@ + + + + + + + +-0.5 + +0.5 + +-2.513 + +-1.885 + +-1.257 + +-0.6283 + +0 + +0.6283 + +1.257 + +1.885 + +2.513 + +3.142 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc/images/von_mises_ulps_pdf_4d.svg b/doc/images/von_mises_ulps_pdf_4d.svg new file mode 100644 index 0000000000..095b4c4d63 --- /dev/null +++ b/doc/images/von_mises_ulps_pdf_4d.svg @@ -0,0 +1,1038 @@ + + + + + + + +-0.5 + +0.5 + +0.3142 + +0.6283 + +0.9425 + +1.257 + +1.571 + +1.885 + +2.199 + +2.513 + +2.827 + +3.142 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/doc/images/von_mises_ulps_pdf_64d.svg b/doc/images/von_mises_ulps_pdf_64d.svg new file mode 100644 index 0000000000..dc81280c2d --- /dev/null +++ b/doc/images/von_mises_ulps_pdf_64d.svg @@ -0,0 +1,1038 @@ + + + + + + + +-0.5 + +0.5 + +0.3142 + +0.6283 + +0.9425 + +1.257 + +1.571 + +1.885 + +2.199 + +2.513 + +2.827 + +3.142 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +