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