special_functions/detail/ibeta_inv_ab.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 beta.hpp
0008 // after it has defined it's definitions.  This inverts the incomplete
0009 // beta functions ibeta and ibetac on the first parameters "a"
0010 // and "b" using a generic root finding algorithm (TOMS Algorithm 748).
0011 //
0012
0013 #ifndef BOOST_MATH_SP_DETAIL_BETA_INV_AB
0014 #define BOOST_MATH_SP_DETAIL_BETA_INV_AB
0015
0016 #ifdef _MSC_VER
0017 #pragma once
0018 #endif
0019
0020 #include <boost/math/tools/config.hpp>
0021 #include <boost/math/tools/toms748_solve.hpp>
0022 #include <boost/math/tools/precision.hpp>
0023 #include <boost/math/tools/tuple.hpp>
0024 #include <boost/math/policies/error_handling.hpp>
0025
0026 namespace boost{ namespace math{ namespace detail{
0027
0028 template <class T, class Policy>
0029 struct beta_inv_ab_t
0030 {
0031    BOOST_MATH_GPU_ENABLED beta_inv_ab_t(T b_, T z_, T p_, bool invert_, bool swap_ab_) : b(b_), z(z_), p(p_), invert(invert_), swap_ab(swap_ab_) {}
0032    BOOST_MATH_GPU_ENABLED T operator()(T a)
0033    {
0034       return invert ?
0035          p - boost::math::ibetac(swap_ab ? b : a, swap_ab ? a : b, z, Policy())
0036          : boost::math::ibeta(swap_ab ? b : a, swap_ab ? a : b, z, Policy()) - p;
0037    }
0038 private:
0039    T b, z, p;
0040    bool invert, swap_ab;
0041 };
0042
0043 template <class T, class Policy>
0044 BOOST_MATH_GPU_ENABLED T inverse_negative_binomial_cornish_fisher(T n, T sf, T sfc, T p, T q, const Policy& pol)
0045 {
0046    BOOST_MATH_STD_USING
0047    // mean:
0048    T m = n * (sfc) / sf;
0049    T t = sqrt(n * (sfc));
0050    // standard deviation:
0051    T sigma = t / sf;
0052    // skewness
0053    T sk = (1 + sfc) / t;
0054    // kurtosis:
0055    T k = (6 - sf * (5+sfc)) / (n * (sfc));
0056    // Get the inverse of a std normal distribution:
0057    T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p, pol) * constants::root_two<T>();
0058    // Set the sign:
0059    if(p < 0.5)
0060       x = -x;
0061    T x2 = x * x;
0062    // w is correction term due to skewness
0063    T w = x + sk * (x2 - 1) / 6;
0064    //
0065    // Add on correction due to kurtosis.
0066    //
0067    if(n >= 10)
0068       w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36;
0069
0070    w = m + sigma * w;
0071    if(w < tools::min_value<T>())
0072       return tools::min_value<T>();
0073    return w;
0074 }
0075
0076 template <class T, class Policy>
0077 BOOST_MATH_GPU_ENABLED T ibeta_inv_ab_imp(const T& b, const T& z, const T& p, const T& q, bool swap_ab, const Policy& pol)
0078 {
0079    BOOST_MATH_STD_USING  // for ADL of std lib math functions
0080    //
0081    // Special cases first:
0082    //
0083    BOOST_MATH_INSTRUMENT_CODE("b = " << b << " z = " << z << " p = " << p << " q = " << " swap = " << swap_ab);
0084    if(p == 0)
0085    {
0086       return swap_ab ? tools::min_value<T>() : tools::max_value<T>();
0087    }
0088    if(q == 0)
0089    {
0090       return swap_ab ? tools::max_value<T>() : tools::min_value<T>();
0091    }
0092    //
0093    // Function object, this is the functor whose root
0094    // we have to solve:
0095    //
0096    beta_inv_ab_t<T, Policy> f(b, z, (p < q) ? p : q, (p < q) ? false : true, swap_ab);
0097    //
0098    // Tolerance: full precision.
0099    //
0100    tools::eps_tolerance<T> tol(policies::digits<T, Policy>());
0101    //
0102    // Now figure out a starting guess for what a may be,
0103    // we'll start out with a value that'll put p or q
0104    // right bang in the middle of their range, the functions
0105    // are quite sensitive so we should need too many steps
0106    // to bracket the root from there:
0107    //
0108    T guess = 0;
0109    T factor = 5;
0110    //
0111    // Convert variables to parameters of a negative binomial distribution:
0112    //
0113    T n = b;
0114    T sf = swap_ab ? z : 1-z;
0115    T sfc = swap_ab ? 1-z : z;
0116    T u = swap_ab ? p : q;
0117    T v = swap_ab ? q : p;
0118    if(u <= pow(sf, n))
0119    {
0120       //
0121       // Result is less than 1, negative binomial approximation
0122       // is useless....
0123       //
0124       if((p < q) != swap_ab)
0125       {
0126          guess = BOOST_MATH_GPU_SAFE_MIN(T(b * 2), T(1));
0127       }
0128       else
0129       {
0130          guess = BOOST_MATH_GPU_SAFE_MIN(T(b / 2), T(1));
0131       }
0132    }
0133    if(n * n * n * u * sf > 0.005)
0134       guess = 1 + inverse_negative_binomial_cornish_fisher(n, sf, sfc, u, v, pol);
0135
0136    if(guess < 10)
0137    {
0138       //
0139       // Negative binomial approximation not accurate in this area:
0140       //
0141       if((p < q) != swap_ab)
0142       {
0143          guess = BOOST_MATH_GPU_SAFE_MIN(T(b * 2), T(10));
0144       }
0145       else
0146       {
0147          guess = BOOST_MATH_GPU_SAFE_MIN(T(b / 2), T(10));
0148       }
0149    }
0150    else
0151       factor = (v < sqrt(tools::epsilon<T>())) ? 2 : (guess < 20 ? 1.2f : 1.1f);
0152    BOOST_MATH_INSTRUMENT_CODE("guess = " << guess);
0153    //
0154    // Max iterations permitted:
0155    //
0156    boost::math::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
0157    boost::math::pair<T, T> r = bracket_and_solve_root(f, guess, factor, swap_ab ? true : false, tol, max_iter, pol);
0158    if(max_iter >= policies::get_max_root_iterations<Policy>())
0159       return policies::raise_evaluation_error<T>("boost::math::ibeta_invab_imp<%1%>(%1%,%1%,%1%)", "Unable to locate the root within a reasonable number of iterations, closest approximation so far was %1%", r.first, pol);
0160    return (r.first + r.second) / 2;
0161 }
0162
0163 } // namespace detail
0164
0165 template <class RT1, class RT2, class RT3, class Policy>
0166 BOOST_MATH_GPU_ENABLED typename tools::promote_args<RT1, RT2, RT3>::type
0167       ibeta_inva(RT1 b, RT2 x, RT3 p, const Policy& pol)
0168 {
0169    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
0170    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0171    typedef typename policies::normalise<
0172       Policy,
0173       policies::promote_float<false>,
0174       policies::promote_double<false>,
0175       policies::discrete_quantile<>,
0176       policies::assert_undefined<> >::type forwarding_policy;
0177
0178    constexpr auto function = "boost::math::ibeta_inva<%1%>(%1%,%1%,%1%)";
0179    if(p == 0)
0180    {
0181       return policies::raise_overflow_error<result_type>(function, 0, Policy());
0182    }
0183    if(p == 1)
0184    {
0185       return tools::min_value<result_type>();
0186    }
0187
0188    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
0189       detail::ibeta_inv_ab_imp(
0190          static_cast<value_type>(b),
0191          static_cast<value_type>(x),
0192          static_cast<value_type>(p),
0193          static_cast<value_type>(1 - static_cast<value_type>(p)),
0194          false, pol),
0195       function);
0196 }
0197
0198 template <class RT1, class RT2, class RT3, class Policy>
0199 BOOST_MATH_GPU_ENABLED typename tools::promote_args<RT1, RT2, RT3>::type
0200       ibetac_inva(RT1 b, RT2 x, RT3 q, const Policy& pol)
0201 {
0202    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
0203    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0204    typedef typename policies::normalise<
0205       Policy,
0206       policies::promote_float<false>,
0207       policies::promote_double<false>,
0208       policies::discrete_quantile<>,
0209       policies::assert_undefined<> >::type forwarding_policy;
0210
0211    constexpr auto function = "boost::math::ibetac_inva<%1%>(%1%,%1%,%1%)";
0212    if(q == 1)
0213    {
0214       return policies::raise_overflow_error<result_type>(function, 0, Policy());
0215    }
0216    if(q == 0)
0217    {
0218       return tools::min_value<result_type>();
0219    }
0220
0221    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
0222       detail::ibeta_inv_ab_imp(
0223          static_cast<value_type>(b),
0224          static_cast<value_type>(x),
0225          static_cast<value_type>(1 - static_cast<value_type>(q)),
0226          static_cast<value_type>(q),
0227          false, pol),
0228       function);
0229 }
0230
0231 template <class RT1, class RT2, class RT3, class Policy>
0232 BOOST_MATH_GPU_ENABLED typename tools::promote_args<RT1, RT2, RT3>::type
0233       ibeta_invb(RT1 a, RT2 x, RT3 p, const Policy& pol)
0234 {
0235    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
0236    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0237    typedef typename policies::normalise<
0238       Policy,
0239       policies::promote_float<false>,
0240       policies::promote_double<false>,
0241       policies::discrete_quantile<>,
0242       policies::assert_undefined<> >::type forwarding_policy;
0243
0244    constexpr auto function = "boost::math::ibeta_invb<%1%>(%1%,%1%,%1%)";
0245    if(p == 0)
0246    {
0247       return tools::min_value<result_type>();
0248    }
0249    if(p == 1)
0250    {
0251       return policies::raise_overflow_error<result_type>(function, 0, Policy());
0252    }
0253
0254    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
0255       detail::ibeta_inv_ab_imp(
0256          static_cast<value_type>(a),
0257          static_cast<value_type>(x),
0258          static_cast<value_type>(p),
0259          static_cast<value_type>(1 - static_cast<value_type>(p)),
0260          true, pol),
0261       function);
0262 }
0263
0264 template <class RT1, class RT2, class RT3, class Policy>
0265 BOOST_MATH_GPU_ENABLED typename tools::promote_args<RT1, RT2, RT3>::type
0266       ibetac_invb(RT1 a, RT2 x, RT3 q, const Policy& pol)
0267 {
0268    constexpr auto function = "boost::math::ibeta_invb<%1%>(%1%, %1%, %1%)";
0269    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
0270    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0271    typedef typename policies::normalise<
0272       Policy,
0273       policies::promote_float<false>,
0274       policies::promote_double<false>,
0275       policies::discrete_quantile<>,
0276       policies::assert_undefined<> >::type forwarding_policy;
0277
0278    if(q == 1)
0279    {
0280       return tools::min_value<result_type>();
0281    }
0282    if(q == 0)
0283    {
0284       return policies::raise_overflow_error<result_type>(function, 0, Policy());
0285    }
0286
0287    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
0288       detail::ibeta_inv_ab_imp(
0289          static_cast<value_type>(a),
0290          static_cast<value_type>(x),
0291          static_cast<value_type>(1 - static_cast<value_type>(q)),
0292          static_cast<value_type>(q),
0293          true, pol),
0294          function);
0295 }
0296
0297 template <class RT1, class RT2, class RT3>
0298 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type
0299          ibeta_inva(RT1 b, RT2 x, RT3 p)
0300 {
0301    return boost::math::ibeta_inva(b, x, p, policies::policy<>());
0302 }
0303
0304 template <class RT1, class RT2, class RT3>
0305 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type
0306          ibetac_inva(RT1 b, RT2 x, RT3 q)
0307 {
0308    return boost::math::ibetac_inva(b, x, q, policies::policy<>());
0309 }
0310
0311 template <class RT1, class RT2, class RT3>
0312 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type
0313          ibeta_invb(RT1 a, RT2 x, RT3 p)
0314 {
0315    return boost::math::ibeta_invb(a, x, p, policies::policy<>());
0316 }
0317
0318 template <class RT1, class RT2, class RT3>
0319 BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<RT1, RT2, RT3>::type
0320          ibetac_invb(RT1 a, RT2 x, RT3 q)
0321 {
0322    return boost::math::ibetac_invb(a, x, q, policies::policy<>());
0323 }
0324
0325 } // namespace math
0326 } // namespace boost
0327
0328 #endif // BOOST_MATH_SP_DETAIL_BETA_INV_AB
0329
0330
0331