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