242 lines · plain
1// (C) Copyright John Maddock 2006.2// (C) Copyright Matt Borland 2024.3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 7//8// This is not a complete header file, it is included by gamma.hpp9// after it has defined it's definitions. This inverts the incomplete10// gamma functions P and Q on the first parameter "a" using a generic11// root finding algorithm (TOMS Algorithm 748).12//13 14#ifndef BOOST_MATH_SP_DETAIL_GAMMA_INVA15#define BOOST_MATH_SP_DETAIL_GAMMA_INVA16 17#ifdef _MSC_VER18#pragma once19#endif20 21#include <boost/math/tools/config.hpp>22#include <boost/math/tools/toms748_solve.hpp>23 24namespace boost{ namespace math{ 25 26#ifdef BOOST_MATH_HAS_NVRTC27template <typename T, typename Policy>28BOOST_MATH_GPU_ENABLED auto erfc_inv(T x, const Policy&);29#endif30 31namespace detail{32 33template <class T, class Policy>34struct gamma_inva_t35{36 BOOST_MATH_GPU_ENABLED gamma_inva_t(T z_, T p_, bool invert_) : z(z_), p(p_), invert(invert_) {}37 BOOST_MATH_GPU_ENABLED T operator()(T a)38 {39 return invert ? p - boost::math::gamma_q(a, z, Policy()) : boost::math::gamma_p(a, z, Policy()) - p;40 }41private:42 T z, p;43 bool invert;44};45 46template <class T, class Policy>47BOOST_MATH_GPU_ENABLED T inverse_poisson_cornish_fisher(T lambda, T p, T q, const Policy& pol)48{49 BOOST_MATH_STD_USING50 // mean:51 T m = lambda;52 // standard deviation:53 T sigma = sqrt(lambda);54 // skewness55 T sk = 1 / sigma;56 // kurtosis:57 // T k = 1/lambda;58 // Get the inverse of a std normal distribution:59 T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p, pol) * constants::root_two<T>();60 // Set the sign:61 if(p < 0.5)62 x = -x;63 T x2 = x * x;64 // w is correction term due to skewness65 T w = x + sk * (x2 - 1) / 6;66 /*67 // Add on correction due to kurtosis.68 // Disabled for now, seems to make things worse?69 //70 if(lambda >= 10)71 w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36;72 */73 w = m + sigma * w;74 return w > tools::min_value<T>() ? w : tools::min_value<T>();75}76 77template <class T, class Policy>78BOOST_MATH_GPU_ENABLED T gamma_inva_imp(const T& z, const T& p, const T& q, const Policy& pol)79{80 BOOST_MATH_STD_USING // for ADL of std lib math functions81 //82 // Special cases first:83 //84 if(p == 0)85 {86 return policies::raise_overflow_error<T>("boost::math::gamma_p_inva<%1%>(%1%, %1%)", nullptr, Policy());87 }88 if(q == 0)89 {90 return tools::min_value<T>();91 }92 //93 // Function object, this is the functor whose root94 // we have to solve:95 //96 gamma_inva_t<T, Policy> f(z, (p < q) ? p : q, (p < q) ? false : true);97 //98 // Tolerance: full precision.99 //100 tools::eps_tolerance<T> tol(policies::digits<T, Policy>());101 //102 // Now figure out a starting guess for what a may be,103 // we'll start out with a value that'll put p or q104 // right bang in the middle of their range, the functions105 // are quite sensitive so we should need too many steps106 // to bracket the root from there:107 //108 T guess;109 T factor = 8;110 if(z >= 1)111 {112 //113 // We can use the relationship between the incomplete114 // gamma function and the poisson distribution to115 // calculate an approximate inverse, for large z116 // this is actually pretty accurate, but it fails badly117 // when z is very small. Also set our step-factor according118 // to how accurate we think the result is likely to be:119 //120 guess = 1 + inverse_poisson_cornish_fisher(z, q, p, pol);121 if(z > 5)122 {123 if(z > 1000)124 factor = 1.01f;125 else if(z > 50)126 factor = 1.1f;127 else if(guess > 10)128 factor = 1.25f;129 else130 factor = 2;131 if(guess < 1.1)132 factor = 8;133 }134 }135 else if(z > 0.5)136 {137 guess = z * 1.2f;138 }139 else140 {141 guess = -0.4f / log(z);142 }143 //144 // Max iterations permitted:145 //146 std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();147 //148 // Use our generic derivative-free root finding procedure.149 // We could use Newton steps here, taking the PDF of the150 // Poisson distribution as our derivative, but that's151 // even worse performance-wise than the generic method :-(152 //153 std::pair<T, T> r = bracket_and_solve_root(f, guess, factor, false, tol, max_iter, pol);154 if(max_iter >= policies::get_max_root_iterations<Policy>())155 return policies::raise_evaluation_error<T>("boost::math::gamma_p_inva<%1%>(%1%, %1%)", "Unable to locate the root within a reasonable number of iterations, closest approximation so far was %1%", r.first, pol);156 return (r.first + r.second) / 2;157}158 159} // namespace detail160 161template <class T1, class T2, class Policy>162BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2>::type163 gamma_p_inva(T1 x, T2 p, const Policy& pol)164{165 typedef typename tools::promote_args<T1, T2>::type result_type;166 typedef typename policies::evaluation<result_type, Policy>::type value_type;167 typedef typename policies::normalise<168 Policy,169 policies::promote_float<false>,170 policies::promote_double<false>,171 policies::discrete_quantile<>,172 policies::assert_undefined<> >::type forwarding_policy;173 174 if(p == 0)175 {176 policies::raise_overflow_error<result_type>("boost::math::gamma_p_inva<%1%>(%1%, %1%)", nullptr, Policy());177 }178 if(p == 1)179 {180 return tools::min_value<result_type>();181 }182 183 return policies::checked_narrowing_cast<result_type, forwarding_policy>(184 detail::gamma_inva_imp(185 static_cast<value_type>(x),186 static_cast<value_type>(p),187 static_cast<value_type>(1 - static_cast<value_type>(p)),188 pol), "boost::math::gamma_p_inva<%1%>(%1%, %1%)");189}190 191template <class T1, class T2, class Policy>192BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2>::type193 gamma_q_inva(T1 x, T2 q, const Policy& pol)194{195 typedef typename tools::promote_args<T1, T2>::type result_type;196 typedef typename policies::evaluation<result_type, Policy>::type value_type;197 typedef typename policies::normalise<198 Policy,199 policies::promote_float<false>,200 policies::promote_double<false>,201 policies::discrete_quantile<>,202 policies::assert_undefined<> >::type forwarding_policy;203 204 if(q == 1)205 {206 policies::raise_overflow_error<result_type>("boost::math::gamma_q_inva<%1%>(%1%, %1%)", nullptr, Policy());207 }208 if(q == 0)209 {210 return tools::min_value<result_type>();211 }212 213 return policies::checked_narrowing_cast<result_type, forwarding_policy>(214 detail::gamma_inva_imp(215 static_cast<value_type>(x),216 static_cast<value_type>(1 - static_cast<value_type>(q)),217 static_cast<value_type>(q),218 pol), "boost::math::gamma_q_inva<%1%>(%1%, %1%)");219}220 221template <class T1, class T2>222BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2>::type223 gamma_p_inva(T1 x, T2 p)224{225 return boost::math::gamma_p_inva(x, p, policies::policy<>());226}227 228template <class T1, class T2>229BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T1, T2>::type230 gamma_q_inva(T1 x, T2 q)231{232 return boost::math::gamma_q_inva(x, q, policies::policy<>());233}234 235} // namespace math236} // namespace boost237 238#endif // BOOST_MATH_SP_DETAIL_GAMMA_INVA239 240 241 242