From de6a562c6d5bea8bcb2f8c10d41d312ae4690d80 Mon Sep 17 00:00:00 2001 From: dschmitz89 Date: Sun, 12 Apr 2026 15:10:26 +0200 Subject: [PATCH 01/18] ENH: add parameter finder for degrees or freedom --- doc/distributions/students_t.qbk | 13 ++ .../boost/math/distributions/students_t.hpp | 124 ++++++++++++++++++ test/test_students_t.cpp | 78 +++++++++++ 3 files changed, 215 insertions(+) diff --git a/doc/distributions/students_t.qbk b/doc/distributions/students_t.qbk index 3396048c5f..0a95080bdb 100644 --- a/doc/distributions/students_t.qbk +++ b/doc/distributions/students_t.qbk @@ -29,6 +29,11 @@ RealType beta, RealType sd, RealType hint = 100); + + // degrees of freedom inversion from a quantile and probability: + BOOST_MATH_GPU_ENABLED static RealType invert_probability_with_respect_to_degrees_of_freedom( + RealType x, + RealType p); }; }} // namespaces @@ -106,6 +111,13 @@ For more information on this function see the [@http://www.itl.nist.gov/div898/handbook/prc/section2/prc222.htm NIST Engineering Statistics Handbook]. + BOOST_MATH_GPU_ENABLED static RealType invert_probability_with_respect_to_degrees_of_freedom( + RealType x, + RealType p); + +Returns the degrees of freedom [nu] such that CDF(/x/; [nu]) = /p/. +Requires 0 < /p/ < 1 and /x/ [ne] 0, otherwise calls __domain_error. + [h4 Non-member Accessors] All the [link math_toolkit.dist_ref.nmp usual non-member accessor functions] that are generic to all @@ -164,6 +176,7 @@ without the subtraction implied above.]] [[skewness][if (v > 3) 0 else NaN ]] [[kurtosis][if (v > 4) 3 * (v - 2) \/ (v - 4) else NaN]] [[kurtosis excess][if (v > 4) 6 \/ (df - 4) else NaN]] +[[invert_probability_with_respect_to_degrees_of_freedom][Uses a 2nd-order Edgeworth expansion as initial guess for a numerical root finder]] ] If the moment index /k/ is less than /v/, then the moment is undefined. diff --git a/include/boost/math/distributions/students_t.hpp b/include/boost/math/distributions/students_t.hpp index 39f20d6e41..830a181cf5 100644 --- a/include/boost/math/distributions/students_t.hpp +++ b/include/boost/math/distributions/students_t.hpp @@ -58,6 +58,10 @@ class students_t_distribution RealType sd, RealType hint = 100); + BOOST_MATH_GPU_ENABLED static RealType invert_probability_with_respect_to_degrees_of_freedom( + RealType x, + RealType p); + private: // Data member: RealType df_; // degrees of freedom is a real number > 0 or +infinity. @@ -308,6 +312,23 @@ struct sample_size_func RealType alpha, beta, ratio; }; +template +struct cdf_to_df_func +{ + BOOST_MATH_GPU_ENABLED cdf_to_df_func(RealType x_val, RealType p_val, bool c) + : x(x_val), p(p_val), comp(c) {} + + BOOST_MATH_GPU_ENABLED RealType operator()(const RealType& df) + { + students_t_distribution t(df); + return comp ? + RealType(p - cdf(complement(t, x))) + : RealType(cdf(t, x) - p); + } + RealType x, p; + bool comp; +}; + } // namespace detail template @@ -344,6 +365,109 @@ BOOST_MATH_GPU_ENABLED RealType students_t_distribution::find_ return result; } +template +BOOST_MATH_GPU_ENABLED RealType students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + RealType x, + RealType p) +{ + BOOST_MATH_STD_USING // for ADL of std functions + constexpr auto function = "boost::math::students_t_distribution<%1%>::invert_probability_with_respect_to_degrees_of_freedom"; + // + // Check for domain errors: + // + RealType error_result; + if (false == detail::check_probability(function, p, &error_result, Policy())) + return error_result; + + // Use symmetry: CDF(-x; df) = 1 - CDF(x; df), so always work with x >= 0. + RealType x_abs = (x < 0) ? -x : x; + RealType p_adj = (x < 0) ? (1 - p) : p; + + if (x_abs == 0) + { + // CDF(0; df) = 0.5 for all df; only solvable when p == 0.5. + if (p_adj == static_cast(0.5)) + return boost::math::numeric_limits::infinity(); + return policies::raise_domain_error( + function, + "No degrees of freedom can satisfy CDF(0; df) == %1% (must be 0.5).", + p, Policy()); + } + + // + // Step 1: Edgeworth warm start. + // F(x; nu) ~ Phi(x) - phi(x)(x + x^3)/(4*nu) + phi(x)(3x + 5x^3 + 7x^5 - 3x^7)/(96*nu^2) + // Substituting u = 1/nu and setting F(x; nu) = p gives: + // B*u^2 - A*u + (Phi(x) - p) = 0 + // where: + // A = phi(x) * (x + x^3) / 4 + // B = phi(x) * (3x + 5x^3 + 7x^5 - 3x^7) / 96 + // + normal_distribution std_normal(0, 1); + RealType phi = pdf(std_normal, x_abs); + RealType Phi = cdf(std_normal, x_abs); + + RealType x2 = x_abs * x_abs; + RealType x3 = x2 * x_abs; + RealType x5 = x3 * x2; + RealType x7 = x5 * x2; + + RealType A = phi * (x_abs + x3) / 4; + RealType B = phi * (3 * x_abs + 5 * x3 + 7 * x5 - 3 * x7) / 96; + RealType C = Phi - p_adj; + + RealType hint = static_cast(0.01); + RealType discriminant = A * A - 4 * B * C; + if (discriminant >= 0 && B != 0) + { + RealType sqrt_disc = sqrt(discriminant); + // Two roots of B*u^2 - A*u + C = 0; pick the smallest positive u (largest nu = 1/u). + RealType u1 = (A - sqrt_disc) / (2 * B); + RealType u2 = (A + sqrt_disc) / (2 * B); + RealType u = -1; + if (u1 > 0 && u2 > 0) + u = (u1 < u2) ? u1 : u2; + else if (u1 > 0) + u = u1; + else if (u2 > 0) + u = u2; + + if (u > 0) + { + RealType nu_hat = 1 / u; + // Step 2: validate by checking relative residual of the exact CDF. + students_t_distribution t_hat(nu_hat); + RealType exact_cdf = cdf(t_hat, x_abs); + RealType residual = (exact_cdf > p_adj) ? (exact_cdf - p_adj) : (p_adj - exact_cdf); + RealType relative_residual = (p_adj != 0) ? residual / p_adj : residual; + if (relative_residual <= static_cast(0.1)) + hint = nu_hat; + } + } + + // + // Step 3: root-find with bracket_and_solve_root on f(df) = CDF(x; df) - p. + // For x > 0, CDF(x; df) is strictly increasing in df, so f is a rising function. + // Pass min(p_adj, q_adj) always and use the complement CDF when p_adj > q_adj + // to avoid catastrophic cancellation near 1 (mirrors non_centrality_finder_f). + // + RealType q_adj = 1 - p_adj; + detail::cdf_to_df_func f(x_abs, p_adj < q_adj ? p_adj : q_adj, p_adj >= q_adj); + tools::eps_tolerance tol(policies::digits()); + boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); + boost::math::pair r = tools::bracket_and_solve_root( + f, hint, RealType(2), true, tol, max_iter, Policy()); + RealType result = r.first + (r.second - r.first) / 2; + if (max_iter >= policies::get_max_root_iterations()) + { + return policies::raise_evaluation_error(function, + "Unable to locate solution in a reasonable time: either there is no answer to " + "the degrees of freedom or the answer is infinite. Current best guess is %1%", + result, Policy()); + } + return result; +} + template BOOST_MATH_GPU_ENABLED inline RealType mode(const students_t_distribution& /*dist*/) { diff --git a/test/test_students_t.cpp b/test/test_students_t.cpp index a08cf75a01..fa113ca602 100644 --- a/test/test_students_t.cpp +++ b/test/test_students_t.cpp @@ -548,6 +548,84 @@ void test_spots(RealType) static_cast(1.0))), 9); + // Tests for invert_probability_with_respect_to_degrees_of_freedom. + // Each case is derived from the CDF spot tests above: the exact df is known, + // so we verify that inverting CDF(x; df) = p recovers df to tight tolerance. + { + RealType tol_inv = static_cast(1000) * boost::math::tools::epsilon(); + + // df=2, x=-6.96455673428326, p=0.01 + BOOST_CHECK_CLOSE_FRACTION( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(-6.96455673428326), + static_cast(0.01)), + static_cast(2), + tol_inv); + + // df=5, x=-3.36492999890721, p=0.01 + BOOST_CHECK_CLOSE_FRACTION( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(-3.36492999890721), + static_cast(0.01)), + static_cast(5), + tol_inv); + + // df=5, x=-0.559429644, p=0.3 + BOOST_CHECK_CLOSE_FRACTION( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(-0.559429644), + static_cast(0.3)), + static_cast(5), + tol_inv); + + // df=5, x=1.475884049, p=0.9 + BOOST_CHECK_CLOSE_FRACTION( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(1.475884049), + static_cast(0.9)), + static_cast(5), + tol_inv); + + // df=5, x=-1.475884049, p=0.1 + BOOST_CHECK_CLOSE_FRACTION( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(-1.475884049), + static_cast(0.1)), + static_cast(5), + tol_inv); + + // df=25, x=-5.2410429995425, p=0.00001 + BOOST_CHECK_CLOSE_FRACTION( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(-5.2410429995425), + static_cast(0.00001)), + static_cast(25), + tol_inv); + + // Symmetry: positive x, p > 0.5 should give the same df. + // df=2, x=+6.96455673428326, p=0.99 + BOOST_CHECK_CLOSE_FRACTION( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(6.96455673428326), + static_cast(0.99)), + static_cast(2), + tol_inv); + + // Domain error: p outside (0,1) +#ifndef BOOST_NO_EXCEPTIONS + BOOST_MATH_CHECK_THROW( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(2), + static_cast(-0.1)), + std::domain_error); + BOOST_MATH_CHECK_THROW( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(2), + static_cast(1.1)), + std::domain_error); +#endif + } // invert_probability_with_respect_to_degrees_of_freedom + // Test for large degrees of freedom when should be same as normal. RealType inf = std::numeric_limits::infinity(); RealType nan = std::numeric_limits::quiet_NaN(); From 7b015ecbf902fca00419ddb5a44b5d6ea212b898 Mon Sep 17 00:00:00 2001 From: dschmitz89 Date: Sun, 12 Apr 2026 16:19:28 +0200 Subject: [PATCH 02/18] Add early exit if approximation is precise to machine precision --- include/boost/math/distributions/students_t.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/include/boost/math/distributions/students_t.hpp b/include/boost/math/distributions/students_t.hpp index 830a181cf5..8ab0d2a28d 100644 --- a/include/boost/math/distributions/students_t.hpp +++ b/include/boost/math/distributions/students_t.hpp @@ -440,6 +440,8 @@ BOOST_MATH_GPU_ENABLED RealType students_t_distribution::inver RealType exact_cdf = cdf(t_hat, x_abs); RealType residual = (exact_cdf > p_adj) ? (exact_cdf - p_adj) : (p_adj - exact_cdf); RealType relative_residual = (p_adj != 0) ? residual / p_adj : residual; + if (relative_residual <= tools::epsilon() * 4) + return nu_hat; // Edgeworth estimate is already exact to machine precision. if (relative_residual <= static_cast(0.1)) hint = nu_hat; } From a1609a9d816bd48aa03d30635b60a8b64bde42dd Mon Sep 17 00:00:00 2001 From: dschmitz89 Date: Sun, 12 Apr 2026 16:19:45 +0200 Subject: [PATCH 03/18] Loosen tolerance and add tests in the tails --- test/test_students_t.cpp | 65 ++++++++++++++++++++++++++-------------- 1 file changed, 42 insertions(+), 23 deletions(-) diff --git a/test/test_students_t.cpp b/test/test_students_t.cpp index fa113ca602..2810566260 100644 --- a/test/test_students_t.cpp +++ b/test/test_students_t.cpp @@ -552,64 +552,83 @@ void test_spots(RealType) // Each case is derived from the CDF spot tests above: the exact df is known, // so we verify that inverting CDF(x; df) = p recovers df to tight tolerance. { - RealType tol_inv = static_cast(1000) * boost::math::tools::epsilon(); - - // df=2, x=-6.96455673428326, p=0.01 - BOOST_CHECK_CLOSE_FRACTION( + // 100 times higher tolerance as in the quantile spot tests above as + // we use root finding + RealType tol_inv = boost::math::tools::epsilon() * 500000; + if (boost::math::tools::digits() > 100) + tol_inv *= 50000; + BOOST_CHECK_CLOSE( students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( static_cast(-6.96455673428326), static_cast(0.01)), static_cast(2), tol_inv); - - // df=5, x=-3.36492999890721, p=0.01 - BOOST_CHECK_CLOSE_FRACTION( + BOOST_CHECK_CLOSE( students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( static_cast(-3.36492999890721), static_cast(0.01)), static_cast(5), tol_inv); - - // df=5, x=-0.559429644, p=0.3 - BOOST_CHECK_CLOSE_FRACTION( + BOOST_CHECK_CLOSE( students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( static_cast(-0.559429644), static_cast(0.3)), static_cast(5), tol_inv); - - // df=5, x=1.475884049, p=0.9 - BOOST_CHECK_CLOSE_FRACTION( + BOOST_CHECK_CLOSE( students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( static_cast(1.475884049), static_cast(0.9)), static_cast(5), tol_inv); - - // df=5, x=-1.475884049, p=0.1 - BOOST_CHECK_CLOSE_FRACTION( + BOOST_CHECK_CLOSE( students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( static_cast(-1.475884049), static_cast(0.1)), static_cast(5), tol_inv); - - // df=25, x=-5.2410429995425, p=0.00001 - BOOST_CHECK_CLOSE_FRACTION( + BOOST_CHECK_CLOSE( students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( static_cast(-5.2410429995425), static_cast(0.00001)), static_cast(25), tol_inv); - - // Symmetry: positive x, p > 0.5 should give the same df. - // df=2, x=+6.96455673428326, p=0.99 - BOOST_CHECK_CLOSE_FRACTION( + BOOST_CHECK_CLOSE( students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( static_cast(6.96455673428326), static_cast(0.99)), static_cast(2), tol_inv); + BOOST_CHECK_CLOSE( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(2), + static_cast(0.610822886098362)), + static_cast(0.1), + tol_inv); + BOOST_CHECK_CLOSE( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(2), + static_cast(0.777242554908434)), + static_cast(0.5), + tol_inv); + BOOST_CHECK_CLOSE( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(2), + static_cast(0.822925875908677)), + static_cast(0.75), + tol_inv); + BOOST_CHECK_CLOSE( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(2), + static_cast(0.977114826753374)), + static_cast(1000), + tol_inv); + BOOST_CHECK_CLOSE( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(2), + static_cast(0.97724973307434)), + static_cast(1e6), + tol_inv); // Domain error: p outside (0,1) #ifndef BOOST_NO_EXCEPTIONS From 1cf9ea2a0c155733e18e1660a1eb80a776fc138e Mon Sep 17 00:00:00 2001 From: dschmitz89 Date: Thu, 16 Apr 2026 21:26:54 +0200 Subject: [PATCH 04/18] Modularize and simplify --- .../boost/math/distributions/students_t.hpp | 163 +++++++++--------- 1 file changed, 82 insertions(+), 81 deletions(-) diff --git a/include/boost/math/distributions/students_t.hpp b/include/boost/math/distributions/students_t.hpp index 8ab0d2a28d..cdf2d44853 100644 --- a/include/boost/math/distributions/students_t.hpp +++ b/include/boost/math/distributions/students_t.hpp @@ -329,6 +329,81 @@ struct cdf_to_df_func bool comp; }; +// +// Shared root-finding helper used by both find_degrees_of_freedom and +// invert_probability_with_respect_to_degrees_of_freedom. +// +template +BOOST_MATH_GPU_ENABLED RealType solve_for_degrees_of_freedom( + Func f, + RealType hint, + bool rising, + const char* function, + Policy const&) +{ + tools::eps_tolerance tol(policies::digits()); + boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); + boost::math::pair r = tools::bracket_and_solve_root( + f, hint, RealType(2), rising, tol, max_iter, Policy()); + RealType result = r.first + (r.second - r.first) / 2; + if (max_iter >= policies::get_max_root_iterations()) + { + return policies::raise_evaluation_error(function, // LCOV_EXCL_LINE + "Unable to locate solution in a reasonable time: either there is no answer to " // LCOV_EXCL_LINE + "the degrees of freedom or the answer is infinite. Current best guess is %1%", // LCOV_EXCL_LINE + result, Policy()); // LCOV_EXCL_LINE + } + return result; +} + +// +// Edgeworth warm-start for invert_probability_with_respect_to_degrees_of_freedom. +// Returns the best df estimate from the quadratic Edgeworth approximation, +// or a small fallback value (0.01) when no positive root is found. +// +template +BOOST_MATH_GPU_ENABLED RealType approximate_df_with_edgeworth_expansion(RealType x_abs, RealType p_adj) +{ + BOOST_MATH_STD_USING + // F(x; nu) ~ cdf_normal(x) - pdf_normal(x)*(x + x^3)/(4*nu) + pdf_normal(x)*(3x + 5x^3 + 7x^5 - 3x^7)/(96*nu^2) + // Substituting u = 1/nu and setting F(x; nu) = p gives b*u^2 - a*u + c = 0, where: + // a = pdf_normal(x) * (x + x^3) / 4 + // b = pdf_normal(x) * (3x + 5x^3 + 7x^5 - 3x^7) / 96 + // c = cdf_normal(x) - p + normal_distribution std_normal(0, 1); + RealType pdf_val = pdf(std_normal, x_abs); + RealType cdf_val = cdf(std_normal, x_abs); + + RealType x2 = x_abs * x_abs; + RealType x3 = x2 * x_abs; + RealType x5 = x3 * x2; + RealType x7 = x5 * x2; + + RealType a = pdf_val * (x_abs + x3) / 4; + RealType b = pdf_val * (3 * x_abs + 5 * x3 + 7 * x5 - 3 * x7) / 96; + RealType c = cdf_val - p_adj; + + RealType discriminant = a * a - 4 * b * c; + if (discriminant >= 0 && b != 0) + { + RealType sqrt_disc = sqrt(discriminant); + // Pick the smallest positive u (= largest df = 1/u). + RealType u1 = (a - sqrt_disc) / (2 * b); + RealType u2 = (a + sqrt_disc) / (2 * b); + RealType u = -1; + if (u1 > 0 && u2 > 0) + u = (u1 < u2) ? u1 : u2; + else if (u1 > 0) + u = u1; + else if (u2 > 0) + u = u2; + + if (u > 0) + return 1 / u; + } + return static_cast(0.01); // fallback +} + } // namespace detail template @@ -353,16 +428,7 @@ BOOST_MATH_GPU_ENABLED RealType students_t_distribution::find_ hint = 1; detail::sample_size_func f(alpha, beta, sd, difference_from_mean); - tools::eps_tolerance tol(policies::digits()); - boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); - boost::math::pair r = tools::bracket_and_solve_root(f, hint, RealType(2), false, tol, max_iter, Policy()); - RealType result = r.first + (r.second - r.first) / 2; - if(max_iter >= policies::get_max_root_iterations()) - { - return policies::raise_evaluation_error(function, "Unable to locate solution in a reasonable time: either there is no answer to how many degrees of freedom are required" // LCOV_EXCL_LINE - " or the answer is infinite. Current best guess is %1%", result, Policy()); // LCOV_EXCL_LINE - } - return result; + return detail::solve_for_degrees_of_freedom(f, hint, false, function, Policy()); } template @@ -394,80 +460,15 @@ BOOST_MATH_GPU_ENABLED RealType students_t_distribution::inver p, Policy()); } - // - // Step 1: Edgeworth warm start. - // F(x; nu) ~ Phi(x) - phi(x)(x + x^3)/(4*nu) + phi(x)(3x + 5x^3 + 7x^5 - 3x^7)/(96*nu^2) - // Substituting u = 1/nu and setting F(x; nu) = p gives: - // B*u^2 - A*u + (Phi(x) - p) = 0 - // where: - // A = phi(x) * (x + x^3) / 4 - // B = phi(x) * (3x + 5x^3 + 7x^5 - 3x^7) / 96 - // - normal_distribution std_normal(0, 1); - RealType phi = pdf(std_normal, x_abs); - RealType Phi = cdf(std_normal, x_abs); - - RealType x2 = x_abs * x_abs; - RealType x3 = x2 * x_abs; - RealType x5 = x3 * x2; - RealType x7 = x5 * x2; - - RealType A = phi * (x_abs + x3) / 4; - RealType B = phi * (3 * x_abs + 5 * x3 + 7 * x5 - 3 * x7) / 96; - RealType C = Phi - p_adj; + // Edgeworth warm start: falls back to 0.01 if no positive root was found. + RealType hint = detail::approximate_df_with_edgeworth_expansion(x_abs, p_adj); - RealType hint = static_cast(0.01); - RealType discriminant = A * A - 4 * B * C; - if (discriminant >= 0 && B != 0) - { - RealType sqrt_disc = sqrt(discriminant); - // Two roots of B*u^2 - A*u + C = 0; pick the smallest positive u (largest nu = 1/u). - RealType u1 = (A - sqrt_disc) / (2 * B); - RealType u2 = (A + sqrt_disc) / (2 * B); - RealType u = -1; - if (u1 > 0 && u2 > 0) - u = (u1 < u2) ? u1 : u2; - else if (u1 > 0) - u = u1; - else if (u2 > 0) - u = u2; - - if (u > 0) - { - RealType nu_hat = 1 / u; - // Step 2: validate by checking relative residual of the exact CDF. - students_t_distribution t_hat(nu_hat); - RealType exact_cdf = cdf(t_hat, x_abs); - RealType residual = (exact_cdf > p_adj) ? (exact_cdf - p_adj) : (p_adj - exact_cdf); - RealType relative_residual = (p_adj != 0) ? residual / p_adj : residual; - if (relative_residual <= tools::epsilon() * 4) - return nu_hat; // Edgeworth estimate is already exact to machine precision. - if (relative_residual <= static_cast(0.1)) - hint = nu_hat; - } - } - - // - // Step 3: root-find with bracket_and_solve_root on f(df) = CDF(x; df) - p. - // For x > 0, CDF(x; df) is strictly increasing in df, so f is a rising function. - // Pass min(p_adj, q_adj) always and use the complement CDF when p_adj > q_adj - // to avoid catastrophic cancellation near 1 (mirrors non_centrality_finder_f). - // + // Root-find on f(df) = CDF(x; df) - p. For x > 0, CDF(x; df) is strictly + // increasing in df. Pass min(p_adj, q_adj) and use the complement CDF when + // p_adj > q_adj to avoid cancellation near 1. RealType q_adj = 1 - p_adj; detail::cdf_to_df_func f(x_abs, p_adj < q_adj ? p_adj : q_adj, p_adj >= q_adj); - tools::eps_tolerance tol(policies::digits()); - boost::math::uintmax_t max_iter = policies::get_max_root_iterations(); - boost::math::pair r = tools::bracket_and_solve_root( - f, hint, RealType(2), true, tol, max_iter, Policy()); - RealType result = r.first + (r.second - r.first) / 2; - if (max_iter >= policies::get_max_root_iterations()) - { - return policies::raise_evaluation_error(function, - "Unable to locate solution in a reasonable time: either there is no answer to " - "the degrees of freedom or the answer is infinite. Current best guess is %1%", - result, Policy()); - } - return result; + return detail::solve_for_degrees_of_freedom(f, hint, true, function, Policy()); } template From 113224160926589acef9e04f80942090f2c921b5 Mon Sep 17 00:00:00 2001 From: dschmitz89 Date: Thu, 16 Apr 2026 21:43:33 +0200 Subject: [PATCH 05/18] Use overflow error instead of inf --- include/boost/math/distributions/students_t.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/boost/math/distributions/students_t.hpp b/include/boost/math/distributions/students_t.hpp index cdf2d44853..982a2f1f42 100644 --- a/include/boost/math/distributions/students_t.hpp +++ b/include/boost/math/distributions/students_t.hpp @@ -453,7 +453,7 @@ BOOST_MATH_GPU_ENABLED RealType students_t_distribution::inver { // CDF(0; df) = 0.5 for all df; only solvable when p == 0.5. if (p_adj == static_cast(0.5)) - return boost::math::numeric_limits::infinity(); + return policies::raise_overflow_error(function, 0, Policy()); return policies::raise_domain_error( function, "No degrees of freedom can satisfy CDF(0; df) == %1% (must be 0.5).", From 123fdd117b497c4d99e7bbca20db8890b0bd9779 Mon Sep 17 00:00:00 2001 From: dschmitz89 Date: Thu, 16 Apr 2026 21:43:41 +0200 Subject: [PATCH 06/18] Relax all test tolerances --- test/test_students_t.cpp | 33 ++++++++++++++------------------- 1 file changed, 14 insertions(+), 19 deletions(-) diff --git a/test/test_students_t.cpp b/test/test_students_t.cpp index 2810566260..9803f714f2 100644 --- a/test/test_students_t.cpp +++ b/test/test_students_t.cpp @@ -551,84 +551,79 @@ void test_spots(RealType) // Tests for invert_probability_with_respect_to_degrees_of_freedom. // Each case is derived from the CDF spot tests above: the exact df is known, // so we verify that inverting CDF(x; df) = p recovers df to tight tolerance. - { - // 100 times higher tolerance as in the quantile spot tests above as - // we use root finding - RealType tol_inv = boost::math::tools::epsilon() * 500000; - if (boost::math::tools::digits() > 100) - tol_inv *= 50000; - BOOST_CHECK_CLOSE( + RealType tol_inv_df = static_cast(0.01); // 0.01% relative error + BOOST_CHECK_CLOSE( students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( static_cast(-6.96455673428326), static_cast(0.01)), static_cast(2), - tol_inv); + tol_inv_df); BOOST_CHECK_CLOSE( students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( static_cast(-3.36492999890721), static_cast(0.01)), static_cast(5), - tol_inv); + tol_inv_df); BOOST_CHECK_CLOSE( students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( static_cast(-0.559429644), static_cast(0.3)), static_cast(5), - tol_inv); + tol_inv_df); BOOST_CHECK_CLOSE( students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( static_cast(1.475884049), static_cast(0.9)), static_cast(5), - tol_inv); + tol_inv_df); BOOST_CHECK_CLOSE( students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( static_cast(-1.475884049), static_cast(0.1)), static_cast(5), - tol_inv); + tol_inv_df); BOOST_CHECK_CLOSE( students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( static_cast(-5.2410429995425), static_cast(0.00001)), static_cast(25), - tol_inv); + tol_inv_df); BOOST_CHECK_CLOSE( students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( static_cast(6.96455673428326), static_cast(0.99)), static_cast(2), - tol_inv); + tol_inv_df); BOOST_CHECK_CLOSE( students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( static_cast(2), static_cast(0.610822886098362)), static_cast(0.1), - tol_inv); + tol_inv_df); BOOST_CHECK_CLOSE( students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( static_cast(2), static_cast(0.777242554908434)), static_cast(0.5), - tol_inv); + tol_inv_df); BOOST_CHECK_CLOSE( students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( static_cast(2), static_cast(0.822925875908677)), static_cast(0.75), - tol_inv); + tol_inv_df); BOOST_CHECK_CLOSE( students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( static_cast(2), static_cast(0.977114826753374)), static_cast(1000), - tol_inv); + tol_inv_df); BOOST_CHECK_CLOSE( students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( static_cast(2), static_cast(0.97724973307434)), static_cast(1e6), - tol_inv); + tol_inv_df); // Domain error: p outside (0,1) #ifndef BOOST_NO_EXCEPTIONS From 24c88badf59dcfb7d3612f91c129ea293161ac82 Mon Sep 17 00:00:00 2001 From: dschmitz89 Date: Thu, 16 Apr 2026 21:58:31 +0200 Subject: [PATCH 07/18] Should use a linter .. --- test/test_students_t.cpp | 163 +++++++++++++++++++-------------------- 1 file changed, 81 insertions(+), 82 deletions(-) diff --git a/test/test_students_t.cpp b/test/test_students_t.cpp index 9803f714f2..11982f6629 100644 --- a/test/test_students_t.cpp +++ b/test/test_students_t.cpp @@ -553,92 +553,91 @@ void test_spots(RealType) // so we verify that inverting CDF(x; df) = p recovers df to tight tolerance. RealType tol_inv_df = static_cast(0.01); // 0.01% relative error BOOST_CHECK_CLOSE( - students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( - static_cast(-6.96455673428326), - static_cast(0.01)), + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(-6.96455673428326), + static_cast(0.01)), + static_cast(2), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(-3.36492999890721), + static_cast(0.01)), + static_cast(5), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(-0.559429644), + static_cast(0.3)), + static_cast(5), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(1.475884049), + static_cast(0.9)), + static_cast(5), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(-1.475884049), + static_cast(0.1)), + static_cast(5), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(-5.2410429995425), + static_cast(0.00001)), + static_cast(25), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(6.96455673428326), + static_cast(0.99)), + static_cast(2), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( static_cast(2), - tol_inv_df); - BOOST_CHECK_CLOSE( - students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( - static_cast(-3.36492999890721), - static_cast(0.01)), - static_cast(5), - tol_inv_df); - BOOST_CHECK_CLOSE( - students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( - static_cast(-0.559429644), - static_cast(0.3)), - static_cast(5), - tol_inv_df); - BOOST_CHECK_CLOSE( - students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( - static_cast(1.475884049), - static_cast(0.9)), - static_cast(5), - tol_inv_df); - BOOST_CHECK_CLOSE( - students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( - static_cast(-1.475884049), - static_cast(0.1)), - static_cast(5), - tol_inv_df); - BOOST_CHECK_CLOSE( - students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( - static_cast(-5.2410429995425), - static_cast(0.00001)), - static_cast(25), - tol_inv_df); - BOOST_CHECK_CLOSE( - students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( - static_cast(6.96455673428326), - static_cast(0.99)), + static_cast(0.610822886098362)), + static_cast(0.1), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(2), + static_cast(0.777242554908434)), + static_cast(0.5), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( static_cast(2), - tol_inv_df); - BOOST_CHECK_CLOSE( - students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( - static_cast(2), - static_cast(0.610822886098362)), - static_cast(0.1), - tol_inv_df); - BOOST_CHECK_CLOSE( - students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( - static_cast(2), - static_cast(0.777242554908434)), - static_cast(0.5), - tol_inv_df); - BOOST_CHECK_CLOSE( - students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( - static_cast(2), - static_cast(0.822925875908677)), - static_cast(0.75), - tol_inv_df); - BOOST_CHECK_CLOSE( - students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( - static_cast(2), - static_cast(0.977114826753374)), - static_cast(1000), - tol_inv_df); - BOOST_CHECK_CLOSE( - students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( - static_cast(2), - static_cast(0.97724973307434)), - static_cast(1e6), - tol_inv_df); - - // Domain error: p outside (0,1) + static_cast(0.822925875908677)), + static_cast(0.75), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(2), + static_cast(0.977114826753374)), + static_cast(1000), + tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(2), + static_cast(0.97724973307434)), + static_cast(1e6), + tol_inv_df); + + // Domain error: p outside (0,1) #ifndef BOOST_NO_EXCEPTIONS - BOOST_MATH_CHECK_THROW( - students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( - static_cast(2), - static_cast(-0.1)), - std::domain_error); - BOOST_MATH_CHECK_THROW( - students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( - static_cast(2), - static_cast(1.1)), - std::domain_error); + BOOST_MATH_CHECK_THROW( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(2), + static_cast(-0.1)), + std::domain_error); + BOOST_MATH_CHECK_THROW( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(2), + static_cast(1.1)), + std::domain_error); #endif - } // invert_probability_with_respect_to_degrees_of_freedom // Test for large degrees of freedom when should be same as normal. RealType inf = std::numeric_limits::infinity(); From 0c8365a5905b419c00ef3b689c5a46af49215cdf Mon Sep 17 00:00:00 2001 From: dschmitz89 Date: Thu, 16 Apr 2026 22:49:27 +0200 Subject: [PATCH 08/18] Relax test tolerance for float precision and remove ill conditioned test with df=1e6 --- test/test_students_t.cpp | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/test/test_students_t.cpp b/test/test_students_t.cpp index 11982f6629..e7f7324b2b 100644 --- a/test/test_students_t.cpp +++ b/test/test_students_t.cpp @@ -551,7 +551,11 @@ void test_spots(RealType) // Tests for invert_probability_with_respect_to_degrees_of_freedom. // Each case is derived from the CDF spot tests above: the exact df is known, // so we verify that inverting CDF(x; df) = p recovers df to tight tolerance. - RealType tol_inv_df = static_cast(0.01); // 0.01% relative error + // float has ~7 significant digits; large-df inversion is ill-conditioned at that precision, + // so use a looser tolerance for single precision. + RealType tol_inv_df = std::is_same::value + ? static_cast(0.1) // float: 0.1% + : static_cast(0.01); // double+: 0.01% BOOST_CHECK_CLOSE( students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( static_cast(-6.96455673428326), @@ -618,12 +622,12 @@ void test_spots(RealType) static_cast(0.977114826753374)), static_cast(1000), tol_inv_df); - BOOST_CHECK_CLOSE( - students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( - static_cast(2), - static_cast(0.97724973307434)), - static_cast(1e6), - tol_inv_df); + BOOST_CHECK_CLOSE( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(1), + static_cast(0.8413326478347855)), + static_cast(1e4), + tol_inv_df); // Domain error: p outside (0,1) #ifndef BOOST_NO_EXCEPTIONS From 1868814d3d76175f45048a0841bb06f5bce09965 Mon Sep 17 00:00:00 2001 From: dschmitz89 Date: Thu, 16 Apr 2026 22:49:44 +0200 Subject: [PATCH 09/18] Fix indentation --- test/test_students_t.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/test_students_t.cpp b/test/test_students_t.cpp index e7f7324b2b..a617fbc1fc 100644 --- a/test/test_students_t.cpp +++ b/test/test_students_t.cpp @@ -622,12 +622,12 @@ void test_spots(RealType) static_cast(0.977114826753374)), static_cast(1000), tol_inv_df); - BOOST_CHECK_CLOSE( - students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( - static_cast(1), - static_cast(0.8413326478347855)), + BOOST_CHECK_CLOSE( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(1), + static_cast(0.8413326478347855)), static_cast(1e4), - tol_inv_df); + tol_inv_df); // Domain error: p outside (0,1) #ifndef BOOST_NO_EXCEPTIONS From b97dceb189c5a62182f7226fc5420ad12f875ba3 Mon Sep 17 00:00:00 2001 From: dschmitz89 Date: Thu, 16 Apr 2026 23:14:40 +0200 Subject: [PATCH 10/18] Loosen tolerance for very high degrees of freedom --- test/test_students_t.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/test_students_t.cpp b/test/test_students_t.cpp index a617fbc1fc..07b02433b9 100644 --- a/test/test_students_t.cpp +++ b/test/test_students_t.cpp @@ -627,8 +627,7 @@ void test_spots(RealType) static_cast(1), static_cast(0.8413326478347855)), static_cast(1e4), - tol_inv_df); - + tol_inv_df * static_cast(5)); // Looser tolerance for ill-conditioned problem with very large df // Domain error: p outside (0,1) #ifndef BOOST_NO_EXCEPTIONS BOOST_MATH_CHECK_THROW( From 9a26afcb4d83c67eaa579a1ac7b4f969ac35703c Mon Sep 17 00:00:00 2001 From: dschmitz89 Date: Sat, 18 Apr 2026 12:19:23 +0200 Subject: [PATCH 11/18] Test if approximation is a good hint, otherwise fall back to 0.01 --- .../boost/math/distributions/students_t.hpp | 49 ++++++++++++------- 1 file changed, 30 insertions(+), 19 deletions(-) diff --git a/include/boost/math/distributions/students_t.hpp b/include/boost/math/distributions/students_t.hpp index 982a2f1f42..c1154ff927 100644 --- a/include/boost/math/distributions/students_t.hpp +++ b/include/boost/math/distributions/students_t.hpp @@ -285,6 +285,11 @@ BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type -BOOST_MATH_GPU_ENABLED RealType approximate_df_with_edgeworth_expansion(RealType x_abs, RealType p_adj) +BOOST_MATH_GPU_ENABLED bool approximate_df_with_edgeworth_expansion(RealType x_abs, RealType p_adj, RealType& result) { BOOST_MATH_STD_USING // F(x; nu) ~ cdf_normal(x) - pdf_normal(x)*(x + x^3)/(4*nu) + pdf_normal(x)*(3x + 5x^3 + 7x^5 - 3x^7)/(96*nu^2) @@ -387,21 +392,17 @@ BOOST_MATH_GPU_ENABLED RealType approximate_df_with_edgeworth_expansion(RealType if (discriminant >= 0 && b != 0) { RealType sqrt_disc = sqrt(discriminant); - // Pick the smallest positive u (= largest df = 1/u). - RealType u1 = (a - sqrt_disc) / (2 * b); - RealType u2 = (a + sqrt_disc) / (2 * b); - RealType u = -1; - if (u1 > 0 && u2 > 0) - u = (u1 < u2) ? u1 : u2; - else if (u1 > 0) - u = u1; - else if (u2 > 0) - u = u2; - + // The two roots of b*u^2 - a*u + c = 0. Pick the smallest positive u (= largest df = 1/u). + RealType u = (a - sqrt_disc) / (2 * b); + if (u <= 0) + u = (a + sqrt_disc) / (2 * b); if (u > 0) - return 1 / u; + { + result = 1 / u; + return true; + } } - return static_cast(0.01); // fallback + return false; } } // namespace detail @@ -460,9 +461,19 @@ BOOST_MATH_GPU_ENABLED RealType students_t_distribution::inver p, Policy()); } - // Edgeworth warm start: falls back to 0.01 if no positive root was found. - RealType hint = detail::approximate_df_with_edgeworth_expansion(x_abs, p_adj); - + // Edgeworth warm start: compute a df estimate; fall back to df_hint_fallback if it fails + // or is too inaccurate. + RealType hint = static_cast(detail::df_hint_fallback); + if (detail::approximate_df_with_edgeworth_expansion(x_abs, p_adj, hint)) + { + // Check that approximation is at least somewhat close; + // for small degrees of freedom it does not fail but is very inaccurate. + students_t_distribution t_approx(hint); + RealType exact_cdf = cdf(t_approx, x_abs); + RealType relative_error = fabs(exact_cdf - p_adj) / p_adj; + if (relative_error > static_cast(0.1)) + hint = static_cast(detail::df_hint_fallback); + } // Root-find on f(df) = CDF(x; df) - p. For x > 0, CDF(x; df) is strictly // increasing in df. Pass min(p_adj, q_adj) and use the complement CDF when // p_adj > q_adj to avoid cancellation near 1. From eecf3d021e280219e815af228b464f05e683cc1b Mon Sep 17 00:00:00 2001 From: dschmitz89 Date: Sat, 18 Apr 2026 12:19:40 +0200 Subject: [PATCH 12/18] Add test coverage of all cases for approximation --- test/test_students_t.cpp | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/test/test_students_t.cpp b/test/test_students_t.cpp index 07b02433b9..37fec8f772 100644 --- a/test/test_students_t.cpp +++ b/test/test_students_t.cpp @@ -628,6 +628,20 @@ void test_spots(RealType) static_cast(0.8413326478347855)), static_cast(1e4), tol_inv_df * static_cast(5)); // Looser tolerance for ill-conditioned problem with very large df + // Small df case 1: Edgeworth expansion breaks down for small degrees if freedom, use fallback + BOOST_CHECK_CLOSE( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(1e20), + static_cast(0.5244796002843015)), + static_cast(1e-3), + tol_inv_df); + // Small df case 2: Edgeworth expansion succeeds but is inaccurate, use fallback + BOOST_CHECK_CLOSE( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(2.0), + static_cast(0.500000000644961)), + static_cast(1e-10), + tol_inv_df); // Domain error: p outside (0,1) #ifndef BOOST_NO_EXCEPTIONS BOOST_MATH_CHECK_THROW( From 0308d3c678442348a127da06b9cf25d435cc3d7e Mon Sep 17 00:00:00 2001 From: dschmitz89 Date: Sat, 18 Apr 2026 13:01:45 +0200 Subject: [PATCH 13/18] Do not use symmetry for better numerical accuracy in the tails --- .../boost/math/distributions/students_t.hpp | 44 +++++++++---------- 1 file changed, 20 insertions(+), 24 deletions(-) diff --git a/include/boost/math/distributions/students_t.hpp b/include/boost/math/distributions/students_t.hpp index c1154ff927..02f7497fbc 100644 --- a/include/boost/math/distributions/students_t.hpp +++ b/include/boost/math/distributions/students_t.hpp @@ -367,7 +367,7 @@ BOOST_MATH_GPU_ENABLED RealType solve_for_degrees_of_freedom( // Returns false when no positive root is found; 'result' is left unchanged. // template -BOOST_MATH_GPU_ENABLED bool approximate_df_with_edgeworth_expansion(RealType x_abs, RealType p_adj, RealType& result) +BOOST_MATH_GPU_ENABLED bool approximate_df_with_edgeworth_expansion(RealType x, RealType p, RealType& result) { BOOST_MATH_STD_USING // F(x; nu) ~ cdf_normal(x) - pdf_normal(x)*(x + x^3)/(4*nu) + pdf_normal(x)*(3x + 5x^3 + 7x^5 - 3x^7)/(96*nu^2) @@ -376,17 +376,17 @@ BOOST_MATH_GPU_ENABLED bool approximate_df_with_edgeworth_expansion(RealType x_a // b = pdf_normal(x) * (3x + 5x^3 + 7x^5 - 3x^7) / 96 // c = cdf_normal(x) - p normal_distribution std_normal(0, 1); - RealType pdf_val = pdf(std_normal, x_abs); - RealType cdf_val = cdf(std_normal, x_abs); + RealType pdf_val = pdf(std_normal, x); + RealType cdf_val = cdf(std_normal, x); - RealType x2 = x_abs * x_abs; - RealType x3 = x2 * x_abs; + RealType x2 = x * x; + RealType x3 = x2 * x; RealType x5 = x3 * x2; RealType x7 = x5 * x2; - RealType a = pdf_val * (x_abs + x3) / 4; - RealType b = pdf_val * (3 * x_abs + 5 * x3 + 7 * x5 - 3 * x7) / 96; - RealType c = cdf_val - p_adj; + RealType a = pdf_val * (x + x3) / 4; + RealType b = pdf_val * (3 * x + 5 * x3 + 7 * x5 - 3 * x7) / 96; + RealType c = cdf_val - p; RealType discriminant = a * a - 4 * b * c; if (discriminant >= 0 && b != 0) @@ -446,14 +446,10 @@ BOOST_MATH_GPU_ENABLED RealType students_t_distribution::inver if (false == detail::check_probability(function, p, &error_result, Policy())) return error_result; - // Use symmetry: CDF(-x; df) = 1 - CDF(x; df), so always work with x >= 0. - RealType x_abs = (x < 0) ? -x : x; - RealType p_adj = (x < 0) ? (1 - p) : p; - - if (x_abs == 0) + if (x == 0) { // CDF(0; df) = 0.5 for all df; only solvable when p == 0.5. - if (p_adj == static_cast(0.5)) + if (p == static_cast(0.5)) return policies::raise_overflow_error(function, 0, Policy()); return policies::raise_domain_error( function, @@ -462,24 +458,24 @@ BOOST_MATH_GPU_ENABLED RealType students_t_distribution::inver } // Edgeworth warm start: compute a df estimate; fall back to df_hint_fallback if it fails - // or is too inaccurate. + // or is too inaccurate RealType hint = static_cast(detail::df_hint_fallback); - if (detail::approximate_df_with_edgeworth_expansion(x_abs, p_adj, hint)) + if (detail::approximate_df_with_edgeworth_expansion(x, p, hint)) { // Check that approximation is at least somewhat close; // for small degrees of freedom it does not fail but is very inaccurate. students_t_distribution t_approx(hint); - RealType exact_cdf = cdf(t_approx, x_abs); - RealType relative_error = fabs(exact_cdf - p_adj) / p_adj; + RealType exact_cdf = cdf(t_approx, x); + RealType relative_error = fabs(exact_cdf - p) / p; if (relative_error > static_cast(0.1)) hint = static_cast(detail::df_hint_fallback); } - // Root-find on f(df) = CDF(x; df) - p. For x > 0, CDF(x; df) is strictly - // increasing in df. Pass min(p_adj, q_adj) and use the complement CDF when - // p_adj > q_adj to avoid cancellation near 1. - RealType q_adj = 1 - p_adj; - detail::cdf_to_df_func f(x_abs, p_adj < q_adj ? p_adj : q_adj, p_adj >= q_adj); - return detail::solve_for_degrees_of_freedom(f, hint, true, function, Policy()); + // Root-find on f(df) = CDF(x; df) - p. + // CDF(x; df) is strictly increasing in df for x > 0, decreasing for x < 0. + // Use the smaller of p and q=1-p to avoid cancellation near 1. + RealType q = 1 - p; + detail::cdf_to_df_func f(x, p < q ? p : q, p < q ? false : true); + return detail::solve_for_degrees_of_freedom(f, hint, x > 0, function, Policy()); } template From 4957f2e8bfba424052db063b466acdca6a0fe49e Mon Sep 17 00:00:00 2001 From: dschmitz89 Date: Sat, 18 Apr 2026 13:04:43 +0200 Subject: [PATCH 14/18] Add edge case tests --- test/test_students_t.cpp | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/test/test_students_t.cpp b/test/test_students_t.cpp index 37fec8f772..4d5b42ba84 100644 --- a/test/test_students_t.cpp +++ b/test/test_students_t.cpp @@ -654,6 +654,18 @@ void test_spots(RealType) static_cast(2), static_cast(1.1)), std::domain_error); + // x == 0, p != 0.5: no df can satisfy CDF(0; df) == p → domain error + BOOST_MATH_CHECK_THROW( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(0), + static_cast(0.3)), + std::domain_error); + // x == 0, p == 0.5: CDF(0; df) == 0.5 for all df → overflow (infinite solutions) + BOOST_MATH_CHECK_THROW( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(0), + static_cast(0.5)), + std::overflow_error); #endif // Test for large degrees of freedom when should be same as normal. From 9c865a3181a2a0bcd93fa40f0096a3cd00c4ff88 Mon Sep 17 00:00:00 2001 From: dschmitz89 Date: Sun, 19 Apr 2026 09:41:43 +0200 Subject: [PATCH 15/18] Exclude float from testing for extreme values --- test/test_students_t.cpp | 32 ++++++++++++++++++-------------- 1 file changed, 18 insertions(+), 14 deletions(-) diff --git a/test/test_students_t.cpp b/test/test_students_t.cpp index 4d5b42ba84..e85fccd729 100644 --- a/test/test_students_t.cpp +++ b/test/test_students_t.cpp @@ -628,20 +628,24 @@ void test_spots(RealType) static_cast(0.8413326478347855)), static_cast(1e4), tol_inv_df * static_cast(5)); // Looser tolerance for ill-conditioned problem with very large df - // Small df case 1: Edgeworth expansion breaks down for small degrees if freedom, use fallback - BOOST_CHECK_CLOSE( - students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( - static_cast(1e20), - static_cast(0.5244796002843015)), - static_cast(1e-3), - tol_inv_df); - // Small df case 2: Edgeworth expansion succeeds but is inaccurate, use fallback - BOOST_CHECK_CLOSE( - students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( - static_cast(2.0), - static_cast(0.500000000644961)), - static_cast(1e-10), - tol_inv_df); + // Small df edge cases where approximation becomes problematic: only meaningful for double and above. + if (!std::is_same::value) + { + // Small df case 1: Edgeworth expansion breaks down for small degrees of freedom, use fallback + BOOST_CHECK_CLOSE( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(1e20), + static_cast(0.5244796002843015)), + static_cast(1e-3), + tol_inv_df); + // Small df case 2: Edgeworth expansion succeeds but is inaccurate, use fallback + BOOST_CHECK_CLOSE( + students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( + static_cast(2.0), + static_cast(0.500000000644961)), + static_cast(1e-10), + tol_inv_df); + } // Domain error: p outside (0,1) #ifndef BOOST_NO_EXCEPTIONS BOOST_MATH_CHECK_THROW( From 80de92258b098d651020403eaa4a7cd9e5042459 Mon Sep 17 00:00:00 2001 From: dschmitz89 Date: Sun, 19 Apr 2026 09:45:20 +0200 Subject: [PATCH 16/18] Update test comment --- test/test_students_t.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/test_students_t.cpp b/test/test_students_t.cpp index e85fccd729..b90e6b72e9 100644 --- a/test/test_students_t.cpp +++ b/test/test_students_t.cpp @@ -628,7 +628,8 @@ void test_spots(RealType) static_cast(0.8413326478347855)), static_cast(1e4), tol_inv_df * static_cast(5)); // Looser tolerance for ill-conditioned problem with very large df - // Small df edge cases where approximation becomes problematic: only meaningful for double and above. + // Small df edge cases where approximation becomes problematic: would need different values for float + // For now just test double and long double to have test overage of these code paths if (!std::is_same::value) { // Small df case 1: Edgeworth expansion breaks down for small degrees of freedom, use fallback From 2540e5f17f9aabd3ea34c6c577b31f0cb42f886f Mon Sep 17 00:00:00 2001 From: dschmitz89 Date: Sun, 19 Apr 2026 10:51:47 +0200 Subject: [PATCH 17/18] Replace non ASCII character.. --- test/test_students_t.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_students_t.cpp b/test/test_students_t.cpp index b90e6b72e9..90d4935931 100644 --- a/test/test_students_t.cpp +++ b/test/test_students_t.cpp @@ -659,13 +659,13 @@ void test_spots(RealType) static_cast(2), static_cast(1.1)), std::domain_error); - // x == 0, p != 0.5: no df can satisfy CDF(0; df) == p → domain error + // x == 0, p != 0.5: no df can satisfy CDF(0; df) == p -> domain error BOOST_MATH_CHECK_THROW( students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( static_cast(0), static_cast(0.3)), std::domain_error); - // x == 0, p == 0.5: CDF(0; df) == 0.5 for all df → overflow (infinite solutions) + // x == 0, p == 0.5: CDF(0; df) == 0.5 for all df -> overflow (infinite solutions) BOOST_MATH_CHECK_THROW( students_t_distribution::invert_probability_with_respect_to_degrees_of_freedom( static_cast(0), From 739df4b9eca717a70fd4c5a1043455c6cc9b4751 Mon Sep 17 00:00:00 2001 From: dschmitz89 Date: Mon, 20 Apr 2026 18:26:20 +0200 Subject: [PATCH 18/18] Simplify using boost's own relative error function --- include/boost/math/distributions/students_t.hpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/include/boost/math/distributions/students_t.hpp b/include/boost/math/distributions/students_t.hpp index 02f7497fbc..70584a239f 100644 --- a/include/boost/math/distributions/students_t.hpp +++ b/include/boost/math/distributions/students_t.hpp @@ -19,6 +19,7 @@ #include #include // for ibeta(a, b, x). #include +#include #include #include #include @@ -464,10 +465,10 @@ BOOST_MATH_GPU_ENABLED RealType students_t_distribution::inver { // Check that approximation is at least somewhat close; // for small degrees of freedom it does not fail but is very inaccurate. + RealType relative_error_threshold = static_cast(0.1); students_t_distribution t_approx(hint); - RealType exact_cdf = cdf(t_approx, x); - RealType relative_error = fabs(exact_cdf - p) / p; - if (relative_error > static_cast(0.1)) + RealType p_approx = cdf(t_approx, x); + if (relative_difference(p_approx, p) > relative_error_threshold) hint = static_cast(detail::df_hint_fallback); } // Root-find on f(df) = CDF(x; df) - p.