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