X-Git-Url: https://git.donarmstrong.com/?p=rsem.git;a=blobdiff_plain;f=boost%2Fmath%2Fspecial_functions%2Fdetail%2Fibeta_inv_ab.hpp;fp=boost%2Fmath%2Fspecial_functions%2Fdetail%2Fibeta_inv_ab.hpp;h=8318a28454285cb53f57e1f4a623dd5bd1d77363;hp=0000000000000000000000000000000000000000;hb=2d71eb92104693ca9baa5a2e1c23eeca776d8fd3;hpb=da57529b92adbb7ae74a89861cb39fb35ac7c62d

diff --git a/boost/math/special_functions/detail/ibeta_inv_ab.hpp b/boost/math/special_functions/detail/ibeta_inv_ab.hpp
new file mode 100644
index 0000000..8318a28
--- /dev/null
+++ b/boost/math/special_functions/detail/ibeta_inv_ab.hpp
@@ -0,0 +1,324 @@
+//  (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 beta.hpp
+// after it has defined it's definitions.  This inverts the incomplete
+// beta functions ibeta and ibetac on the first parameters "a"
+// and "b" using a generic root finding algorithm (TOMS Algorithm 748).
+//
+
+#ifndef BOOST_MATH_SP_DETAIL_BETA_INV_AB
+#define BOOST_MATH_SP_DETAIL_BETA_INV_AB
+
+#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 beta_inv_ab_t
+{
+   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_) {}
+   T operator()(T a)
+   {
+      return invert ? 
+         p - boost::math::ibetac(swap_ab ? b : a, swap_ab ? a : b, z, Policy()) 
+         : boost::math::ibeta(swap_ab ? b : a, swap_ab ? a : b, z, Policy()) - p;
+   }
+private:
+   T b, z, p;
+   bool invert, swap_ab;
+};
+
+template <class T, class Policy>
+T inverse_negative_binomial_cornish_fisher(T n, T sf, T sfc, T p, T q, const Policy& pol)
+{
+   BOOST_MATH_STD_USING
+   // mean:
+   T m = n * (sfc) / sf;
+   T t = sqrt(n * (sfc));
+   // standard deviation:
+   T sigma = t / sf;
+   // skewness
+   T sk = (1 + sfc) / t;
+   // kurtosis:
+   T k = (6 - sf * (5+sfc)) / (n * (sfc));
+   // 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.
+   //
+   if(n >= 10)
+      w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36;
+
+   w = m + sigma * w;
+   if(w < tools::min_value<T>())
+      return tools::min_value<T>();
+   return w;
+}
+
+template <class T, class Policy>
+T ibeta_inv_ab_imp(const T& b, const T& z, const T& p, const T& q, bool swap_ab, const Policy& pol)
+{
+   BOOST_MATH_STD_USING  // for ADL of std lib math functions
+   //
+   // Special cases first:
+   //
+   BOOST_MATH_INSTRUMENT_CODE("b = " << b << " z = " << z << " p = " << p << " q = " << " swap = " << swap_ab);
+   if(p == 0)
+   {
+      return swap_ab ? tools::min_value<T>() : tools::max_value<T>();
+   }
+   if(q == 0)
+   {
+      return swap_ab ? tools::max_value<T>() : tools::min_value<T>();
+   }
+   //
+   // Function object, this is the functor whose root
+   // we have to solve:
+   //
+   beta_inv_ab_t<T, Policy> f(b, z, (p < q) ? p : q, (p < q) ? false : true, swap_ab);
+   //
+   // 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 = 0;
+   T factor = 5;
+   //
+   // Convert variables to parameters of a negative binomial distribution:
+   //
+   T n = b;
+   T sf = swap_ab ? z : 1-z;
+   T sfc = swap_ab ? 1-z : z;
+   T u = swap_ab ? p : q;
+   T v = swap_ab ? q : p;
+   if(u <= pow(sf, n))
+   {
+      //
+      // Result is less than 1, negative binomial approximation
+      // is useless....
+      //
+      if((p < q) != swap_ab)
+      {
+         guess = (std::min)(T(b * 2), T(1));
+      }
+      else
+      {
+         guess = (std::min)(T(b / 2), T(1));
+      }
+   }
+   if(n * n * n * u * sf > 0.005)
+      guess = 1 + inverse_negative_binomial_cornish_fisher(n, sf, sfc, u, v, pol);
+
+   if(guess < 10)
+   {
+      //
+      // Negative binomial approximation not accurate in this area:
+      //
+      if((p < q) != swap_ab)
+      {
+         guess = (std::min)(T(b * 2), T(10));
+      }
+      else
+      {
+         guess = (std::min)(T(b / 2), T(10));
+      }
+   }
+   else
+      factor = (v < sqrt(tools::epsilon<T>())) ? 2 : (guess < 20 ? 1.2f : 1.1f);
+   BOOST_MATH_INSTRUMENT_CODE("guess = " << guess);
+   //
+   // Max iterations permitted:
+   //
+   boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
+   std::pair<T, T> r = bracket_and_solve_root(f, guess, factor, swap_ab ? true : false, tol, max_iter, pol);
+   if(max_iter >= policies::get_max_root_iterations<Policy>())
+      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);
+   return (r.first + r.second) / 2;
+}
+
+} // namespace detail
+
+template <class RT1, class RT2, class RT3, class Policy>
+typename tools::promote_args<RT1, RT2, RT3>::type 
+      ibeta_inva(RT1 b, RT2 x, RT3 p, const Policy& pol)
+{
+   typedef typename tools::promote_args<RT1, RT2, RT3>::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::ibeta_inv_ab_imp(
+         static_cast<value_type>(b), 
+         static_cast<value_type>(x), 
+         static_cast<value_type>(p), 
+         static_cast<value_type>(1 - static_cast<value_type>(p)), 
+         false, pol), 
+      "boost::math::ibeta_inva<%1%>(%1%,%1%,%1%)");
+}
+
+template <class RT1, class RT2, class RT3, class Policy>
+typename tools::promote_args<RT1, RT2, RT3>::type 
+      ibetac_inva(RT1 b, RT2 x, RT3 q, const Policy& pol)
+{
+   typedef typename tools::promote_args<RT1, RT2, RT3>::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::ibeta_inv_ab_imp(
+         static_cast<value_type>(b), 
+         static_cast<value_type>(x), 
+         static_cast<value_type>(1 - static_cast<value_type>(q)), 
+         static_cast<value_type>(q), 
+         false, pol),
+      "boost::math::ibetac_inva<%1%>(%1%,%1%,%1%)");
+}
+
+template <class RT1, class RT2, class RT3, class Policy>
+typename tools::promote_args<RT1, RT2, RT3>::type 
+      ibeta_invb(RT1 a, RT2 x, RT3 p, const Policy& pol)
+{
+   typedef typename tools::promote_args<RT1, RT2, RT3>::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::min_value<result_type>();
+   }
+   if(p == 1)
+   {
+      return tools::max_value<result_type>();
+   }
+
+   return policies::checked_narrowing_cast<result_type, forwarding_policy>(
+      detail::ibeta_inv_ab_imp(
+         static_cast<value_type>(a), 
+         static_cast<value_type>(x), 
+         static_cast<value_type>(p), 
+         static_cast<value_type>(1 - static_cast<value_type>(p)), 
+         true, pol),
+      "boost::math::ibeta_invb<%1%>(%1%,%1%,%1%)");
+}
+
+template <class RT1, class RT2, class RT3, class Policy>
+typename tools::promote_args<RT1, RT2, RT3>::type 
+      ibetac_invb(RT1 a, RT2 x, RT3 q, const Policy& pol)
+{
+   typedef typename tools::promote_args<RT1, RT2, RT3>::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::min_value<result_type>();
+   }
+   if(q == 0)
+   {
+      return tools::max_value<result_type>();
+   }
+
+   return policies::checked_narrowing_cast<result_type, forwarding_policy>(
+      detail::ibeta_inv_ab_imp(
+         static_cast<value_type>(a), 
+         static_cast<value_type>(x), 
+         static_cast<value_type>(1 - static_cast<value_type>(q)), 
+         static_cast<value_type>(q), 
+         true, pol),
+         "boost::math::ibetac_invb<%1%>(%1%,%1%,%1%)");
+}
+
+template <class RT1, class RT2, class RT3>
+inline typename tools::promote_args<RT1, RT2, RT3>::type 
+         ibeta_inva(RT1 b, RT2 x, RT3 p)
+{
+   return boost::math::ibeta_inva(b, x, p, policies::policy<>());
+}
+
+template <class RT1, class RT2, class RT3>
+inline typename tools::promote_args<RT1, RT2, RT3>::type 
+         ibetac_inva(RT1 b, RT2 x, RT3 q)
+{
+   return boost::math::ibetac_inva(b, x, q, policies::policy<>());
+}
+
+template <class RT1, class RT2, class RT3>
+inline typename tools::promote_args<RT1, RT2, RT3>::type 
+         ibeta_invb(RT1 a, RT2 x, RT3 p)
+{
+   return boost::math::ibeta_invb(a, x, p, policies::policy<>());
+}
+
+template <class RT1, class RT2, class RT3>
+inline typename tools::promote_args<RT1, RT2, RT3>::type 
+         ibetac_invb(RT1 a, RT2 x, RT3 q)
+{
+   return boost::math::ibetac_invb(a, x, q, policies::policy<>());
+}
+
+} // namespace math
+} // namespace boost
+
+#endif // BOOST_MATH_SP_DETAIL_BETA_INV_AB
+
+
+