diff --git a/doc/math.qbk b/doc/math.qbk index 385c93a5e8..d56fef2626 100644 --- a/doc/math.qbk +++ b/doc/math.qbk @@ -339,6 +339,10 @@ and use the function's name as the link text.] [def __hermite [link math_toolkit.sf_poly.hermite hermite]] [def __cardinal_b_splines [link math_toolkit.sf_poly.cardinal_b_splines cardinal_b_splines]] +[/logistic functions] +[def __logit [link math_toolkit.logistic.logit logit]] +[def __logistic_sigmoid [link math_toolkit.logistic.logistic_sigmoid logistic_sigmoid]] + [/Misc] [def __expint [link math_toolkit.expint.expint_i expint]] [def __spherical_harmonic [link math_toolkit.sf_poly.sph_harm spherical_harmonic]] @@ -666,6 +670,11 @@ and as a CD ISBN 0-9504833-2-X 978-0-9504833-2-0, Classification 519.2-dc22. [include sf/jacobi.qbk] [endsect] [/section:sf_poly Polynomials] +[section:logistic Logistic Functions] +[include sf/logit.qbk] +[include sf/logistic_sigmoid.qbk] +[endsect] [/section:logistic Logistic Functions] + [section:bessel Bessel Functions] [include sf/bessel_introduction.qbk] [include sf/bessel_jy.qbk] diff --git a/doc/sf/logistic_sigmoid.qbk b/doc/sf/logistic_sigmoid.qbk new file mode 100644 index 0000000000..bf1c6c37fd --- /dev/null +++ b/doc/sf/logistic_sigmoid.qbk @@ -0,0 +1,30 @@ +[/ + Copyright Matt Borland 2025 + Distributed under the Boost Software License, Version 1.0. + (See accompanying file LICENSE_1_0.txt or copy at + http://www.boost.org/LICENSE_1_0.txt). +] + +[section:logistic_sigmoid logistic_sigmoid] + +[h4 Synopsis] + + #include + + namespace boost { namespace math { + + template + RealType logistic_sigmoid(RealType x, const Policy&); + + template + RealType logistic_sigmoid(RealType x) + + }} // namespaces + +[h4 Description] + +Returns the [@https://en.wikipedia.org/wiki/Logistic_function logistic sigmoid function] +This function is broadly useful, and is used to compute the CDF of the logistic distribution. +It is also sometimes referred to as expit as it is the inverse of the logit function. + +[endsect] diff --git a/doc/sf/logit.qbk b/doc/sf/logit.qbk new file mode 100644 index 0000000000..7daebe667f --- /dev/null +++ b/doc/sf/logit.qbk @@ -0,0 +1,30 @@ +[/ + Copyright Matt Borland 2025 + Distributed under the Boost Software License, Version 1.0. + (See accompanying file LICENSE_1_0.txt or copy at + http://www.boost.org/LICENSE_1_0.txt). +] + +[section:logit logit] + +[h4 Synopsis] + + #include + + namespace boost { namespace math { + + template + RealType logit(RealType x, const Policy&); + + template + RealType logit(RealType x) + + }} // namespaces + +[h4 Description] + +Returns the [@https://en.wikipedia.org/wiki/Logit logit function] +This function is broadly useful, and is used to compute the Quantile of the logistic distribution. +The inverse of this function is the logistic_sigmoid. + +[endsect] diff --git a/include/boost/math/distributions/logistic.hpp b/include/boost/math/distributions/logistic.hpp index 69ec818abe..bc2dd1ff5a 100644 --- a/include/boost/math/distributions/logistic.hpp +++ b/include/boost/math/distributions/logistic.hpp @@ -19,6 +19,8 @@ #include #include #include +#include +#include namespace boost { namespace math { @@ -144,13 +146,10 @@ namespace boost { namespace math { { return result; } - BOOST_MATH_STD_USING - RealType power = (location - x) / scale; - if(power > tools::log_max_value()) - return 0; - if(power < -tools::log_max_value()) - return 1; - return 1 / (1 + exp(power)); + + using promoted_real_type = typename policies::evaluation::type; + promoted_real_type power = (static_cast(x) - static_cast(location)) / static_cast(scale); + return logistic_sigmoid(power, policies::make_forwarding_policy_t()); } template @@ -199,7 +198,6 @@ namespace boost { namespace math { template BOOST_MATH_GPU_ENABLED inline RealType quantile(const logistic_distribution& dist, const RealType& p) { - BOOST_MATH_STD_USING RealType location = dist.location(); RealType scale = dist.scale(); @@ -221,21 +219,13 @@ namespace boost { namespace math { { return policies::raise_overflow_error(function,"probability argument is 1, must be >0 and <1",Policy()); } - //Expressions to try - //return location+scale*log(p/(1-p)); - //return location+scale*log1p((2*p-1)/(1-p)); - - //return location - scale*log( (1-p)/p); - //return location - scale*log1p((1-2*p)/p); - //return -scale*log(1/p-1) + location; - return location - scale * log((1 - p) / p); + return location + scale * logit(p, Policy()); } // RealType quantile(const logistic_distribution& dist, const RealType& p) template BOOST_MATH_GPU_ENABLED inline RealType cdf(const complemented2_type, RealType>& c) { - BOOST_MATH_STD_USING RealType location = c.dist.location(); RealType scale = c.dist.scale(); RealType x = c.param; @@ -259,12 +249,10 @@ namespace boost { namespace math { { return result; } - RealType power = (x - location) / scale; - if(power > tools::log_max_value()) - return 0; - if(power < -tools::log_max_value()) - return 1; - return 1 / (1 + exp(power)); + + using promoted_real_type = typename policies::evaluation::type; + promoted_real_type power = (static_cast(location) - static_cast(x)) / static_cast(scale); + return logistic_sigmoid(power, policies::make_forwarding_policy_t()); } template @@ -306,7 +294,6 @@ namespace boost { namespace math { template BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type, RealType>& c) { - BOOST_MATH_STD_USING RealType scale = c.dist.scale(); RealType location = c.dist.location(); constexpr auto function = "boost::math::quantile(const complement(logistic_distribution<%1%>&), %1%)"; @@ -328,15 +315,8 @@ namespace boost { namespace math { { return policies::raise_overflow_error(function,"probability argument is 0, but must be >0 and <1",Policy()); } - //Expressions to try - //return location+scale*log((1-q)/q); - return location + scale * log((1 - q) / q); - - //return location-scale*log(q/(1-q)); - //return location-scale*log1p((2*q-1)/(1-q)); - //return location+scale*log(1/q-1); - //return location+scale*log1p(1/q-2); + return location - scale * logit(q, Policy()); } template diff --git a/include/boost/math/policies/policy.hpp b/include/boost/math/policies/policy.hpp index bec8886408..5fb41bcd1a 100644 --- a/include/boost/math/policies/policy.hpp +++ b/include/boost/math/policies/policy.hpp @@ -1003,6 +1003,25 @@ struct is_noexcept_error_policy && (t8::value != throw_on_error) && (t8::value != user_error)); }; +// Generate a forwarding policy to stop further promotion from occurring +// For example if a special function for float promotes to double, we don't want the next +// function in the call chain to then promote to long double +template +struct make_forwarding_policy +{ + using type = + typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> + >::type; +}; + +template +using make_forwarding_policy_t = typename make_forwarding_policy::type; + }}} // namespaces #endif // BOOST_MATH_POLICY_HPP diff --git a/include/boost/math/special_functions/logistic_sigmoid.hpp b/include/boost/math/special_functions/logistic_sigmoid.hpp new file mode 100644 index 0000000000..925b576ca3 --- /dev/null +++ b/include/boost/math/special_functions/logistic_sigmoid.hpp @@ -0,0 +1,46 @@ +// Copyright Matt Borland 2025. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_SF_EXPIT_HPP +#define BOOST_MATH_SF_EXPIT_HPP + +#include +#include +#include + +namespace boost { +namespace math { + +template +RealType logistic_sigmoid(RealType x, const Policy&) +{ + BOOST_MATH_STD_USING + + using promoted_real_type = typename policies::evaluation::type; + + if(-x >= tools::log_max_value()) + { + return static_cast(0); + } + if(-x <= -tools::log_max_value()) + { + return static_cast(1); + } + + const auto res {static_cast(1 / (1 + exp(static_cast(-x))))}; + return res; +} + +template +RealType logistic_sigmoid(RealType x) +{ + return logistic_sigmoid(x, policies::policy<>()); +} + +} // namespace math +} // namespace boost + +#endif // BOOST_MATH_SF_EXPIT_HPP diff --git a/include/boost/math/special_functions/logit.hpp b/include/boost/math/special_functions/logit.hpp new file mode 100644 index 0000000000..8a2ddf68dd --- /dev/null +++ b/include/boost/math/special_functions/logit.hpp @@ -0,0 +1,56 @@ +// Copyright Matt Borland 2025. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_SF_LOGIT_HPP +#define BOOST_MATH_SF_LOGIT_HPP + +#include +#include +#include +#include +#include + +namespace boost { +namespace math { + +template +RealType logit(RealType p, const Policy&) +{ + BOOST_MATH_STD_USING + using std::atanh; + + using promoted_real_type = typename policies::evaluation::type; + + if (p < tools::min_value()) + { + return -policies::raise_overflow_error("logit", "sub-normals will overflow ln(x/(1-x))", Policy()); + } + + static const RealType crossover {RealType{1}/4}; + const auto promoted_p {static_cast(p)}; + RealType result {}; + if (p > crossover) + { + result = static_cast(2 * atanh(2 * promoted_p - 1)); + } + else + { + result = static_cast(log(promoted_p / (1 - promoted_p))); + } + + return result; +} + +template +RealType logit(RealType p) +{ + return logit(p, policies::policy<>()); +} + +} // namespace math +} // namespace boost + +#endif // BOOST_MATH_SF_LOGIT_HPP diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 7a7be3f5d9..3843264429 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -590,6 +590,8 @@ test-suite special_fun : [ run test_sinc.cpp /boost/test//boost_unit_test_framework pch_light ] [ run test_fibonacci.cpp /boost/test//boost_unit_test_framework ] [ run test_prime.cpp /boost/test//boost_unit_test_framework ] + [ run test_logistic_sigmoid.cpp ] + [ run test_logit.cpp ] ; test-suite distribution_tests : diff --git a/test/git_issue_1294.cpp b/test/git_issue_1294.cpp new file mode 100644 index 0000000000..047a8ea7b3 --- /dev/null +++ b/test/git_issue_1294.cpp @@ -0,0 +1,33 @@ +// Copyright 2025 Matt Borland +// 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) +// +// See: https://github.com/boostorg/math/issues/1294 + +#include +#include "math_unit_test.hpp" + +int main() +{ + using namespace boost::math::policies; + using boost::math::logistic_distribution; + + typedef policy< + promote_float, + promote_double + > with_promotion; + + constexpr double p = 2049.0/4096; + constexpr double ref = 9.76562577610225755e-04; + + logistic_distribution dist_promote; + const double x = quantile(dist_promote, p); + + // Previously we had: 9.76562577610170027e-04 + // Which is an ULP distance of 256 + CHECK_ULP_CLOSE(x, ref, 1); + + return boost::math::test::report_errors(); +} diff --git a/test/math_unit_test.hpp b/test/math_unit_test.hpp index 648f7de9c5..bd85caa86b 100644 --- a/test/math_unit_test.hpp +++ b/test/math_unit_test.hpp @@ -58,7 +58,7 @@ bool check_mollified_close(Real expected, Real computed, Real tol, std::string c } using std::max; using std::abs; - Real denom = (max)(abs(expected), Real(1)); + Real denom = (max)(Real(abs(expected)), Real(1)); Real mollified_relative_error = abs(expected - computed)/denom; if (mollified_relative_error > tol) { diff --git a/test/test_logistic_sigmoid.cpp b/test/test_logistic_sigmoid.cpp new file mode 100644 index 0000000000..c0c8a04920 --- /dev/null +++ b/test/test_logistic_sigmoid.cpp @@ -0,0 +1,91 @@ +// Copyright Matt Borland 2025. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include +#include "math_unit_test.hpp" +#include +#include +#include + +#pragma STDC FENV_ACCESS ON + +template +void test() +{ + const std::array x_values = { + 0, + 1, + 1000, + 0.5, + 0.75 + }; + const std::array y_values = { + static_cast(1) / 2, + static_cast(0.73105857863000487925115924182183627436514464016505651927636590791904045307), + static_cast(1), + static_cast(0.62245933120185456463890056574550847875327936530891016305943716265854500), + static_cast(0.6791786991753929731596801157765790212342212482195760219829517436805) + }; + + for (std::size_t i = 0; i < x_values.size(); ++i) + { + const RealType test_value {boost::math::logistic_sigmoid(x_values[i])}; + BOOST_MATH_IF_CONSTEXPR (std::is_same::value || std::is_same::value) + { + CHECK_ULP_CLOSE(test_value, y_values[i], 1); + } + else + { + RealType comparison_value = y_values[i]; + CHECK_MOLLIFIED_CLOSE(test_value, comparison_value, static_cast(1e-15)); + } + + bool fe {false}; + if (std::fetestexcept(FE_OVERFLOW)) + { + fe = true; // LCOV_EXCL_LINE + std::cerr << "FE_OVERFLOW" << std::endl; // LCOV_EXCL_LINE + } + if (std::fetestexcept(FE_UNDERFLOW)) + { + fe = true; // LCOV_EXCL_LINE + std::cerr << "FE_UNDERFLOW" << std::endl; // LCOV_EXCL_LINE + } + if (std::fetestexcept(FE_DIVBYZERO)) + { + fe = true; // LCOV_EXCL_LINE + std::cerr << "FE_DIVBYZERO" << std::endl; // LCOV_EXCL_LINE + } + if (std::fetestexcept(FE_INVALID)) + { + fe = true; // LCOV_EXCL_LINE + std::cerr << "FE_INVALID" << std::endl; // LCOV_EXCL_LINE + } + + CHECK_EQUAL(fe, false); + } +} + +int main() +{ + std::feclearexcept(FE_ALL_EXCEPT); + test(); + + std::feclearexcept(FE_ALL_EXCEPT); + test(); + + std::feclearexcept(FE_ALL_EXCEPT); + test(); + + std::feclearexcept(FE_ALL_EXCEPT); + test(); + + test(); + + return boost::math::test::report_errors(); +} diff --git a/test/test_logit.cpp b/test/test_logit.cpp new file mode 100644 index 0000000000..bef357e87a --- /dev/null +++ b/test/test_logit.cpp @@ -0,0 +1,107 @@ +// Copyright Matt Borland 2025. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt +// or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include +#include "math_unit_test.hpp" +#include +#include +#include + +#pragma STDC FENV_ACCESS ON + +template +void test() +{ + const std::array x_values = { + 0.01, + 0.24, + 0.5, + 0.75, + 0.995, + }; + const std::array y_values = { + RealType{-4.595119850134589926852434051810180709116687969582916078687956376405}, + RealType{-1.15267950993838545919655007350715126451438856911612411268258589327840479}, + RealType{0}, + RealType{1.09861228866810969139524523692252570464749055782274945173469433363749429}, + RealType{5.2933048247244923954101212918685372018911052805694724989064609879440992} + }; + + for (std::size_t i = 0; i < x_values.size(); ++i) + { + const RealType test_value {boost::math::logit(x_values[i])}; + + BOOST_MATH_IF_CONSTEXPR (std::is_same::value || std::is_same::value) + { + CHECK_ULP_CLOSE(test_value, y_values[i], 5); + } + else + { + CHECK_MOLLIFIED_CLOSE(test_value, y_values[i], 1e-15); + } + + bool fe {false}; + if (std::fetestexcept(FE_OVERFLOW)) + { + fe = true; // LCOV_EXCL_LINE + std::cerr << "FE_OVERFLOW" << std::endl; // LCOV_EXCL_LINE + } + if (std::fetestexcept(FE_UNDERFLOW)) + { + fe = true; // LCOV_EXCL_LINE + std::cerr << "FE_UNDERFLOW" << std::endl; // LCOV_EXCL_LINE + } + if (std::fetestexcept(FE_DIVBYZERO)) + { + fe = true; // LCOV_EXCL_LINE + std::cerr << "FE_DIVBYZERO" << std::endl; // LCOV_EXCL_LINE + } + if (std::fetestexcept(FE_INVALID)) + { + fe = true; // LCOV_EXCL_LINE + std::cerr << "FE_INVALID" << std::endl; // LCOV_EXCL_LINE + } + + CHECK_EQUAL(fe, false); + } + + #if defined(_CPPUNWIND) || defined(__EXCEPTIONS) + + BOOST_MATH_IF_CONSTEXPR (std::is_arithmetic::value) + { + bool thrown {false}; + try + { + boost::math::logit(std::numeric_limits::denorm_min()); + } + catch (...) + { + thrown = true; + } + + CHECK_EQUAL(thrown, true); + } + + #endif // Exceptional environments +} + +int main() +{ + std::feclearexcept(FE_ALL_EXCEPT); + test(); + + std::feclearexcept(FE_ALL_EXCEPT); + test(); + + std::feclearexcept(FE_ALL_EXCEPT); + test(); + + std::feclearexcept(FE_ALL_EXCEPT); + test(); + + return boost::math::test::report_errors(); +}