special_functions/detail/gamma_inva.hpp

0001 //  (C) Copyright John Maddock 2006.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 //
0007 // This is not a complete header file, it is included by gamma.hpp
0008 // after it has defined it's definitions.  This inverts the incomplete
0009 // gamma functions P and Q on the first parameter "a" using a generic
0010 // root finding algorithm (TOMS Algorithm 748).
0011 //
0012
0013 #ifndef BOOST_MATH_SP_DETAIL_GAMMA_INVA
0014 #define BOOST_MATH_SP_DETAIL_GAMMA_INVA
0015
0016 #ifdef _MSC_VER
0017 #pragma once
0018 #endif
0019
0020 #include <cstdint>
0021 #include <boost/math/tools/toms748_solve.hpp>
0022
0023 namespace boost{ namespace math{ namespace detail{
0024
0025 template <class T, class Policy>
0026 struct gamma_inva_t
0027 {
0028    gamma_inva_t(T z_, T p_, bool invert_) : z(z_), p(p_), invert(invert_) {}
0029    T operator()(T a)
0030    {
0031       return invert ? p - boost::math::gamma_q(a, z, Policy()) : boost::math::gamma_p(a, z, Policy()) - p;
0032    }
0033 private:
0034    T z, p;
0035    bool invert;
0036 };
0037
0038 template <class T, class Policy>
0039 T inverse_poisson_cornish_fisher(T lambda, T p, T q, const Policy& pol)
0040 {
0041    BOOST_MATH_STD_USING
0042    // mean:
0043    T m = lambda;
0044    // standard deviation:
0045    T sigma = sqrt(lambda);
0046    // skewness
0047    T sk = 1 / sigma;
0048    // kurtosis:
0049    // T k = 1/lambda;
0050    // Get the inverse of a std normal distribution:
0051    T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p, pol) * constants::root_two<T>();
0052    // Set the sign:
0053    if(p < 0.5)
0054       x = -x;
0055    T x2 = x * x;
0056    // w is correction term due to skewness
0057    T w = x + sk * (x2 - 1) / 6;
0058    /*
0059    // Add on correction due to kurtosis.
0060    // Disabled for now, seems to make things worse?
0061    //
0062    if(lambda >= 10)
0063       w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36;
0064    */
0065    w = m + sigma * w;
0066    return w > tools::min_value<T>() ? w : tools::min_value<T>();
0067 }
0068
0069 template <class T, class Policy>
0070 T gamma_inva_imp(const T& z, const T& p, const T& q, const Policy& pol)
0071 {
0072    BOOST_MATH_STD_USING  // for ADL of std lib math functions
0073    //
0074    // Special cases first:
0075    //
0076    if(p == 0)
0077    {
0078       return policies::raise_overflow_error<T>("boost::math::gamma_p_inva<%1%>(%1%, %1%)", nullptr, Policy());
0079    }
0080    if(q == 0)
0081    {
0082       return tools::min_value<T>();
0083    }
0084    //
0085    // Function object, this is the functor whose root
0086    // we have to solve:
0087    //
0088    gamma_inva_t<T, Policy> f(z, (p < q) ? p : q, (p < q) ? false : true);
0089    //
0090    // Tolerance: full precision.
0091    //
0092    tools::eps_tolerance<T> tol(policies::digits<T, Policy>());
0093    //
0094    // Now figure out a starting guess for what a may be,
0095    // we'll start out with a value that'll put p or q
0096    // right bang in the middle of their range, the functions
0097    // are quite sensitive so we should need too many steps
0098    // to bracket the root from there:
0099    //
0100    T guess;
0101    T factor = 8;
0102    if(z >= 1)
0103    {
0104       //
0105       // We can use the relationship between the incomplete
0106       // gamma function and the poisson distribution to
0107       // calculate an approximate inverse, for large z
0108       // this is actually pretty accurate, but it fails badly
0109       // when z is very small.  Also set our step-factor according
0110       // to how accurate we think the result is likely to be:
0111       //
0112       guess = 1 + inverse_poisson_cornish_fisher(z, q, p, pol);
0113       if(z > 5)
0114       {
0115          if(z > 1000)
0116             factor = 1.01f;
0117          else if(z > 50)
0118             factor = 1.1f;
0119          else if(guess > 10)
0120             factor = 1.25f;
0121          else
0122             factor = 2;
0123          if(guess < 1.1)
0124             factor = 8;
0125       }
0126    }
0127    else if(z > 0.5)
0128    {
0129       guess = z * 1.2f;
0130    }
0131    else
0132    {
0133       guess = -0.4f / log(z);
0134    }
0135    //
0136    // Max iterations permitted:
0137    //
0138    std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
0139    //
0140    // Use our generic derivative-free root finding procedure.
0141    // We could use Newton steps here, taking the PDF of the
0142    // Poisson distribution as our derivative, but that's
0143    // even worse performance-wise than the generic method :-(
0144    //
0145    std::pair<T, T> r = bracket_and_solve_root(f, guess, factor, false, tol, max_iter, pol);
0146    if(max_iter >= policies::get_max_root_iterations<Policy>())
0147       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);
0148    return (r.first + r.second) / 2;
0149 }
0150
0151 } // namespace detail
0152
0153 template <class T1, class T2, class Policy>
0154 inline typename tools::promote_args<T1, T2>::type
0155    gamma_p_inva(T1 x, T2 p, const Policy& pol)
0156 {
0157    typedef typename tools::promote_args<T1, T2>::type result_type;
0158    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0159    typedef typename policies::normalise<
0160       Policy,
0161       policies::promote_float<false>,
0162       policies::promote_double<false>,
0163       policies::discrete_quantile<>,
0164       policies::assert_undefined<> >::type forwarding_policy;
0165
0166    if(p == 0)
0167    {
0168       policies::raise_overflow_error<result_type>("boost::math::gamma_p_inva<%1%>(%1%, %1%)", nullptr, Policy());
0169    }
0170    if(p == 1)
0171    {
0172       return tools::min_value<result_type>();
0173    }
0174
0175    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
0176       detail::gamma_inva_imp(
0177          static_cast<value_type>(x),
0178          static_cast<value_type>(p),
0179          static_cast<value_type>(1 - static_cast<value_type>(p)),
0180          pol), "boost::math::gamma_p_inva<%1%>(%1%, %1%)");
0181 }
0182
0183 template <class T1, class T2, class Policy>
0184 inline typename tools::promote_args<T1, T2>::type
0185    gamma_q_inva(T1 x, T2 q, const Policy& pol)
0186 {
0187    typedef typename tools::promote_args<T1, T2>::type result_type;
0188    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0189    typedef typename policies::normalise<
0190       Policy,
0191       policies::promote_float<false>,
0192       policies::promote_double<false>,
0193       policies::discrete_quantile<>,
0194       policies::assert_undefined<> >::type forwarding_policy;
0195
0196    if(q == 1)
0197    {
0198       policies::raise_overflow_error<result_type>("boost::math::gamma_q_inva<%1%>(%1%, %1%)", nullptr, Policy());
0199    }
0200    if(q == 0)
0201    {
0202       return tools::min_value<result_type>();
0203    }
0204
0205    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
0206       detail::gamma_inva_imp(
0207          static_cast<value_type>(x),
0208          static_cast<value_type>(1 - static_cast<value_type>(q)),
0209          static_cast<value_type>(q),
0210          pol), "boost::math::gamma_q_inva<%1%>(%1%, %1%)");
0211 }
0212
0213 template <class T1, class T2>
0214 inline typename tools::promote_args<T1, T2>::type
0215    gamma_p_inva(T1 x, T2 p)
0216 {
0217    return boost::math::gamma_p_inva(x, p, policies::policy<>());
0218 }
0219
0220 template <class T1, class T2>
0221 inline typename tools::promote_args<T1, T2>::type
0222    gamma_q_inva(T1 x, T2 q)
0223 {
0224    return boost::math::gamma_q_inva(x, q, policies::policy<>());
0225 }
0226
0227 } // namespace math
0228 } // namespace boost
0229
0230 #endif // BOOST_MATH_SP_DETAIL_GAMMA_INVA
0231
0232
0233