From e33f6b5360c0141f03b05644b40eedaa2b1a6b8b Mon Sep 17 00:00:00 2001 From: pcjmuenster Date: Sun, 26 Jan 2020 15:40:53 +0100 Subject: [PATCH 1/8] Initial commit for Von Mises Distribution Add class and test file; tests for PDF pass for single precision --- .../detail/common_error_handling.hpp | 20 + .../boost/math/distributions/von_mises.hpp | 440 ++++++++++++++++++ test/Jamfile.v2 | 1 + test/test_von_mises.cpp | 397 ++++++++++++++++ 4 files changed, 858 insertions(+) create mode 100644 include/boost/math/distributions/von_mises.hpp 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 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/von_mises.hpp b/include/boost/math/distributions/von_mises.hpp new file mode 100644 index 0000000000..f2f8598591 --- /dev/null +++ b/include/boost/math/distributions/von_mises.hpp @@ -0,0 +1,440 @@ +// Copyright John Maddock 2006, 2007. +// Copyright Paul A. Bristow 2006, 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) + +#ifndef BOOST_STATS_VON_MISES_HPP +#define BOOST_STATS_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 erf/erfc. +#include +#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(-constants::pi(), + +constants::pi()); // -pi to +pi +} + +#ifdef BOOST_MSVC +#pragma warning(pop) +#endif + +// float version of log_bessel_i0 +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 < 7.75) + { + RealType bessel_i0 = cyl_bessel_i(0, conc, Policy()); + return exp(conc * cos(x - mean)) / bessel_i0 / boost::math::constants::two_pi(); + } + else if(conc < 50) + { + // Max error in interpolated form: 5.195e-08 + // Max Error found at float precision = Poly: 8.502534e-08 + static const float P[] = { + 3.98942651588301770e-01f, + 4.98327234176892844e-02f, + 2.91866904423115499e-02f, + 1.35614940793742178e-02f, + 1.31409251787866793e-01f + }; + RealType exponent = conc * (cos(x - mean) - 1.f); + exponent -= std::log(boost::math::tools::evaluate_polynomial(P, RealType(1 / conc)) / sqrt(conc)); + return exp(exponent) / boost::math::constants::two_pi(); + } + else + { + // Max error in interpolated form: 1.782e-09 + // Max Error found at float precision = Poly: 6.473568e-08 + static const float P[] = { + 3.98942280401432677e-01f, + 4.98677850501790847e-02f, + 2.80506290907257351e-02f, + 2.92194053028393074e-02f, + 4.47422143699726895e-02f, + }; + RealType exponent = conc * (cos(x - mean) - 1.f); + exponent -= std::log(boost::math::tools::evaluate_polynomial(P, RealType(1 / conc)) / sqrt(conc)); + return exp(exponent) / boost::math::constants::two_pi(); + } +} + +// double version of log_bessel_i0 +template +inline RealType log_bessel_i0_impl(RealType x, const boost::integral_constant&) +{ + BOOST_MATH_STD_USING + if(x < 7.75) + { + static const double P[] = { + 1.00000000000000000e+00, + 2.49999999999999909e-01, + 2.77777777777782257e-02, + 1.73611111111023792e-03, + 6.94444444453352521e-05, + 1.92901234513219920e-06, + 3.93675991102510739e-08, + 6.15118672704439289e-10, + 7.59407002058973446e-12, + 7.59389793369836367e-14, + 6.27767773636292611e-16, + 4.34709704153272287e-18, + 2.63417742690109154e-20, + 1.13943037744822825e-22, + 9.07926920085624812e-25 + }; + RealType a = x * x / 4; + return std::log(a * boost::math::tools::evaluate_polynomial(P, a) + 1); + } + else if(x < 500) + { + static const double P[] = { + 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 + }; + return x + std::log(boost::math::tools::evaluate_polynomial(P, RealType(1 / x)) / sqrt(x)); + } + else + { + static const double P[] = { + 3.98942280401432905e-01, + 4.98677850491434560e-02, + 2.80506308916506102e-02, + 2.92179096853915176e-02, + 4.53371208762579442e-02 + }; + RealType ex = x / 2; + RealType result = ex + std::log(boost::math::tools::evaluate_polynomial(P, RealType(1 / x)) / sqrt(x)); + result += ex; + return result; + } +} + +template +inline RealType log_bessel_i0(RealType x) +{ + 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 log_bessel_i0_impl(x, tag_type()); +} + +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_x(function, x, &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 pdf_impl(dist, x, tag_type()); +} // pdf + +template +inline RealType cdf(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::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, &result, Policy())) + { + return result; + } + + RealType diff = (x - mean) / (conc * constants::root_two()); + result = boost::math::erfc(-diff, Policy()) / 2; + return result; +} // 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; + + // result = boost::math::erfc_inv(2 * p, Policy()); + // result = -result; + // result *= sd * constants::root_two(); + // result += mean; + // return result; +// } // quantile + +template +inline RealType cdf(const complemented2_type, RealType>& c) +{ + BOOST_MATH_STD_USING // for ADL of std functions + + 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, &result, Policy())) + return result; + + RealType diff = (x - mean) / (conc * constants::root_two()); + result = boost::math::erfc(diff, Policy()) / 2; + return result; +} // 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; + // result = boost::math::erfc_inv(2 * q, Policy()); + // result *= conc * constants::root_two(); + // result += mean; + // return result; +// } // quantile + +template +inline RealType mean(const von_mises_distribution& dist) +{ + return dist.mean(); +} + +template +inline RealType standard_deviation(const von_mises_distribution& dist) +{ + return dist.concentration(); +} + +template +inline RealType mode(const von_mises_distribution& dist) +{ + return dist.mean(); +} + +template +inline RealType median(const von_mises_distribution& dist) +{ + return dist.mean(); +} + +template +inline RealType skewness(const von_mises_distribution& /*dist*/) +{ + return 0; +} + +// template +// inline RealType kurtosis(const von_mises_distribution& /*dist*/) +// { +// return 3; +// } + +// template +// inline RealType kurtosis_excess(const von_mises_distribution& /*dist*/) +// { + // return 0; +// } + +template +inline RealType entropy(const von_mises_distribution & dist) +{ + using std::log; + RealType arg = constants::two_pi() * cyl_bessel_i(0, dist.concentration(), Policy()); + return log(arg) - dist.concentration() * cyl_bessel_i(1, dist.concentration(), Policy()) / cyl_bessel_i(0, dist.concentration(), Policy()); +} + +} // 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_STATS_VON_MISES_HPP + + diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index a1ce0e7daa..db6c4bc282 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -793,6 +793,7 @@ 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 ] [ 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..65e099a17f --- /dev/null +++ b/test/test_von_mises.cpp @@ -0,0 +1,397 @@ +// 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 + +#include // include directory /libs/math/src/tr1/ is needed. + +#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 "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 + RealType tolerance = 1e-2f; // 1e-4 (as %) + // Some tests only pass at 1e-4 because values generated by + // http://faculty.vassar.edu/lowry/VassarStats.html + // give only 5 or 6 *fixed* places, so small values have fewer digits. + + // 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) + // + 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), 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); + // Tests with larger numbers whose bessel function value exceeds the representable range + BOOST_CHECK_CLOSE( + pdf(von_mises_distribution(3, 25), static_cast(3)), + static_cast(1.98455543847726689510475504795539869409664L), // e^25/(2*pi*I0(25) + tolerance); + BOOST_CHECK_CLOSE( + pdf(von_mises_distribution(3, 125), static_cast(3)), + static_cast(4.45583423571762164664231477348923797078210L), // e^125/(2*pi*I0(125) + tolerance); + BOOST_CHECK_CLOSE( + pdf(von_mises_distribution(3, 1250), static_cast(3)), + static_cast(14.10332862065253024470686549675132142170213L), // e^1250/(2*pi*I0(1250) + tolerance); + + return; + + cout << "Tolerance for type " << typeid(RealType).name() << " is " << tolerance << " %" << endl; + + check_von_mises( + static_cast(5), + static_cast(2), + static_cast(4.8), + static_cast(0.46017), + static_cast(1 - 0.46017), + tolerance); + + check_von_mises( + static_cast(5), + static_cast(2), + static_cast(5.2), + static_cast(1 - 0.46017), + static_cast(0.46017), + tolerance); + + check_von_mises( + static_cast(5), + static_cast(2), + static_cast(2.2), + static_cast(0.08076), + static_cast(1 - 0.08076), + tolerance); + + check_von_mises( + static_cast(5), + static_cast(2), + static_cast(7.8), + static_cast(1 - 0.08076), + static_cast(0.08076), + tolerance); + + check_von_mises( + static_cast(-3), + static_cast(5), + static_cast(-4.5), + static_cast(0.38209), + static_cast(1 - 0.38209), + tolerance); + + check_von_mises( + static_cast(-3), + static_cast(5), + static_cast(-1.5), + static_cast(1 - 0.38209), + static_cast(0.38209), + tolerance); + + check_von_mises( + static_cast(-3), + static_cast(5), + static_cast(-8.5), + static_cast(0.13567), + static_cast(1 - 0.13567), + tolerance); + + check_von_mises( + static_cast(-3), + static_cast(5), + static_cast(2.5), + static_cast(1 - 0.13567), + static_cast(0.13567), + tolerance); + + // + // Spot checks for mean = -5, conc = 6: + // +// for(RealType x = -15; x < 5; x += 0.125) +// { +// BOOST_CHECK_CLOSE( +// pdf(von_mises_distribution(-5, 6), x), +// NaivePDF(RealType(-5), RealType(6), x), +// tolerance); +// } + + RealType tol2 = boost::math::tools::epsilon() * 5; + von_mises_distribution dist(8, 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(8), tol2); + // variance: + BOOST_CHECK_CLOSE( + variance(dist) + , static_cast(9), tol2); + // std deviation: + BOOST_CHECK_CLOSE( + standard_deviation(dist) + , static_cast(3), 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(8), tol2); + + BOOST_CHECK_CLOSE( + median(dist) + , static_cast(8), tol2); + + // skewness: + BOOST_CHECK_CLOSE( + skewness(dist) + , static_cast(0), tol2); + // kurtosis: + // BOOST_CHECK_CLOSE( + // kurtosis(dist) + // , static_cast(3), tol2); + // kurtosis excess: + // BOOST_CHECK_CLOSE( + // kurtosis_excess(dist) + // , static_cast(0), tol2); + + RealType expected_entropy = log(boost::math::constants::two_pi()*boost::math::constants::e()*9)/2; + BOOST_CHECK_CLOSE( + entropy(dist) + ,expected_entropy, 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 + + BOOST_CHECK_CLOSE( + standard_deviation(def_norm01), + static_cast(1), 0); // Mean == zero + + // Error tests: + check_out_of_range >(0, 1); // (All) valid constructor parameter values. + + BOOST_MATH_CHECK_THROW(pdf(von_mises_distribution(0, 0), 0), std::domain_error); + BOOST_MATH_CHECK_THROW(pdf(von_mises_distribution(0, -1), 0), std::domain_error); + 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) + +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 % +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + //test_spots(0.0L); // Test long double. +#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x0582)) + //test_spots(boost::math::concepts::real_concept(0.)); // Test real concept. +#endif +#else + std::cout << "The long double tests have been disabled on this platform " + "either because the long double overloads of the usual math functions are " + "not available at all, or because they are too inaccurate for these tests " + "to pass." << std::endl; +#endif + + +} // BOOST_AUTO_TEST_CASE( test_main ) + +/* + +Output: + +Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_von_mises.exe" +Running 1 test case... +Tolerance for type float is 0.01 % +Tolerance for type double is 0.01 % +Tolerance for type long double is 0.01 % +Tolerance for type class boost::math::concepts::real_concept is 0.01 % +*** No errors detected + + + + +------ Build started: Project: test_von_mises, Configuration: Release Win32 ------ + test_von_mises.cpp + Generating code + Finished generating code + test_von_mises.vcxproj -> J:\Cpp\MathToolkit\test\Math_test\Release\test_von_mises.exe + Running 1 test case... + Tolerance for type float is 0.01 % + Tolerance for type double is 0.01 % + Tolerance for type long double is 0.01 % + Tolerance for type class boost::math::concepts::real_concept is 0.01 % + + *** No errors detected + Detected memory leaks! + Dumping objects -> + {2413} normal block at 0x00321190, 42 bytes long. + Data: 63 6C 61 73 73 20 62 6F 6F 73 74 3A 3A 6D 61 74 + {2412} normal block at 0x003231F0, 8 bytes long. + Data: < 2 22 > 90 11 32 00 98 32 32 00 + {1824} normal block at 0x00323180, 12 bytes long. + Data: 6C 6F 6E 67 20 64 6F 75 62 6C 65 00 + {1823} normal block at 0x00323298, 8 bytes long. + Data: < 12 `22 > 80 31 32 00 60 32 32 00 + {1227} normal block at 0x00323148, 7 bytes long. + Data: 64 6F 75 62 6C 65 00 + {1226} normal block at 0x00323260, 8 bytes long. + Data: 48 31 32 00 A0 30 32 00 + {633} normal block at 0x003230D8, 6 bytes long. + Data: 66 6C 6F 61 74 00 + {632} normal block at 0x003230A0, 8 bytes long. + Data: < 02 > D8 30 32 00 00 00 00 00 + Object dump complete. +========== Build: 1 succeeded, 0 failed, 0 up-to-date, 0 skipped ========== + + +*/ + + From 828f64d7b95de4e00b3142771213ef8c1eecfd54 Mon Sep 17 00:00:00 2001 From: pcjmuenster Date: Sun, 16 Feb 2020 14:18:11 +0100 Subject: [PATCH 2/8] Add support for von Mises CDF Redefine support after discussion on GitHub (#296). Use an established algorithm from the literature for CDF. Add tests for various concentration values. Some tests fail for very large values. --- .../boost/math/distributions/von_mises.hpp | 280 ++++++++++++------ test/test_von_mises.cpp | 246 +++++++-------- 2 files changed, 294 insertions(+), 232 deletions(-) diff --git a/include/boost/math/distributions/von_mises.hpp b/include/boost/math/distributions/von_mises.hpp index f2f8598591..10c9970689 100644 --- a/include/boost/math/distributions/von_mises.hpp +++ b/include/boost/math/distributions/von_mises.hpp @@ -79,7 +79,7 @@ inline const std::pair range(const von_mises_distribution::has_infinity) { - return std::pair(-std::numeric_limits::infinity(), + return std::pair(-std::numeric_limits::infinity(), +std::numeric_limits::infinity()); // - to + infinity. } else @@ -91,17 +91,17 @@ inline const std::pair range(const von_mises_distribution -inline const std::pair support(const von_mises_distribution& /*dist*/) +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(-constants::pi(), - +constants::pi()); // -pi to +pi + return std::pair(dist.mean() - constants::pi(), + dist.mean() + constants::pi()); // [µ-π, µ+π) } #ifdef BOOST_MSVC #pragma warning(pop) #endif -// float version of log_bessel_i0 +// float version of pdf_impl template inline RealType pdf_impl(const von_mises_distribution& dist, RealType x, const boost::integral_constant&) @@ -117,8 +117,9 @@ inline RealType pdf_impl(const von_mises_distribution& dist, R } else if(conc < 50) { - // Max error in interpolated form: 5.195e-08 - // Max Error found at float precision = Poly: 8.502534e-08 + // Polynomial coefficients from boost/math/special_functions/detail/bessel_i0.hpp + // > Max error in interpolated form: 5.195e-08 + // > Max Error found at float precision = Poly: 8.502534e-08 static const float P[] = { 3.98942651588301770e-01f, 4.98327234176892844e-02f, @@ -126,14 +127,16 @@ inline RealType pdf_impl(const von_mises_distribution& dist, R 1.35614940793742178e-02f, 1.31409251787866793e-01f }; - RealType exponent = conc * (cos(x - mean) - 1.f); - exponent -= std::log(boost::math::tools::evaluate_polynomial(P, RealType(1 / conc)) / sqrt(conc)); - return exp(exponent) / boost::math::constants::two_pi(); + 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; } else { - // Max error in interpolated form: 1.782e-09 - // Max Error found at float precision = Poly: 6.473568e-08 + // Polynomial coefficients from boost/math/special_functions/detail/bessel_i0.hpp + // > Max error in interpolated form: 1.782e-09 + // > Max Error found at float precision = Poly: 6.473568e-08 static const float P[] = { 3.98942280401432677e-01f, 4.98677850501790847e-02f, @@ -141,42 +144,33 @@ inline RealType pdf_impl(const von_mises_distribution& dist, R 2.92194053028393074e-02f, 4.47422143699726895e-02f, }; - RealType exponent = conc * (cos(x - mean) - 1.f); - exponent -= std::log(boost::math::tools::evaluate_polynomial(P, RealType(1 / conc)) / sqrt(conc)); - return exp(exponent) / boost::math::constants::two_pi(); + 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 log_bessel_i0 -template -inline RealType log_bessel_i0_impl(RealType x, const boost::integral_constant&) +// double version of pdf_impl +template +inline RealType pdf_impl(const von_mises_distribution& dist, RealType x, + const boost::integral_constant&) { - BOOST_MATH_STD_USING - if(x < 7.75) - { - static const double P[] = { - 1.00000000000000000e+00, - 2.49999999999999909e-01, - 2.77777777777782257e-02, - 1.73611111111023792e-03, - 6.94444444453352521e-05, - 1.92901234513219920e-06, - 3.93675991102510739e-08, - 6.15118672704439289e-10, - 7.59407002058973446e-12, - 7.59389793369836367e-14, - 6.27767773636292611e-16, - 4.34709704153272287e-18, - 2.63417742690109154e-20, - 1.13943037744822825e-22, - 9.07926920085624812e-25 - }; - RealType a = x * x / 4; - return std::log(a * boost::math::tools::evaluate_polynomial(P, a) + 1); - } - else if(x < 500) - { - static const double P[] = { + RealType mean = dist.mean(); + RealType conc = dist.concentration(); + + BOOST_MATH_STD_USING + if(conc < 7.75) + { + RealType bessel_i0 = cyl_bessel_i(0, conc, Policy()); + return exp(conc * cos(x - mean)) / bessel_i0 / boost::math::constants::two_pi(); + } + else if(conc < 500) + { + // Polynomial coefficients from boost/math/special_functions/detail/bessel_i0.hpp + // > Max error in interpolated form : 1.685e-16 + // > Max Error found at double precision = Poly : 2.575063e-16 Cheb : 2.247615e+00 + static const double P[] = { 3.98942280401425088e-01, 4.98677850604961985e-02, 2.80506233928312623e-02, @@ -200,41 +194,103 @@ inline RealType log_bessel_i0_impl(RealType x, const boost::integral_constant(); + return result; + } + else + { + // Polynomial coefficients from boost/math/special_functions/detail/bessel_i0.hpp + // > Max error in interpolated form : 2.437e-18 + // > Max Error found at double precision = Poly : 1.216719e-16 + 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; + } } -template -inline RealType log_bessel_i0(RealType x) +// long double version of pdf_impl +template +inline RealType pdf_impl(const von_mises_distribution& dist, RealType x, + const boost::integral_constant&) { - 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 log_bessel_i0_impl(x, tag_type()); + RealType mean = dist.mean(); + RealType conc = dist.concentration(); + + BOOST_MATH_STD_USING + if(conc < 15) + { + RealType bessel_i0 = cyl_bessel_i(0, conc, Policy()); + return exp(conc * cos(x - mean)) / bessel_i0 / boost::math::constants::two_pi(); + } + else if(x < 50) + { + // Max error in interpolated form: 1.035e-21 + // Max Error found at float80 precision = Poly: 1.885872e-21 + static const RealType Y = 4.011702537536621093750e-01f; + static const RealType P[] = { + BOOST_MATH_BIG_CONSTANT(RealType, 64, -2.227973351806078464328e-03), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.986778486088017419036e-02), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 2.805066823812285310011e-02), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 2.921443721160964964623e-02), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.517504941996594744052e-02), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 6.316922639868793684401e-02), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.535891099168810015433e+00), + BOOST_MATH_BIG_CONSTANT(RealType, 64, -4.706078229522448308087e+01), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.351015763079160914632e+03), + BOOST_MATH_BIG_CONSTANT(RealType, 64, -2.948809013999277355098e+04), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.967598958582595361757e+05), + BOOST_MATH_BIG_CONSTANT(RealType, 64, -6.346924657995383019558e+06), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 5.998794574259956613472e+07), + BOOST_MATH_BIG_CONSTANT(RealType, 64, -4.016371355801690142095e+08), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.768791455631826490838e+09), + BOOST_MATH_BIG_CONSTANT(RealType, 64, -4.441995678177349895640e+09), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.482292669974971387738e+09) + }; + RealType result = exp(conc * (cos(x - mean) - 1.0)); + result /= (boost::math::tools::evaluate_polynomial(P, RealType(1.0 / conc)) + Y) / sqrt(conc); + * boost::math::constants::two_pi(); + return result; + } + 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 @@ -256,7 +312,7 @@ inline RealType pdf(const von_mises_distribution& dist, const { return result; } - if(false == detail::check_x(function, x, &result, Policy())) + if(false == detail::check_angle(function, x - mean, &result, Policy())) { return result; } @@ -284,28 +340,58 @@ inline RealType pdf(const von_mises_distribution& dist, const template inline RealType cdf(const von_mises_distribution& dist, const RealType& x) { - BOOST_MATH_STD_USING // for ADL of std functions + BOOST_MATH_STD_USING // for ADL of std functions - 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, &result, Policy())) - { - return result; - } + 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; + } - RealType diff = (x - mean) / (conc * constants::root_two()); - result = boost::math::erfc(-diff, Policy()) / 2; - return result; + RealType u = x - mean; + if (u <= -boost::math::constants::pi()) + return 0; + if (u >= +boost::math::constants::pi()) + return 1; + + if (conc > RealType{10.5}) { + auto c = 24.0 * conc; + auto v = c - 56.0; + auto r = sqrt((54 / (347 / v + 26 - c) - 6 + c) / 12); + auto z = sin(u / 2) * r; + auto s = z * z * 2; + v = v - s + 3; + auto y = (c - s - s - 16) / 3; + 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) { + int iterations = static_cast(ceil(conc * 0.8 - 8 / (conc + 1) + 12)); + 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) / 2 + v) / boost::math::constants::pi() + 0.5; + } + + return result <= 0 ? 0 : (1 <= result ? 1 : result); } // cdf // template diff --git a/test/test_von_mises.cpp b/test/test_von_mises.cpp index 65e099a17f..6fdd7ada12 100644 --- a/test/test_von_mises.cpp +++ b/test/test_von_mises.cpp @@ -48,37 +48,32 @@ void check_von_mises(RealType mean, RealType conc, RealType x, RealType p, RealT 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. + //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 - RealType tolerance = 1e-2f; // 1e-4 (as %) - // Some tests only pass at 1e-4 because values generated by - // http://faculty.vassar.edu/lowry/VassarStats.html - // give only 5 or 6 *fixed* places, so small values have fewer digits. - + // 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 @@ -90,10 +85,10 @@ void test_spots(RealType) 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(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 @@ -121,111 +116,133 @@ void test_spots(RealType) // // Tests for PDF: we know that the peak value is at e/(2*pi*I0(1) // - tolerance = boost::math::tools::epsilon() * 5 * 100; // 5 eps as a percentage + 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) + 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) + 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) + static_cast(0.86713652854235200257351846969777045343907L), // e^5/(2*pi*I0(5)) tolerance); // Tests with larger numbers whose bessel function value exceeds the representable range BOOST_CHECK_CLOSE( pdf(von_mises_distribution(3, 25), static_cast(3)), - static_cast(1.98455543847726689510475504795539869409664L), // e^25/(2*pi*I0(25) + static_cast(1.98455543847726689510475504795539869409664L), // e^25/(2*pi*I0(25)) tolerance); BOOST_CHECK_CLOSE( pdf(von_mises_distribution(3, 125), static_cast(3)), - static_cast(4.45583423571762164664231477348923797078210L), // e^125/(2*pi*I0(125) + static_cast(4.45583423571762164664231477348923797078210L), // e^125/(2*pi*I0(125)) tolerance); BOOST_CHECK_CLOSE( pdf(von_mises_distribution(3, 1250), static_cast(3)), - static_cast(14.10332862065253024470686549675132142170213L), // e^1250/(2*pi*I0(1250) + static_cast(14.1033286206525302447068654967513214217021L), // e^1250/(2*pi*I0(1250)) tolerance); - return; + tolerance = 1e-3f; // 1e-5 (as %) + // Some tests only pass at 1e-5 because values generated by cout << "Tolerance for type " << typeid(RealType).name() << " is " << tolerance << " %" << endl; + // test CDF for mean and interval edges check_von_mises( - static_cast(5), - static_cast(2), - static_cast(4.8), - static_cast(0.46017), - static_cast(1 - 0.46017), + static_cast(0), + static_cast(1), + static_cast(0), + static_cast(0.5L), + static_cast(0.5L), tolerance); - + check_von_mises( - static_cast(5), - static_cast(2), - static_cast(5.2), - static_cast(1 - 0.46017), - static_cast(0.46017), + 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(5), - static_cast(2), - static_cast(2.2), - static_cast(0.08076), - static_cast(1 - 0.08076), + 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(5), - static_cast(2), - static_cast(7.8), - static_cast(1 - 0.08076), - static_cast(0.08076), + 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(-3), + static_cast(0), static_cast(5), - static_cast(-4.5), - static_cast(0.38209), - static_cast(1 - 0.38209), - tolerance); - + static_cast(1), + static_cast(0.98096204546814689054581384796251763480020360394758184271L), + static_cast(0.01903795453185310945418615203748236519979639605241815729L), + tolerance); + check_von_mises( - static_cast(-3), + static_cast(0), static_cast(5), - static_cast(-1.5), - static_cast(1 - 0.38209), - static_cast(0.38209), + static_cast(-1), + static_cast(0.01903795453185310945418615203748236519979639605241815729L), + static_cast(0.98096204546814689054581384796251763480020360394758184271L), tolerance); - + + // test CDF for high concentrations check_von_mises( - static_cast(-3), - static_cast(5), - static_cast(-8.5), - static_cast(0.13567), - static_cast(1 - 0.13567), + 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(-3), - static_cast(5), - static_cast(2.5), - static_cast(1 - 0.13567), - static_cast(0.13567), + 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); - - // - // Spot checks for mean = -5, conc = 6: - // -// for(RealType x = -15; x < 5; x += 0.125) -// { -// BOOST_CHECK_CLOSE( -// pdf(von_mises_distribution(-5, 6), x), -// NaivePDF(RealType(-5), RealType(6), x), -// tolerance); -// } + return; RealType tol2 = boost::math::tools::epsilon() * 5; von_mises_distribution dist(8, 3); @@ -327,7 +344,7 @@ BOOST_AUTO_TEST_CASE( test_main ) // 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.0); // Test double. OK at decdigits 7, tolerance = 1e07 % #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS //test_spots(0.0L); // Test long double. #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x0582)) @@ -353,45 +370,4 @@ Tolerance for type float is 0.01 % Tolerance for type double is 0.01 % Tolerance for type long double is 0.01 % Tolerance for type class boost::math::concepts::real_concept is 0.01 % -*** No errors detected - - - - ------- Build started: Project: test_von_mises, Configuration: Release Win32 ------ - test_von_mises.cpp - Generating code - Finished generating code - test_von_mises.vcxproj -> J:\Cpp\MathToolkit\test\Math_test\Release\test_von_mises.exe - Running 1 test case... - Tolerance for type float is 0.01 % - Tolerance for type double is 0.01 % - Tolerance for type long double is 0.01 % - Tolerance for type class boost::math::concepts::real_concept is 0.01 % - - *** No errors detected - Detected memory leaks! - Dumping objects -> - {2413} normal block at 0x00321190, 42 bytes long. - Data: 63 6C 61 73 73 20 62 6F 6F 73 74 3A 3A 6D 61 74 - {2412} normal block at 0x003231F0, 8 bytes long. - Data: < 2 22 > 90 11 32 00 98 32 32 00 - {1824} normal block at 0x00323180, 12 bytes long. - Data: 6C 6F 6E 67 20 64 6F 75 62 6C 65 00 - {1823} normal block at 0x00323298, 8 bytes long. - Data: < 12 `22 > 80 31 32 00 60 32 32 00 - {1227} normal block at 0x00323148, 7 bytes long. - Data: 64 6F 75 62 6C 65 00 - {1226} normal block at 0x00323260, 8 bytes long. - Data: 48 31 32 00 A0 30 32 00 - {633} normal block at 0x003230D8, 6 bytes long. - Data: 66 6C 6F 61 74 00 - {632} normal block at 0x003230A0, 8 bytes long. - Data: < 02 > D8 30 32 00 00 00 00 00 - Object dump complete. -========== Build: 1 succeeded, 0 failed, 0 up-to-date, 0 skipped ========== - - -*/ - - +*** No errors detected */ From a59c12ea824bb12d7a0cbff262fc3f05da41fcd5 Mon Sep 17 00:00:00 2001 From: pcjmuenster Date: Sun, 23 Feb 2020 16:28:09 +0100 Subject: [PATCH 3/8] Add quantile function for von Mises distribution MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Use Newton–Raphson iteration to calculate zeros of CDF - quantile. --- include/boost/math/distributions/fwd.hpp | 4 + .../boost/math/distributions/von_mises.hpp | 385 ++++++++++-------- test/test_von_mises.cpp | 64 +-- 3 files changed, 242 insertions(+), 211 deletions(-) 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 index 10c9970689..3ce27c8325 100644 --- a/include/boost/math/distributions/von_mises.hpp +++ b/include/boost/math/distributions/von_mises.hpp @@ -13,10 +13,11 @@ // http://mathworld.wolfram.com/VonMisesDistribution.html #include -#include -#include // for erf/erfc. #include #include +#include // for besseli0 and besseli1 +#include // for erf +#include #include @@ -26,45 +27,45 @@ 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; - } + 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. + // + // 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; @@ -93,14 +94,15 @@ inline const std::pair range(const von_mises_distribution 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()); // [µ-π, µ+π) + 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, @@ -130,7 +132,7 @@ inline RealType pdf_impl(const von_mises_distribution& dist, R 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; + return result; } else { @@ -147,7 +149,7 @@ inline RealType pdf_impl(const von_mises_distribution& dist, R 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; + return result; } } @@ -195,8 +197,8 @@ inline RealType pdf_impl(const von_mises_distribution& dist, R 2.17587543863819074e+15 }; 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(); + result /= boost::math::tools::evaluate_polynomial(P, RealType(1.0 / conc)) / sqrt(conc) + * boost::math::constants::two_pi(); return result; } else @@ -227,34 +229,34 @@ inline RealType pdf_impl(const von_mises_distribution& dist, R RealType conc = dist.concentration(); BOOST_MATH_STD_USING - if(conc < 15) + if (conc < 15) { - RealType bessel_i0 = cyl_bessel_i(0, conc, Policy()); - return exp(conc * cos(x - mean)) / bessel_i0 / boost::math::constants::two_pi(); + RealType bessel_i0 = cyl_bessel_i(0, conc, Policy()); + return exp(conc * cos(x - mean)) / bessel_i0 / boost::math::constants::two_pi(); } - else if(x < 50) + else if (x < 50) { // Max error in interpolated form: 1.035e-21 // Max Error found at float80 precision = Poly: 1.885872e-21 static const RealType Y = 4.011702537536621093750e-01f; static const RealType P[] = { - BOOST_MATH_BIG_CONSTANT(RealType, 64, -2.227973351806078464328e-03), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.986778486088017419036e-02), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 2.805066823812285310011e-02), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 2.921443721160964964623e-02), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.517504941996594744052e-02), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 6.316922639868793684401e-02), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.535891099168810015433e+00), - BOOST_MATH_BIG_CONSTANT(RealType, 64, -4.706078229522448308087e+01), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.351015763079160914632e+03), - BOOST_MATH_BIG_CONSTANT(RealType, 64, -2.948809013999277355098e+04), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.967598958582595361757e+05), - BOOST_MATH_BIG_CONSTANT(RealType, 64, -6.346924657995383019558e+06), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 5.998794574259956613472e+07), - BOOST_MATH_BIG_CONSTANT(RealType, 64, -4.016371355801690142095e+08), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.768791455631826490838e+09), - BOOST_MATH_BIG_CONSTANT(RealType, 64, -4.441995678177349895640e+09), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.482292669974971387738e+09) + BOOST_MATH_BIG_CONSTANT(RealType, 64, -2.227973351806078464328e-03), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.986778486088017419036e-02), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 2.805066823812285310011e-02), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 2.921443721160964964623e-02), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.517504941996594744052e-02), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 6.316922639868793684401e-02), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.535891099168810015433e+00), + BOOST_MATH_BIG_CONSTANT(RealType, 64, -4.706078229522448308087e+01), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.351015763079160914632e+03), + BOOST_MATH_BIG_CONSTANT(RealType, 64, -2.948809013999277355098e+04), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.967598958582595361757e+05), + BOOST_MATH_BIG_CONSTANT(RealType, 64, -6.346924657995383019558e+06), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 5.998794574259956613472e+07), + BOOST_MATH_BIG_CONSTANT(RealType, 64, -4.016371355801690142095e+08), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.768791455631826490838e+09), + BOOST_MATH_BIG_CONSTANT(RealType, 64, -4.441995678177349895640e+09), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.482292669974971387738e+09) }; RealType result = exp(conc * (cos(x - mean) - 1.0)); result /= (boost::math::tools::evaluate_polynomial(P, RealType(1.0 / conc)) + Y) / sqrt(conc); @@ -267,24 +269,24 @@ inline RealType pdf_impl(const von_mises_distribution& dist, R // 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) + 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) @@ -292,79 +294,66 @@ inline RealType pdf_impl(const von_mises_distribution& dist, R return result; } } +} // 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 pdf_impl(dist, x, tag_type()); -} // pdf - -template -inline RealType cdf(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::cdf(const von_mises_distribution<%1%>&, %1%)"; + + 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())) - { + if (false == detail::check_positive_x(function, conc, &result, Policy())) return result; - } - if(false == detail::check_angle(function, mean, &result, Policy())) - { + if (false == detail::check_angle(function, mean, &result, Policy())) return result; - } - if(false == detail::check_angle(function, x - mean, &result, Policy())) - { + 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; if (conc > RealType{10.5}) { auto c = 24.0 * conc; auto v = c - 56.0; @@ -388,57 +377,95 @@ inline RealType cdf(const von_mises_distribution& dist, const v = (sj / j + v) * r; } } - result = ((x - mean) / 2 + v) / boost::math::constants::pi() + 0.5; + result = (x - mean + boost::math::constants::pi()) / 2; + result = (result + v) / boost::math::constants::pi(); } return result <= 0 ? 0 : (1 <= result ? 1 : result); -} // cdf +} +} // namespace detail -// template -// inline RealType quantile(const von_mises_distribution& dist, const RealType& p) -// { - // BOOST_MATH_STD_USING // for ADL of std functions +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 - // RealType conc = dist.concentration(); - // RealType mean = dist.mean(); - // static const char* function = "boost::math::quantile(const von_mises_distribution<%1%>&, %1%)"; +template +inline RealType quantile(const von_mises_distribution& dist, const RealType& p) +{ + BOOST_MATH_STD_USING // for ADL of std functions - // 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; + RealType conc = dist.concentration(); + RealType mean = dist.mean(); + static const char* function = "boost::math::quantile(const von_mises_distribution<%1%>&, %1%)"; - // result = boost::math::erfc_inv(2 * p, Policy()); - // result = -result; - // result *= sd * constants::root_two(); - // result += mean; - // return result; -// } // quantile + 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; + + auto step_func = [&](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, mean, lower, upper, 15 /* digits */); + + return zero; +} // quantile template inline RealType cdf(const complemented2_type, RealType>& c) { - BOOST_MATH_STD_USING // for ADL of std functions - - 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, &result, Policy())) - return result; - - RealType diff = (x - mean) / (conc * constants::root_two()); - result = boost::math::erfc(diff, Policy()) / 2; - return result; + 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 diff --git a/test/test_von_mises.cpp b/test/test_von_mises.cpp index 6fdd7ada12..d8a2c0e95e 100644 --- a/test/test_von_mises.cpp +++ b/test/test_von_mises.cpp @@ -42,25 +42,25 @@ template void check_von_mises(RealType mean, RealType conc, RealType x, RealType p, RealType q, RealType tol) { - BOOST_CHECK_CLOSE( + 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. + 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( @@ -227,21 +227,21 @@ void test_spots(RealType) 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(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); + //~ check_von_mises( + //~ static_cast(0), + //~ static_cast(125), + //~ static_cast(-1), + //~ static_cast(3.3546365686506670890489976669185096109e-20L), + //~ static_cast(0.99999999999999999996645363431349332910951002333081490389L), + //~ tolerance); return; RealType tol2 = boost::math::tools::epsilon() * 5; From f76081c2bde791b7bd98b674bfb1de806fdfdff8 Mon Sep 17 00:00:00 2001 From: pcjmuenster Date: Sun, 15 Mar 2020 16:47:56 +0100 Subject: [PATCH 4/8] [von_mises] add complementary quantile and variance --- .../boost/math/distributions/von_mises.hpp | 368 ++++++++++++++++-- test/test_von_mises.cpp | 78 ++-- 2 files changed, 372 insertions(+), 74 deletions(-) diff --git a/include/boost/math/distributions/von_mises.hpp b/include/boost/math/distributions/von_mises.hpp index 3ce27c8325..d1ea74a406 100644 --- a/include/boost/math/distributions/von_mises.hpp +++ b/include/boost/math/distributions/von_mises.hpp @@ -354,21 +354,30 @@ inline RealType cdf_impl(const von_mises_distribution& dist, c // "Algorithm 518: Incomplete Bessel Function I0. The Von Mises Distribution", 1977, ACM // DOI: 10.1145/355744.355753 RealType result = 0; - if (conc > RealType{10.5}) { - auto c = 24.0 * conc; - auto v = c - 56.0; - auto r = sqrt((54 / (347 / v + 26 - c) - 6 + c) / 12); - auto z = sin(u / 2) * r; - auto s = z * z * 2; + + 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; - auto y = (c - s - s - 16) / 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) { - int iterations = static_cast(ceil(conc * 0.8 - 8 / (conc + 1) + 12)); + // 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) { @@ -409,7 +418,8 @@ inline RealType quantile(const von_mises_distribution& dist, c RealType conc = dist.concentration(); RealType mean = dist.mean(); - static const char* function = "boost::math::quantile(const von_mises_distribution<%1%>&, %1%)"; + 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())) @@ -425,17 +435,14 @@ inline RealType quantile(const von_mises_distribution& dist, c 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; + ((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 + > tag_type; auto step_func = [&](RealType x) { return std::make_pair(detail::cdf_impl(dist, x) - p, // f(x) @@ -455,7 +462,8 @@ inline RealType cdf(const complemented2_type -// 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; - // result = boost::math::erfc_inv(2 * q, Policy()); - // result *= conc * constants::root_two(); - // result += mean; - // return result; -// } // quantile +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; + + auto step_func = [&](RealType x) { + RealType xc = 2 * mean - x; + return std::make_pair(detail::cdf_impl(c.dist, xc) - q , // 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, mean, lower, upper, 15 /* digits */); + + return zero; +} // quantile template inline RealType mean(const von_mises_distribution& dist) @@ -513,6 +546,257 @@ 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(RealType conc, const boost::integral_constant&) +//~ { + //~ BOOST_MATH_STD_USING + //~ if (conc < 15) + //~ { + //~ 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 (x < 50) + //~ { + //~ // Max error in interpolated form: 1.035e-21 + //~ // Max Error found at float80 precision = Poly: 1.885872e-21 + //~ static const RealType Y = 4.011702537536621093750e-01f; + //~ static const RealType P[] = { + //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, -2.227973351806078464328e-03), + //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.986778486088017419036e-02), + //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, 2.805066823812285310011e-02), + //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, 2.921443721160964964623e-02), + //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.517504941996594744052e-02), + //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, 6.316922639868793684401e-02), + //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.535891099168810015433e+00), + //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, -4.706078229522448308087e+01), + //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.351015763079160914632e+03), + //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, -2.948809013999277355098e+04), + //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.967598958582595361757e+05), + //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, -6.346924657995383019558e+06), + //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, 5.998794574259956613472e+07), + //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, -4.016371355801690142095e+08), + //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.768791455631826490838e+09), + //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, -4.441995678177349895640e+09), + //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.482292669974971387738e+09) + //~ }; + //~ RealType result = exp(conc * (cos(x - mean) - 1.0)); + //~ result /= (boost::math::tools::evaluate_polynomial(P, RealType(1.0 / conc)) + Y) / sqrt(conc); + //~ * boost::math::constants::two_pi(); + //~ return result; + //~ } + //~ 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; + //~ } +//~ } +} // 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*/) diff --git a/test/test_von_mises.cpp b/test/test_von_mises.cpp index d8a2c0e95e..cd2e146c7d 100644 --- a/test/test_von_mises.cpp +++ b/test/test_von_mises.cpp @@ -61,13 +61,13 @@ void check_von_mises(RealType mean, RealType conc, RealType x, RealType p, RealT 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. + BOOST_CHECK_CLOSE( + ::boost::math::quantile( + complement( + von_mises_distribution(mean, conc), // distribution. + q)), // probability complement. + x, // random variable. + tol); // %tolerance. } template @@ -107,10 +107,10 @@ void test_spots(RealType) { // 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 + 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 } // @@ -211,21 +211,21 @@ void test_spots(RealType) 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(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(25), + //~ static_cast(-1), + //~ static_cast(9.37595535559547766495510700217223833735532789427234e-7L), + //~ static_cast(0.999999062404464440452233504489299782776166264467210572765L), + //~ tolerance); //~ check_von_mises( //~ static_cast(0), @@ -242,10 +242,10 @@ void test_spots(RealType) //~ static_cast(3.3546365686506670890489976669185096109e-20L), //~ static_cast(0.99999999999999999996645363431349332910951002333081490389L), //~ tolerance); - return; + - RealType tol2 = boost::math::tools::epsilon() * 5; - von_mises_distribution dist(8, 3); + 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 @@ -253,11 +253,25 @@ void test_spots(RealType) // mean: BOOST_CHECK_CLOSE( mean(dist) - , static_cast(8), tol2); - // variance: + , static_cast(2), tol2); + // variance: BOOST_CHECK_CLOSE( - variance(dist) - , static_cast(9), tol2); + 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); + return; + // std deviation: BOOST_CHECK_CLOSE( standard_deviation(dist) From a385e69ff358e1e0b6af0357b1ea8b7e916a50ad Mon Sep 17 00:00:00 2001 From: pcjmuenster Date: Fri, 27 Mar 2020 15:00:31 +0100 Subject: [PATCH 5/8] clean up code before PR --- include/boost/math/distributions/von_mises.hpp | 9 +++++++-- test/test_von_mises.cpp | 16 ++++------------ 2 files changed, 11 insertions(+), 14 deletions(-) diff --git a/include/boost/math/distributions/von_mises.hpp b/include/boost/math/distributions/von_mises.hpp index d1ea74a406..eecb57c3e5 100644 --- a/include/boost/math/distributions/von_mises.hpp +++ b/include/boost/math/distributions/von_mises.hpp @@ -532,7 +532,9 @@ inline RealType mean(const von_mises_distribution& dist) template inline RealType standard_deviation(const von_mises_distribution& dist) { - return dist.concentration(); + RealType bessel_quot = cyl_bessel_i(1, dist.concentration(), Policy()) + / cyl_bessel_i(0, dist.concentration(), Policy()); + return std::sqrt(-2 * std::log(bessel_quot)); } template @@ -546,6 +548,7 @@ inline RealType median(const von_mises_distribution& dist) { return dist.mean(); } + namespace detail { // float version of variance_impl template @@ -821,7 +824,9 @@ inline RealType entropy(const von_mises_distribution & dist) { using std::log; RealType arg = constants::two_pi() * cyl_bessel_i(0, dist.concentration(), Policy()); - return log(arg) - dist.concentration() * cyl_bessel_i(1, dist.concentration(), Policy()) / 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 diff --git a/test/test_von_mises.cpp b/test/test_von_mises.cpp index cd2e146c7d..cbcdaf39b3 100644 --- a/test/test_von_mises.cpp +++ b/test/test_von_mises.cpp @@ -270,12 +270,11 @@ void test_spots(RealType) //~ BOOST_CHECK_CLOSE( //~ variance(von_mises_distribution(2, 125)) //~ , static_cast(0.004008064813593637057963834031301840442156903531539609019L), tol2); - return; // std deviation: BOOST_CHECK_CLOSE( standard_deviation(dist) - , static_cast(3), tol2); + , static_cast(0.649213658343262740252572725779410261163523458085389547123L), tol2); // hazard: BOOST_CHECK_CLOSE( hazard(dist, x) @@ -291,11 +290,11 @@ void test_spots(RealType) // mode: BOOST_CHECK_CLOSE( mode(dist) - , static_cast(8), tol2); + , static_cast(2), tol2); BOOST_CHECK_CLOSE( median(dist) - , static_cast(8), tol2); + , static_cast(2), tol2); // skewness: BOOST_CHECK_CLOSE( @@ -310,10 +309,9 @@ void test_spots(RealType) // kurtosis_excess(dist) // , static_cast(0), tol2); - RealType expected_entropy = log(boost::math::constants::two_pi()*boost::math::constants::e()*9)/2; BOOST_CHECK_CLOSE( entropy(dist) - ,expected_entropy, tol2); + , 0.993228806353252817873771736349142645476889275244864748502L, tol2); von_mises_distribution norm01(0, 1); // Test default (0, 1) BOOST_CHECK_CLOSE( @@ -330,15 +328,9 @@ void test_spots(RealType) mean(def_norm01), static_cast(0), 0); // Mean == zero - BOOST_CHECK_CLOSE( - standard_deviation(def_norm01), - static_cast(1), 0); // Mean == zero - // Error tests: check_out_of_range >(0, 1); // (All) valid constructor parameter values. - BOOST_MATH_CHECK_THROW(pdf(von_mises_distribution(0, 0), 0), std::domain_error); - BOOST_MATH_CHECK_THROW(pdf(von_mises_distribution(0, -1), 0), std::domain_error); 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) From f3d324c550f6cabc9b5fdadfdfed9586df0355d7 Mon Sep 17 00:00:00 2001 From: pcjmuenster Date: Wed, 13 May 2020 18:14:19 +0200 Subject: [PATCH 6/8] Reduce the ULP error of PDF and CDF, minor corrections --- .../boost/math/distributions/von_mises.hpp | 349 +++++++----------- test/test_von_mises.cpp | 106 ++++-- 2 files changed, 207 insertions(+), 248 deletions(-) diff --git a/include/boost/math/distributions/von_mises.hpp b/include/boost/math/distributions/von_mises.hpp index eecb57c3e5..da218eada4 100644 --- a/include/boost/math/distributions/von_mises.hpp +++ b/include/boost/math/distributions/von_mises.hpp @@ -1,12 +1,13 @@ // 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_STATS_VON_MISES_HPP -#define BOOST_STATS_VON_MISES_HPP +#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. @@ -15,8 +16,8 @@ #include #include #include -#include // for besseli0 and besseli1 -#include // for erf +#include // for besseli0 and besseli1 +#include // for erf #include #include @@ -64,7 +65,7 @@ class von_mises_distribution // // Data members: // - RealType m_mean; // distribution mean or location. + RealType m_mean; // distribution mean or location. RealType m_concentration; // distribution standard deviation or scale. }; // class von_mises_distribution @@ -79,9 +80,9 @@ 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. + +std::numeric_limits::infinity()); // - to + infinity. } else { // Can only use max_value. @@ -112,33 +113,15 @@ inline RealType pdf_impl(const von_mises_distribution& dist, R RealType conc = dist.concentration(); BOOST_MATH_STD_USING - if(conc < 7.75) + 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 if(conc < 50) - { - // Polynomial coefficients from boost/math/special_functions/detail/bessel_i0.hpp - // > Max error in interpolated form: 5.195e-08 - // > Max Error found at float precision = Poly: 8.502534e-08 - static const float P[] = { - 3.98942651588301770e-01f, - 4.98327234176892844e-02f, - 2.91866904423115499e-02f, - 1.35614940793742178e-02f, - 1.31409251787866793e-01f - }; - 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; - } - else + else // exp(88) > MAX_FLOAT { - // Polynomial coefficients from boost/math/special_functions/detail/bessel_i0.hpp - // > Max error in interpolated form: 1.782e-09 - // > Max Error found at float precision = Poly: 6.473568e-08 + // 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, @@ -162,50 +145,15 @@ inline RealType pdf_impl(const von_mises_distribution& dist, R RealType conc = dist.concentration(); BOOST_MATH_STD_USING - if(conc < 7.75) + 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 if(conc < 500) - { - // Polynomial coefficients from boost/math/special_functions/detail/bessel_i0.hpp - // > Max error in interpolated form : 1.685e-16 - // > Max Error found at double precision = Poly : 2.575063e-16 Cheb : 2.247615e+00 - static const double P[] = { - 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 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; - } - else + else // exp(709) > MAX_DOUBLE { - // Polynomial coefficients from boost/math/special_functions/detail/bessel_i0.hpp - // > Max error in interpolated form : 2.437e-18 - // > Max Error found at double precision = Poly : 1.216719e-16 + // 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, @@ -234,7 +182,7 @@ inline RealType pdf_impl(const von_mises_distribution& dist, R RealType bessel_i0 = cyl_bessel_i(0, conc, Policy()); return exp(conc * cos(x - mean)) / bessel_i0 / boost::math::constants::two_pi(); } - else if (x < 50) + else if (conc < 50) { // Max error in interpolated form: 1.035e-21 // Max Error found at float80 precision = Poly: 1.885872e-21 @@ -259,7 +207,7 @@ inline RealType pdf_impl(const von_mises_distribution& dist, R BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.482292669974971387738e+09) }; RealType result = exp(conc * (cos(x - mean) - 1.0)); - result /= (boost::math::tools::evaluate_polynomial(P, RealType(1.0 / conc)) + Y) / sqrt(conc); + result /= (boost::math::tools::evaluate_polynomial(P, RealType(1.0 / conc)) + Y) / sqrt(conc) * boost::math::constants::two_pi(); return result; } @@ -303,9 +251,9 @@ inline RealType pdf(const von_mises_distribution& dist, const 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; @@ -337,30 +285,30 @@ inline RealType pdf(const von_mises_distribution& dist, const namespace detail { template -inline RealType cdf_impl(const von_mises_distribution& dist, const RealType& x) +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; + 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 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; @@ -371,16 +319,16 @@ inline RealType cdf_impl(const von_mises_distribution& dist, c else { RealType v = 0; if(conc > 0) { - // extrapolation of the tables given in the paper + // 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) { + for (int j = iterations - 1; j > 0; --j) { RealType sj = sin(j * u); r = 1 / (j * z + r); v = (sj / j + v) * r; @@ -389,7 +337,7 @@ inline RealType cdf_impl(const von_mises_distribution& dist, c 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 @@ -407,7 +355,7 @@ inline RealType cdf(const von_mises_distribution& dist, const return result; if (false == detail::check_angle(function, x - mean, &result, Policy())) return result; - + return detail::cdf_impl(dist, x); } // cdf @@ -418,7 +366,7 @@ inline RealType quantile(const von_mises_distribution& dist, c RealType conc = dist.concentration(); RealType mean = dist.mean(); - static const char* function + static const char* function = "boost::math::quantile(const von_mises_distribution<%1%>&, %1%)"; RealType result = 0; @@ -428,31 +376,35 @@ inline RealType quantile(const von_mises_distribution& dist, c 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::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 > tag_type; - - auto step_func = [&](RealType x) { + + 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, mean, lower, upper, 15 /* digits */); - + step_func{dist, p}, mean, lower, upper, 15 /* digits */); + return zero; } // quantile @@ -462,7 +414,7 @@ inline RealType cdf(const complemented2_type(); if (q >= 1) - return -boost::math::constants::pi(); - + return -boost::math::constants::pi(); + typedef boost::integral_constant::digits == 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 > tag_type; - - auto step_func = [&](RealType x) { - RealType xc = 2 * mean - x; - return std::make_pair(detail::cdf_impl(c.dist, xc) - q , // f(x) - -detail::pdf_impl(c.dist, xc, tag_type())); // f'(x) + + 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, mean, lower, upper, 15 /* digits */); - + step_func{c}, mean, lower, upper, 15 /* digits */); + return zero; } // quantile @@ -532,7 +487,7 @@ inline RealType mean(const von_mises_distribution& dist) template inline RealType standard_deviation(const von_mises_distribution& dist) { - RealType bessel_quot = cyl_bessel_i(1, dist.concentration(), Policy()) + RealType bessel_quot = cyl_bessel_i(1, dist.concentration(), Policy()) / cyl_bessel_i(0, dist.concentration(), Policy()); return std::sqrt(-2 * std::log(bessel_quot)); } @@ -552,7 +507,7 @@ inline RealType median(const von_mises_distribution& dist) namespace detail { // float version of variance_impl template -inline RealType variance_impl(const von_mises_distribution& dist, +inline RealType variance_impl(const von_mises_distribution& dist, const boost::integral_constant&) { RealType conc = dist.concentration(); @@ -565,7 +520,7 @@ inline RealType variance_impl(const von_mises_distribution& di } else if(conc < 50) { - // Polynomial coefficients from + // 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 @@ -583,11 +538,11 @@ inline RealType variance_impl(const von_mises_distribution& di +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 @@ -598,7 +553,7 @@ inline RealType variance_impl(const von_mises_distribution& di 7.569820327857441460e-02f, 5.573513685531774810e-02f, 1.918908150536447035e-01f - }; + }; static const float P_denom[] = { 3.98942280401432677e-01f, 4.98677850501790847e-02f, @@ -609,14 +564,14 @@ inline RealType variance_impl(const von_mises_distribution& di 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, +inline RealType variance_impl(const von_mises_distribution& dist, const boost::integral_constant&) { RealType conc = dist.concentration(); @@ -629,7 +584,7 @@ inline RealType variance_impl(const von_mises_distribution& di } else if(conc < 500) { - // Polynomial coefficients from + // 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 @@ -684,7 +639,7 @@ inline RealType variance_impl(const von_mises_distribution& di 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 @@ -706,82 +661,74 @@ inline RealType variance_impl(const von_mises_distribution& di 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(RealType conc, const boost::integral_constant&) -//~ { - //~ BOOST_MATH_STD_USING - //~ if (conc < 15) - //~ { - //~ 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 (x < 50) - //~ { - //~ // Max error in interpolated form: 1.035e-21 - //~ // Max Error found at float80 precision = Poly: 1.885872e-21 - //~ static const RealType Y = 4.011702537536621093750e-01f; - //~ static const RealType P[] = { - //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, -2.227973351806078464328e-03), - //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.986778486088017419036e-02), - //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, 2.805066823812285310011e-02), - //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, 2.921443721160964964623e-02), - //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.517504941996594744052e-02), - //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, 6.316922639868793684401e-02), - //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.535891099168810015433e+00), - //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, -4.706078229522448308087e+01), - //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.351015763079160914632e+03), - //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, -2.948809013999277355098e+04), - //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.967598958582595361757e+05), - //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, -6.346924657995383019558e+06), - //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, 5.998794574259956613472e+07), - //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, -4.016371355801690142095e+08), - //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.768791455631826490838e+09), - //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, -4.441995678177349895640e+09), - //~ BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.482292669974971387738e+09) - //~ }; - //~ RealType result = exp(conc * (cos(x - mean) - 1.0)); - //~ result /= (boost::math::tools::evaluate_polynomial(P, RealType(1.0 / conc)) + Y) / sqrt(conc); - //~ * boost::math::constants::two_pi(); - //~ return result; - //~ } - //~ 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; - //~ } -//~ } +//~ // 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 < 15) + { + 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) + { + static const RealType Y = 4.011702537536621093750e-01f; + static const RealType P[] = { + BOOST_MATH_BIG_CONSTANT(RealType, 64, -2.227973351806078464328e-03), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.986778486088017419036e-02), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 2.805066823812285310011e-02), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 2.921443721160964964623e-02), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.517504941996594744052e-02), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 6.316922639868793684401e-02), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.535891099168810015433e+00), + BOOST_MATH_BIG_CONSTANT(RealType, 64, -4.706078229522448308087e+01), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.351015763079160914632e+03), + BOOST_MATH_BIG_CONSTANT(RealType, 64, -2.948809013999277355098e+04), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.967598958582595361757e+05), + BOOST_MATH_BIG_CONSTANT(RealType, 64, -6.346924657995383019558e+06), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 5.998794574259956613472e+07), + BOOST_MATH_BIG_CONSTANT(RealType, 64, -4.016371355801690142095e+08), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.768791455631826490838e+09), + BOOST_MATH_BIG_CONSTANT(RealType, 64, -4.441995678177349895640e+09), + BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.482292669974971387738e+09) + }; + + return 0; + } + 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; + } +} } // namespace detail template @@ -789,12 +736,12 @@ inline RealType variance(const von_mises_distribution& dist) { typedef boost::integral_constant::digits == 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 : + std::numeric_limits::digits <= 113 ? 113 : -1 > tag_type; @@ -807,24 +754,12 @@ inline RealType skewness(const von_mises_distribution& /*dist* return 0; } -// template -// inline RealType kurtosis(const von_mises_distribution& /*dist*/) -// { -// return 3; -// } - -// template -// inline RealType kurtosis_excess(const von_mises_distribution& /*dist*/) -// { - // return 0; -// } - template inline RealType entropy(const von_mises_distribution & dist) { using std::log; RealType arg = constants::two_pi() * cyl_bessel_i(0, dist.concentration(), Policy()); - RealType bessel_quot = cyl_bessel_i(1, 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; } @@ -837,6 +772,6 @@ inline RealType entropy(const von_mises_distribution & dist) // keep compilers that support two-phase lookup happy. #include -#endif // BOOST_STATS_VON_MISES_HPP +#endif // BOOST_MATH_DISTRIBUTIONS_VON_MISES_HPP diff --git a/test/test_von_mises.cpp b/test/test_von_mises.cpp index cbcdaf39b3..7a9a8d4728 100644 --- a/test/test_von_mises.cpp +++ b/test/test_von_mises.cpp @@ -1,3 +1,4 @@ +// Copyright Philipp C. J. Münster 2020. // Copyright Paul A. Bristow 2010. // Copyright John Maddock 2007. @@ -12,8 +13,6 @@ // From MathWorld--A Wolfram Web Resource. // http://mathworld.wolfram.com/VonMisesDistribution.html -#include // include directory /libs/math/src/tr1/ is needed. - #ifdef _MSC_VER # pragma warning (disable: 4127) // conditional expression is constant // caused by using if(std::numeric_limits::has_infinity) @@ -27,8 +26,9 @@ #include #include - using boost::math::von_mises_distribution; + using boost::math::von_mises_distribution; #include +#include "math_unit_test.hpp" #include "test_out_of_range.hpp" #include @@ -37,7 +37,7 @@ using std::endl; using std::setprecision; #include - using std::numeric_limits; + using std::numeric_limits; template void check_von_mises(RealType mean, RealType conc, RealType x, RealType p, RealType q, RealType tol) @@ -73,9 +73,9 @@ void check_von_mises(RealType mean, RealType conc, RealType x, RealType p, RealT template void test_spots(RealType) { - // Basic sanity checks + // Basic sanity checks // Check some bad parameters to the distribution, -#ifndef BOOST_NO_EXCEPTIONS +#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 @@ -96,7 +96,7 @@ void test_spots(RealType) 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 +#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 @@ -116,7 +116,7 @@ void test_spots(RealType) // // 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 + 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)) @@ -133,24 +133,33 @@ void test_spots(RealType) pdf(von_mises_distribution(3, 5), static_cast(3)), static_cast(0.86713652854235200257351846969777045343907L), // e^5/(2*pi*I0(5)) tolerance); - // Tests with larger numbers whose bessel function value exceeds the representable range 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, 125), static_cast(3)), - static_cast(4.45583423571762164664231477348923797078210L), // e^125/(2*pi*I0(125)) + 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, 1250), static_cast(3)), - static_cast(14.1033286206525302447068654967513214217021L), // e^1250/(2*pi*I0(1250)) + 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 = 1e-3f; // 1e-5 (as %) - // Some tests only pass at 1e-5 because values generated by + tolerance = 2e-3f; // 2e-5 (as %) - cout << "Tolerance for type " << typeid(RealType).name() << " is " << tolerance << " %" << endl; + cout << "Tolerance for type " << typeid(RealType).name() + << " is " << tolerance << " %" << endl; // test CDF for mean and interval edges check_von_mises( @@ -160,7 +169,7 @@ void test_spots(RealType) static_cast(0.5L), static_cast(0.5L), tolerance); - + check_von_mises( static_cast(0), static_cast(1), @@ -184,8 +193,8 @@ void test_spots(RealType) static_cast(1), static_cast(0.794355307434683479987678129735260058645499262629455722769L), static_cast(0.205644692565316520012321870264739941354500737370544277230L), - tolerance); - + tolerance); + check_von_mises( static_cast(0), static_cast(1), @@ -200,8 +209,8 @@ void test_spots(RealType) static_cast(1), static_cast(0.98096204546814689054581384796251763480020360394758184271L), static_cast(0.01903795453185310945418615203748236519979639605241815729L), - tolerance); - + tolerance); + check_von_mises( static_cast(0), static_cast(5), @@ -209,7 +218,7 @@ void test_spots(RealType) static_cast(0.01903795453185310945418615203748236519979639605241815729L), static_cast(0.98096204546814689054581384796251763480020360394758184271L), tolerance); - + // test CDF for high concentrations //~ check_von_mises( //~ static_cast(0), @@ -217,8 +226,8 @@ void test_spots(RealType) //~ static_cast(1), //~ static_cast(0.999999062404464440452233504489299782776166264467210572765L), //~ static_cast(9.37595535559547766495510700217223833735532789427234e-7L), - //~ tolerance); - + //~ tolerance); + //~ check_von_mises( //~ static_cast(0), //~ static_cast(25), @@ -226,15 +235,15 @@ void test_spots(RealType) //~ 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); - + //~ tolerance); + //~ check_von_mises( //~ static_cast(0), //~ static_cast(125), @@ -242,7 +251,7 @@ void test_spots(RealType) //~ static_cast(3.3546365686506670890489976669185096109e-20L), //~ static_cast(0.99999999999999999996645363431349332910951002333081490389L), //~ tolerance); - + RealType tol2 = boost::math::tools::epsilon() * 500; von_mises_distribution dist(2, 3); @@ -254,7 +263,7 @@ void test_spots(RealType) BOOST_CHECK_CLOSE( mean(dist) , static_cast(2), tol2); - // variance: + // variance: BOOST_CHECK_CLOSE( variance(von_mises_distribution(2, 0)) , static_cast(1), tol2); @@ -270,7 +279,7 @@ void test_spots(RealType) //~ BOOST_CHECK_CLOSE( //~ variance(von_mises_distribution(2, 125)) //~ , static_cast(0.004008064813593637057963834031301840442156903531539609019L), tol2); - + // std deviation: BOOST_CHECK_CLOSE( standard_deviation(dist) @@ -300,14 +309,6 @@ void test_spots(RealType) BOOST_CHECK_CLOSE( skewness(dist) , static_cast(0), tol2); - // kurtosis: - // BOOST_CHECK_CLOSE( - // kurtosis(dist) - // , static_cast(3), tol2); - // kurtosis excess: - // BOOST_CHECK_CLOSE( - // kurtosis_excess(dist) - // , static_cast(0), tol2); BOOST_CHECK_CLOSE( entropy(dist) @@ -330,11 +331,29 @@ void test_spots(RealType) // 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: @@ -351,8 +370,13 @@ BOOST_AUTO_TEST_CASE( test_main ) // (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 % + // Check symmetry of PDF and CDF + test_symmetry(0.0F); + test_symmetry(0.0); + #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS - //test_spots(0.0L); // Test long double. + test_spots(0.0L); // Test long double. + //test_symmetry(0.0L); #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x0582)) //test_spots(boost::math::concepts::real_concept(0.)); // Test real concept. #endif @@ -363,7 +387,7 @@ BOOST_AUTO_TEST_CASE( test_main ) "to pass." << std::endl; #endif - + } // BOOST_AUTO_TEST_CASE( test_main ) /* From 8230a6da72b2f3c8d1bb64ba811f83a4190157b1 Mon Sep 17 00:00:00 2001 From: pcjmuenster Date: Sat, 18 Jul 2020 20:05:59 +0200 Subject: [PATCH 7/8] [von Mises] add quad precision test --- .../boost/math/distributions/von_mises.hpp | 141 ++++++++++-------- test/Jamfile.v2 | 4 +- test/test_von_mises.cpp | 38 ++--- 3 files changed, 99 insertions(+), 84 deletions(-) diff --git a/include/boost/math/distributions/von_mises.hpp b/include/boost/math/distributions/von_mises.hpp index da218eada4..00c928543f 100644 --- a/include/boost/math/distributions/von_mises.hpp +++ b/include/boost/math/distributions/von_mises.hpp @@ -118,9 +118,9 @@ inline RealType pdf_impl(const von_mises_distribution& dist, R 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 + else // exp(88) > MAX_FLOAT { - // we make use of I0(conc) = exp(conc) * P(conc), with polynomial P + // 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, @@ -177,40 +177,11 @@ inline RealType pdf_impl(const von_mises_distribution& dist, R RealType conc = dist.concentration(); BOOST_MATH_STD_USING - if (conc < 15) + 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 if (conc < 50) - { - // Max error in interpolated form: 1.035e-21 - // Max Error found at float80 precision = Poly: 1.885872e-21 - static const RealType Y = 4.011702537536621093750e-01f; - static const RealType P[] = { - BOOST_MATH_BIG_CONSTANT(RealType, 64, -2.227973351806078464328e-03), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.986778486088017419036e-02), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 2.805066823812285310011e-02), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 2.921443721160964964623e-02), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.517504941996594744052e-02), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 6.316922639868793684401e-02), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.535891099168810015433e+00), - BOOST_MATH_BIG_CONSTANT(RealType, 64, -4.706078229522448308087e+01), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.351015763079160914632e+03), - BOOST_MATH_BIG_CONSTANT(RealType, 64, -2.948809013999277355098e+04), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.967598958582595361757e+05), - BOOST_MATH_BIG_CONSTANT(RealType, 64, -6.346924657995383019558e+06), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 5.998794574259956613472e+07), - BOOST_MATH_BIG_CONSTANT(RealType, 64, -4.016371355801690142095e+08), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.768791455631826490838e+09), - BOOST_MATH_BIG_CONSTANT(RealType, 64, -4.441995678177349895640e+09), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.482292669974971387738e+09) - }; - RealType result = exp(conc * (cos(x - mean) - 1.0)); - result /= (boost::math::tools::evaluate_polynomial(P, RealType(1.0 / conc)) + Y) / sqrt(conc) - * boost::math::constants::two_pi(); - return result; - } else { // Bessel I0 over[50, INF] @@ -242,6 +213,16 @@ inline RealType pdf_impl(const von_mises_distribution& dist, R 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 @@ -487,9 +468,10 @@ inline RealType mean(const von_mises_distribution& dist) 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 std::sqrt(-2 * std::log(bessel_quot)); + return sqrt(-2 * log(bessel_quot)); } template @@ -666,44 +648,19 @@ inline RealType variance_impl(const von_mises_distribution& di } } -//~ // long double version of variance_impl +// 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 < 15) + 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 if (conc < 50) - { - static const RealType Y = 4.011702537536621093750e-01f; - static const RealType P[] = { - BOOST_MATH_BIG_CONSTANT(RealType, 64, -2.227973351806078464328e-03), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.986778486088017419036e-02), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 2.805066823812285310011e-02), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 2.921443721160964964623e-02), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.517504941996594744052e-02), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 6.316922639868793684401e-02), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.535891099168810015433e+00), - BOOST_MATH_BIG_CONSTANT(RealType, 64, -4.706078229522448308087e+01), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.351015763079160914632e+03), - BOOST_MATH_BIG_CONSTANT(RealType, 64, -2.948809013999277355098e+04), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.967598958582595361757e+05), - BOOST_MATH_BIG_CONSTANT(RealType, 64, -6.346924657995383019558e+06), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 5.998794574259956613472e+07), - BOOST_MATH_BIG_CONSTANT(RealType, 64, -4.016371355801690142095e+08), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 1.768791455631826490838e+09), - BOOST_MATH_BIG_CONSTANT(RealType, 64, -4.441995678177349895640e+09), - BOOST_MATH_BIG_CONSTANT(RealType, 64, 4.482292669974971387738e+09) - }; - - return 0; - } else { static const RealType P[] = { @@ -729,6 +686,68 @@ inline RealType variance_impl(const von_mises_distribution& di 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 @@ -757,7 +776,7 @@ inline RealType skewness(const von_mises_distribution& /*dist* template inline RealType entropy(const von_mises_distribution & dist) { - using std::log; + 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()); diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index db6c4bc282..e92723b8a3 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -793,7 +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 ] + [ 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 index 7a9a8d4728..3a89d17650 100644 --- a/test/test_von_mises.cpp +++ b/test/test_von_mises.cpp @@ -19,8 +19,11 @@ // and if (std::numeric_limits::has_quiet_NaN) #endif -#include +#include +#include + using boost::multiprecision::float128; #include // for real_concept + #define BOOST_TEST_MAIN #include // Boost.Test #include @@ -28,7 +31,9 @@ #include using boost::math::von_mises_distribution; #include + #include "math_unit_test.hpp" +#include "pch.hpp" #include "test_out_of_range.hpp" #include @@ -370,34 +375,23 @@ BOOST_AUTO_TEST_CASE( test_main ) // (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. + test_spots(static_cast(0)); + // Check symmetry of PDF and CDF test_symmetry(0.0F); test_symmetry(0.0); - -#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS - test_spots(0.0L); // Test long double. - //test_symmetry(0.0L); -#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x0582)) - //test_spots(boost::math::concepts::real_concept(0.)); // Test real concept. -#endif -#else - std::cout << "The long double tests have been disabled on this platform " - "either because the long double overloads of the usual math functions are " - "not available at all, or because they are too inaccurate for these tests " - "to pass." << std::endl; -#endif - - + test_symmetry(0.0L); + test_symmetry(static_cast(0)); } // BOOST_AUTO_TEST_CASE( test_main ) /* - +./test_von_mises.exe Output: -Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_von_mises.exe" Running 1 test case... -Tolerance for type float is 0.01 % -Tolerance for type double is 0.01 % -Tolerance for type long double is 0.01 % -Tolerance for type class boost::math::concepts::real_concept is 0.01 % +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 3e4f2eb6de6ed975c12ba95be2a588c7cc88f5a5 Mon Sep 17 00:00:00 2001 From: pcjmuenster Date: Sun, 23 Aug 2020 18:21:01 +0200 Subject: [PATCH 8/8] [von_mises] add ULP plots for CDF and PDF [CI SKIP] --- .../von_mises/von_mises_ulps_cdf_4d.svg | 1048 ++++++++++++++++ .../von_mises/von_mises_ulps_cdf_4f.svg | 1049 ++++++++++++++++ .../von_mises/von_mises_ulps_cdf_64d.svg | 1070 +++++++++++++++++ .../von_mises/von_mises_ulps_cdf_64f.svg | 1066 ++++++++++++++++ .../von_mises/von_mises_ulps_pdf_4d.svg | 1050 ++++++++++++++++ .../von_mises/von_mises_ulps_pdf_4f.svg | 1050 ++++++++++++++++ .../von_mises/von_mises_ulps_pdf_64d.svg | 1050 ++++++++++++++++ .../von_mises/von_mises_ulps_pdf_64f.svg | 1050 ++++++++++++++++ reporting/accuracy/von_mises_ulps.cpp | 81 ++ 9 files changed, 8514 insertions(+) create mode 100644 doc/graphs/von_mises/von_mises_ulps_cdf_4d.svg create mode 100644 doc/graphs/von_mises/von_mises_ulps_cdf_4f.svg create mode 100644 doc/graphs/von_mises/von_mises_ulps_cdf_64d.svg create mode 100644 doc/graphs/von_mises/von_mises_ulps_cdf_64f.svg create mode 100644 doc/graphs/von_mises/von_mises_ulps_pdf_4d.svg create mode 100644 doc/graphs/von_mises/von_mises_ulps_pdf_4f.svg create mode 100644 doc/graphs/von_mises/von_mises_ulps_pdf_64d.svg create mode 100644 doc/graphs/von_mises/von_mises_ulps_pdf_64f.svg create mode 100644 reporting/accuracy/von_mises_ulps.cpp 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 @@ + + + + + + + +-4.393e+004 + +-3.753e+004 + +-3.112e+004 + +-2.471e+004 + +-1.83e+004 + +-1.19e+004 + +-5490 + +916.6 + +-2.513 + +-1.885 + +-1.257 + +-0.6283 + +0 + +0.6283 + +1.257 + +1.885 + +2.513 + +3.142 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 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 @@ + + + + + + + +1.111e+004 + +2.223e+004 + +3.335e+004 + +4.447e+004 + +5.559e+004 + +6.672e+004 + +7.784e+004 + +8.896e+004 + +-2.513 + +-1.885 + +-1.257 + +-0.6283 + +0 + +0.6283 + +1.257 + +1.885 + +2.513 + +3.142 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 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 @@ + + + + + + + +2.682e+306 + +5.365e+306 + +8.047e+306 + +1.073e+307 + +1.341e+307 + +1.609e+307 + +1.878e+307 + +2.146e+307 + +-2.513 + +-1.885 + +-1.257 + +-0.6283 + +0 + +0.6283 + +1.257 + +1.885 + +2.513 + +3.142 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 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 @@ + + + + + + + +-1.362e+007 + +-1.056e+007 + +-7.489e+006 + +-4.423e+006 + +-1.356e+006 + +1.71e+006 + +4.777e+006 + +7.843e+006 + +-2.513 + +-1.885 + +-1.257 + +-0.6283 + +0 + +0.6283 + +1.257 + +1.885 + +2.513 + +3.142 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 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 @@ + + + + + + + +-1.84 + +-1.138 + +-0.4357 + +0.2666 + +0.9689 + +1.671 + +2.373 + +3.076 + +0.3142 + +0.6283 + +0.9425 + +1.257 + +1.571 + +1.885 + +2.199 + +2.513 + +2.827 + +3.142 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 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 @@ + + + + + + + +-2.422 + +-1.721 + +-1.02 + +-0.3184 + +0.3827 + +1.084 + +1.785 + +2.486 + +0.3142 + +0.6283 + +0.9425 + +1.257 + +1.571 + +1.885 + +2.199 + +2.513 + +2.827 + +3.142 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 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 @@ + + + + + + + +-23.42 + +-16.12 + +-8.814 + +-1.513 + +5.789 + +13.09 + +20.39 + +27.69 + +0.3142 + +0.6283 + +0.9425 + +1.257 + +1.571 + +1.885 + +2.199 + +2.513 + +2.827 + +3.142 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 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 @@ + + + + + + + +-21.17 + +-14.22 + +-7.281 + +-0.3367 + +6.607 + +13.55 + +20.5 + +27.44 + +0.3142 + +0.6283 + +0.9425 + +1.257 + +1.571 + +1.885 + +2.199 + +2.513 + +2.827 + +3.142 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 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; +}