boost/math/special_functions/detail/ibeta_inv_ab.hpp

   1 //  (C) Copyright John Maddock 2006.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 //
   7 // This is not a complete header file, it is included by beta.hpp
   8 // after it has defined it's definitions.  This inverts the incomplete
   9 // beta functions ibeta and ibetac on the first parameters "a"
  10 // and "b" using a generic root finding algorithm (TOMS Algorithm 748).
  11 //
  12
  13 #ifndef BOOST_MATH_SP_DETAIL_BETA_INV_AB
  14 #define BOOST_MATH_SP_DETAIL_BETA_INV_AB
  15
  16 #ifdef _MSC_VER
  17 #pragma once
  18 #endif
  19
  20 #include <boost/math/tools/toms748_solve.hpp>
  21 #include <boost/cstdint.hpp>
  22
  23 namespace boost{ namespace math{ namespace detail{
  24
  25 template <class T, class Policy>
  26 struct beta_inv_ab_t
  27 {
  28    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_) {}
  29    T operator()(T a)
  30    {
  31       return invert ?
  32          p - boost::math::ibetac(swap_ab ? b : a, swap_ab ? a : b, z, Policy())
  33          : boost::math::ibeta(swap_ab ? b : a, swap_ab ? a : b, z, Policy()) - p;
  34    }
  35 private:
  36    T b, z, p;
  37    bool invert, swap_ab;
  38 };
  39
  40 template <class T, class Policy>
  41 T inverse_negative_binomial_cornish_fisher(T n, T sf, T sfc, T p, T q, const Policy& pol)
  42 {
  43    BOOST_MATH_STD_USING
  44    // mean:
  45    T m = n * (sfc) / sf;
  46    T t = sqrt(n * (sfc));
  47    // standard deviation:
  48    T sigma = t / sf;
  49    // skewness
  50    T sk = (1 + sfc) / t;
  51    // kurtosis:
  52    T k = (6 - sf * (5+sfc)) / (n * (sfc));
  53    // Get the inverse of a std normal distribution:
  54    T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p, pol) * constants::root_two<T>();
  55    // Set the sign:
  56    if(p < 0.5)
  57       x = -x;
  58    T x2 = x * x;
  59    // w is correction term due to skewness
  60    T w = x + sk * (x2 - 1) / 6;
  61    //
  62    // Add on correction due to kurtosis.
  63    //
  64    if(n >= 10)
  65       w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36;
  66
  67    w = m + sigma * w;
  68    if(w < tools::min_value<T>())
  69       return tools::min_value<T>();
  70    return w;
  71 }
  72
  73 template <class T, class Policy>
  74 T ibeta_inv_ab_imp(const T& b, const T& z, const T& p, const T& q, bool swap_ab, const Policy& pol)
  75 {
  76    BOOST_MATH_STD_USING  // for ADL of std lib math functions
  77    //
  78    // Special cases first:
  79    //
  80    BOOST_MATH_INSTRUMENT_CODE("b = " << b << " z = " << z << " p = " << p << " q = " << " swap = " << swap_ab);
  81    if(p == 0)
  82    {
  83       return swap_ab ? tools::min_value<T>() : tools::max_value<T>();
  84    }
  85    if(q == 0)
  86    {
  87       return swap_ab ? tools::max_value<T>() : tools::min_value<T>();
  88    }
  89    //
  90    // Function object, this is the functor whose root
  91    // we have to solve:
  92    //
  93    beta_inv_ab_t<T, Policy> f(b, z, (p < q) ? p : q, (p < q) ? false : true, swap_ab);
  94    //
  95    // Tolerance: full precision.
  96    //
  97    tools::eps_tolerance<T> tol(policies::digits<T, Policy>());
  98    //
  99    // Now figure out a starting guess for what a may be,
 100    // we'll start out with a value that'll put p or q
 101    // right bang in the middle of their range, the functions
 102    // are quite sensitive so we should need too many steps
 103    // to bracket the root from there:
 104    //
 105    T guess = 0;
 106    T factor = 5;
 107    //
 108    // Convert variables to parameters of a negative binomial distribution:
 109    //
 110    T n = b;
 111    T sf = swap_ab ? z : 1-z;
 112    T sfc = swap_ab ? 1-z : z;
 113    T u = swap_ab ? p : q;
 114    T v = swap_ab ? q : p;
 115    if(u <= pow(sf, n))
 116    {
 117       //
 118       // Result is less than 1, negative binomial approximation
 119       // is useless....
 120       //
 121       if((p < q) != swap_ab)
 122       {
 123          guess = (std::min)(T(b * 2), T(1));
 124       }
 125       else
 126       {
 127          guess = (std::min)(T(b / 2), T(1));
 128       }
 129    }
 130    if(n * n * n * u * sf > 0.005)
 131       guess = 1 + inverse_negative_binomial_cornish_fisher(n, sf, sfc, u, v, pol);
 132
 133    if(guess < 10)
 134    {
 135       //
 136       // Negative binomial approximation not accurate in this area:
 137       //
 138       if((p < q) != swap_ab)
 139       {
 140          guess = (std::min)(T(b * 2), T(10));
 141       }
 142       else
 143       {
 144          guess = (std::min)(T(b / 2), T(10));
 145       }
 146    }
 147    else
 148       factor = (v < sqrt(tools::epsilon<T>())) ? 2 : (guess < 20 ? 1.2f : 1.1f);
 149    BOOST_MATH_INSTRUMENT_CODE("guess = " << guess);
 150    //
 151    // Max iterations permitted:
 152    //
 153    boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
 154    std::pair<T, T> r = bracket_and_solve_root(f, guess, factor, swap_ab ? true : false, tol, max_iter, pol);
 155    if(max_iter >= policies::get_max_root_iterations<Policy>())
 156       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);
 157    return (r.first + r.second) / 2;
 158 }
 159
 160 } // namespace detail
 161
 162 template <class RT1, class RT2, class RT3, class Policy>
 163 typename tools::promote_args<RT1, RT2, RT3>::type
 164       ibeta_inva(RT1 b, RT2 x, RT3 p, const Policy& pol)
 165 {
 166    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
 167    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 168    typedef typename policies::normalise<
 169       Policy,
 170       policies::promote_float<false>,
 171       policies::promote_double<false>,
 172       policies::discrete_quantile<>,
 173       policies::assert_undefined<> >::type forwarding_policy;
 174
 175    if(p == 0)
 176    {
 177       return tools::max_value<result_type>();
 178    }
 179    if(p == 1)
 180    {
 181       return tools::min_value<result_type>();
 182    }
 183
 184    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
 185       detail::ibeta_inv_ab_imp(
 186          static_cast<value_type>(b),
 187          static_cast<value_type>(x),
 188          static_cast<value_type>(p),
 189          static_cast<value_type>(1 - static_cast<value_type>(p)),
 190          false, pol),
 191       "boost::math::ibeta_inva<%1%>(%1%,%1%,%1%)");
 192 }
 193
 194 template <class RT1, class RT2, class RT3, class Policy>
 195 typename tools::promote_args<RT1, RT2, RT3>::type
 196       ibetac_inva(RT1 b, RT2 x, RT3 q, const Policy& pol)
 197 {
 198    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
 199    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 200    typedef typename policies::normalise<
 201       Policy,
 202       policies::promote_float<false>,
 203       policies::promote_double<false>,
 204       policies::discrete_quantile<>,
 205       policies::assert_undefined<> >::type forwarding_policy;
 206
 207    if(q == 1)
 208    {
 209       return tools::max_value<result_type>();
 210    }
 211    if(q == 0)
 212    {
 213       return tools::min_value<result_type>();
 214    }
 215
 216    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
 217       detail::ibeta_inv_ab_imp(
 218          static_cast<value_type>(b),
 219          static_cast<value_type>(x),
 220          static_cast<value_type>(1 - static_cast<value_type>(q)),
 221          static_cast<value_type>(q),
 222          false, pol),
 223       "boost::math::ibetac_inva<%1%>(%1%,%1%,%1%)");
 224 }
 225
 226 template <class RT1, class RT2, class RT3, class Policy>
 227 typename tools::promote_args<RT1, RT2, RT3>::type
 228       ibeta_invb(RT1 a, RT2 x, RT3 p, const Policy& pol)
 229 {
 230    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
 231    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 232    typedef typename policies::normalise<
 233       Policy,
 234       policies::promote_float<false>,
 235       policies::promote_double<false>,
 236       policies::discrete_quantile<>,
 237       policies::assert_undefined<> >::type forwarding_policy;
 238
 239    if(p == 0)
 240    {
 241       return tools::min_value<result_type>();
 242    }
 243    if(p == 1)
 244    {
 245       return tools::max_value<result_type>();
 246    }
 247
 248    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
 249       detail::ibeta_inv_ab_imp(
 250          static_cast<value_type>(a),
 251          static_cast<value_type>(x),
 252          static_cast<value_type>(p),
 253          static_cast<value_type>(1 - static_cast<value_type>(p)),
 254          true, pol),
 255       "boost::math::ibeta_invb<%1%>(%1%,%1%,%1%)");
 256 }
 257
 258 template <class RT1, class RT2, class RT3, class Policy>
 259 typename tools::promote_args<RT1, RT2, RT3>::type
 260       ibetac_invb(RT1 a, RT2 x, RT3 q, const Policy& pol)
 261 {
 262    typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
 263    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 264    typedef typename policies::normalise<
 265       Policy,
 266       policies::promote_float<false>,
 267       policies::promote_double<false>,
 268       policies::discrete_quantile<>,
 269       policies::assert_undefined<> >::type forwarding_policy;
 270
 271    if(q == 1)
 272    {
 273       return tools::min_value<result_type>();
 274    }
 275    if(q == 0)
 276    {
 277       return tools::max_value<result_type>();
 278    }
 279
 280    return policies::checked_narrowing_cast<result_type, forwarding_policy>(
 281       detail::ibeta_inv_ab_imp(
 282          static_cast<value_type>(a),
 283          static_cast<value_type>(x),
 284          static_cast<value_type>(1 - static_cast<value_type>(q)),
 285          static_cast<value_type>(q),
 286          true, pol),
 287          "boost::math::ibetac_invb<%1%>(%1%,%1%,%1%)");
 288 }
 289
 290 template <class RT1, class RT2, class RT3>
 291 inline typename tools::promote_args<RT1, RT2, RT3>::type
 292          ibeta_inva(RT1 b, RT2 x, RT3 p)
 293 {
 294    return boost::math::ibeta_inva(b, x, p, policies::policy<>());
 295 }
 296
 297 template <class RT1, class RT2, class RT3>
 298 inline typename tools::promote_args<RT1, RT2, RT3>::type
 299          ibetac_inva(RT1 b, RT2 x, RT3 q)
 300 {
 301    return boost::math::ibetac_inva(b, x, q, policies::policy<>());
 302 }
 303
 304 template <class RT1, class RT2, class RT3>
 305 inline typename tools::promote_args<RT1, RT2, RT3>::type
 306          ibeta_invb(RT1 a, RT2 x, RT3 p)
 307 {
 308    return boost::math::ibeta_invb(a, x, p, policies::policy<>());
 309 }
 310
 311 template <class RT1, class RT2, class RT3>
 312 inline typename tools::promote_args<RT1, RT2, RT3>::type
 313          ibetac_invb(RT1 a, RT2 x, RT3 q)
 314 {
 315    return boost::math::ibetac_invb(a, x, q, policies::policy<>());
 316 }
 317
 318 } // namespace math
 319 } // namespace boost
 320
 321 #endif // BOOST_MATH_SP_DETAIL_BETA_INV_AB
 322
 323
 324