Chris@16: //  (C) Copyright John Maddock 2006.
Chris@16: //  Use, modification and distribution are subject to the
Chris@16: //  Boost Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: //
Chris@16: // This is not a complete header file, it is included by gamma.hpp
Chris@16: // after it has defined it's definitions.  This inverts the incomplete
Chris@16: // gamma functions P and Q on the first parameter "a" using a generic
Chris@16: // root finding algorithm (TOMS Algorithm 748).
Chris@16: //
Chris@16: 
Chris@16: #ifndef BOOST_MATH_SP_DETAIL_GAMMA_INVA
Chris@16: #define BOOST_MATH_SP_DETAIL_GAMMA_INVA
Chris@16: 
Chris@16: #ifdef _MSC_VER
Chris@16: #pragma once
Chris@16: #endif
Chris@16: 
Chris@16: #include <boost/math/tools/toms748_solve.hpp>
Chris@16: #include <boost/cstdint.hpp>
Chris@16: 
Chris@16: namespace boost{ namespace math{ namespace detail{
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: struct gamma_inva_t
Chris@16: {
Chris@16:    gamma_inva_t(T z_, T p_, bool invert_) : z(z_), p(p_), invert(invert_) {}
Chris@16:    T operator()(T a)
Chris@16:    {
Chris@16:       return invert ? p - boost::math::gamma_q(a, z, Policy()) : boost::math::gamma_p(a, z, Policy()) - p;
Chris@16:    }
Chris@16: private:
Chris@16:    T z, p;
Chris@16:    bool invert;
Chris@16: };
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: T inverse_poisson_cornish_fisher(T lambda, T p, T q, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING
Chris@16:    // mean:
Chris@16:    T m = lambda;
Chris@16:    // standard deviation:
Chris@16:    T sigma = sqrt(lambda);
Chris@16:    // skewness
Chris@16:    T sk = 1 / sigma;
Chris@16:    // kurtosis:
Chris@16:    // T k = 1/lambda;
Chris@16:    // Get the inverse of a std normal distribution:
Chris@16:    T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p, pol) * constants::root_two<T>();
Chris@16:    // Set the sign:
Chris@16:    if(p < 0.5)
Chris@16:       x = -x;
Chris@16:    T x2 = x * x;
Chris@16:    // w is correction term due to skewness
Chris@16:    T w = x + sk * (x2 - 1) / 6;
Chris@16:    /*
Chris@16:    // Add on correction due to kurtosis.
Chris@16:    // Disabled for now, seems to make things worse?
Chris@16:    //
Chris@16:    if(lambda >= 10)
Chris@16:       w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36;
Chris@16:    */
Chris@16:    w = m + sigma * w;
Chris@16:    return w > tools::min_value<T>() ? w : tools::min_value<T>();
Chris@16: }
Chris@16: 
Chris@16: template <class T, class Policy>
Chris@16: T gamma_inva_imp(const T& z, const T& p, const T& q, const Policy& pol)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING  // for ADL of std lib math functions
Chris@16:    //
Chris@16:    // Special cases first:
Chris@16:    //
Chris@16:    if(p == 0)
Chris@16:    {
Chris@101:       return policies::raise_overflow_error<T>("boost::math::gamma_p_inva<%1%>(%1%, %1%)", 0, Policy());
Chris@16:    }
Chris@16:    if(q == 0)
Chris@16:    {
Chris@16:       return tools::min_value<T>();
Chris@16:    }
Chris@16:    //
Chris@16:    // Function object, this is the functor whose root
Chris@16:    // we have to solve:
Chris@16:    //
Chris@16:    gamma_inva_t<T, Policy> f(z, (p < q) ? p : q, (p < q) ? false : true);
Chris@16:    //
Chris@16:    // Tolerance: full precision.
Chris@16:    //
Chris@16:    tools::eps_tolerance<T> tol(policies::digits<T, Policy>());
Chris@16:    //
Chris@16:    // Now figure out a starting guess for what a may be, 
Chris@16:    // we'll start out with a value that'll put p or q
Chris@16:    // right bang in the middle of their range, the functions
Chris@16:    // are quite sensitive so we should need too many steps
Chris@16:    // to bracket the root from there:
Chris@16:    //
Chris@16:    T guess;
Chris@16:    T factor = 8;
Chris@16:    if(z >= 1)
Chris@16:    {
Chris@16:       //
Chris@16:       // We can use the relationship between the incomplete 
Chris@16:       // gamma function and the poisson distribution to
Chris@16:       // calculate an approximate inverse, for large z
Chris@16:       // this is actually pretty accurate, but it fails badly
Chris@16:       // when z is very small.  Also set our step-factor according
Chris@16:       // to how accurate we think the result is likely to be:
Chris@16:       //
Chris@16:       guess = 1 + inverse_poisson_cornish_fisher(z, q, p, pol);
Chris@16:       if(z > 5)
Chris@16:       {
Chris@16:          if(z > 1000)
Chris@16:             factor = 1.01f;
Chris@16:          else if(z > 50)
Chris@16:             factor = 1.1f;
Chris@16:          else if(guess > 10)
Chris@16:             factor = 1.25f;
Chris@16:          else
Chris@16:             factor = 2;
Chris@16:          if(guess < 1.1)
Chris@16:             factor = 8;
Chris@16:       }
Chris@16:    }
Chris@16:    else if(z > 0.5)
Chris@16:    {
Chris@16:       guess = z * 1.2f;
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:       guess = -0.4f / log(z);
Chris@16:    }
Chris@16:    //
Chris@16:    // Max iterations permitted:
Chris@16:    //
Chris@16:    boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
Chris@16:    //
Chris@16:    // Use our generic derivative-free root finding procedure.
Chris@16:    // We could use Newton steps here, taking the PDF of the
Chris@16:    // Poisson distribution as our derivative, but that's
Chris@16:    // even worse performance-wise than the generic method :-(
Chris@16:    //
Chris@16:    std::pair<T, T> r = bracket_and_solve_root(f, guess, factor, false, tol, max_iter, pol);
Chris@16:    if(max_iter >= policies::get_max_root_iterations<Policy>())
Chris@101:       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);
Chris@16:    return (r.first + r.second) / 2;
Chris@16: }
Chris@16: 
Chris@16: } // namespace detail
Chris@16: 
Chris@16: template <class T1, class T2, class Policy>
Chris@16: inline typename tools::promote_args<T1, T2>::type 
Chris@16:    gamma_p_inva(T1 x, T2 p, const Policy& pol)
Chris@16: {
Chris@16:    typedef typename tools::promote_args<T1, T2>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:    if(p == 0)
Chris@16:    {
Chris@101:       policies::raise_overflow_error<result_type>("boost::math::gamma_p_inva<%1%>(%1%, %1%)", 0, Policy());
Chris@16:    }
Chris@16:    if(p == 1)
Chris@16:    {
Chris@16:       return tools::min_value<result_type>();
Chris@16:    }
Chris@16: 
Chris@16:    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
Chris@16:       detail::gamma_inva_imp(
Chris@16:          static_cast<value_type>(x), 
Chris@16:          static_cast<value_type>(p), 
Chris@16:          static_cast<value_type>(1 - static_cast<value_type>(p)), 
Chris@16:          pol), "boost::math::gamma_p_inva<%1%>(%1%, %1%)");
Chris@16: }
Chris@16: 
Chris@16: template <class T1, class T2, class Policy>
Chris@16: inline typename tools::promote_args<T1, T2>::type 
Chris@16:    gamma_q_inva(T1 x, T2 q, const Policy& pol)
Chris@16: {
Chris@16:    typedef typename tools::promote_args<T1, T2>::type result_type;
Chris@16:    typedef typename policies::evaluation<result_type, Policy>::type value_type;
Chris@16:    typedef typename policies::normalise<
Chris@16:       Policy, 
Chris@16:       policies::promote_float<false>, 
Chris@16:       policies::promote_double<false>, 
Chris@16:       policies::discrete_quantile<>,
Chris@16:       policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:    if(q == 1)
Chris@16:    {
Chris@101:       policies::raise_overflow_error<result_type>("boost::math::gamma_q_inva<%1%>(%1%, %1%)", 0, Policy());
Chris@16:    }
Chris@16:    if(q == 0)
Chris@16:    {
Chris@16:       return tools::min_value<result_type>();
Chris@16:    }
Chris@16: 
Chris@16:    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
Chris@16:       detail::gamma_inva_imp(
Chris@16:          static_cast<value_type>(x), 
Chris@16:          static_cast<value_type>(1 - static_cast<value_type>(q)), 
Chris@16:          static_cast<value_type>(q), 
Chris@16:          pol), "boost::math::gamma_q_inva<%1%>(%1%, %1%)");
Chris@16: }
Chris@16: 
Chris@16: template <class T1, class T2>
Chris@16: inline typename tools::promote_args<T1, T2>::type 
Chris@16:    gamma_p_inva(T1 x, T2 p)
Chris@16: {
Chris@16:    return boost::math::gamma_p_inva(x, p, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: template <class T1, class T2>
Chris@16: inline typename tools::promote_args<T1, T2>::type
Chris@16:    gamma_q_inva(T1 x, T2 q)
Chris@16: {
Chris@16:    return boost::math::gamma_q_inva(x, q, policies::policy<>());
Chris@16: }
Chris@16: 
Chris@16: } // namespace math
Chris@16: } // namespace boost
Chris@16: 
Chris@16: #endif // BOOST_MATH_SP_DETAIL_GAMMA_INVA
Chris@16: 
Chris@16: 
Chris@16: