t_distribution_inv.hpp

来自「Boost provides free peer-reviewed portab」· HPP 代码 · 共 517 行 · 第 1/2 页

            //
            T a = 4 * (u - u * u);//1 - 4 * (u - 0.5f) * (u - 0.5f);
            T b = boost::math::cbrt(a);
            static const T c = 0.85498797333834849467655443627193L;
            T p = 6 * (1 + c * (1 / b - 1));
            T p0;
            do{
               T p2 = p * p;
               T p4 = p2 * p2;
               T p5 = p * p4;
               p0 = p;
               // next term is given by Eq 41:
               p = 2 * (8 * a * p5 - 270 * p2 + 2187) / (5 * (4 * a * p4 - 216 * p - 243));
            }while(fabs((p - p0) / p) > tolerance);
            //
            // Use Eq 45 to extract the result:
            //
            p = sqrt(p - df);
            result = (u - 0.5f) < 0 ? -p : p;
            break;
         }
#if 0
         //
         // These are Shaw's "exact" but iterative solutions
         // for even df, the numerical accuracy of these is
         // rather less than Hill's method, so these are disabled
         // for now, which is a shame because they are reasonably
         // quick to evaluate...
         //
      case 8:
         {
            //
            // Newton-Raphson iteration of a polynomial case,
            // choice of seed value is taken from Shaw's online
            // supplement:
            //
            static const T c8 = 0.85994765706259820318168359251872L;
            T a = 4 * (u - u * u); //1 - 4 * (u - 0.5f) * (u - 0.5f);
            T b = pow(a, T(1) / 4);
            T p = 8 * (1 + c8 * (1 / b - 1));
            T p0 = p;
            do{
               T p5 = p * p;
               p5 *= p5 * p;
               p0 = p;
               // Next term is given by Eq 42:
               p = 2 * (3 * p + (640 * (160 + p * (24 + p * (p + 4)))) / (-5120 + p * (-2048 - 960 * p + a * p5))) / 7;
            }while(fabs((p - p0) / p) > tolerance);
            //
            // Use Eq 45 to extract the result:
            //
            p = sqrt(p - df);
            result = (u - 0.5f) < 0 ? -p : p;
            break;
         }
      case 10:
         {
            //
            // Newton-Raphson iteration of a polynomial case,
            // choice of seed value is taken from Shaw's online
            // supplement:
            //
            static const T c10 = 0.86781292867813396759105692122285L;
            T a = 4 * (u - u * u); //1 - 4 * (u - 0.5f) * (u - 0.5f);
            T b = pow(a, T(1) / 5);
            T p = 10 * (1 + c10 * (1 / b - 1));
            T p0;
            do{
               T p6 = p * p;
               p6 *= p6 * p6;
               p0 = p;
               // Next term given by Eq 43:
               p = (8 * p) / 9 + (218750 * (21875 + 4 * p * (625 + p * (75 + 2 * p * (5 + p))))) /
                  (9 * (-68359375 + 8 * p * (-2343750 + p * (-546875 - 175000 * p + 8 * a * p6))));
            }while(fabs((p - p0) / p) > tolerance);
            //
            // Use Eq 45 to extract the result:
            //
            p = sqrt(p - df);
            result = (u - 0.5f) < 0 ? -p : p;
            break;
         }
#endif
      default:
         goto calculate_real;
      }
   }
   else
   {
calculate_real:
      if(df < 3)
      {
         //
         // Use a roughly linear scheme to choose between Shaw's
         // tail series and body series:
         //
         T crossover = 0.2742f - df * 0.0242143f;
         if(u > crossover)
         {
            result = boost::math::detail::inverse_students_t_body_series(df, u, pol);
         }
         else
         {
            result = boost::math::detail::inverse_students_t_tail_series(df, u, pol);
         }
      }
      else
      {
         //
         // Use Hill's method except in the exteme tails
         // where we use Shaw's tail series.
         // The crossover point is roughly exponential in -df:
         //
         T crossover = ldexp(1.0f, iround(df / -0.654f, pol));
         if(u > crossover)
         {
            result = boost::math::detail::inverse_students_t_hill(df, u, pol);
         }
         else
         {
            result = boost::math::detail::inverse_students_t_tail_series(df, u, pol);
         }
      }
   }
   return invert ? -result : result;
}

template <class T, class Policy>
inline T find_ibeta_inv_from_t_dist(T a, T p, T q, T* py, const Policy& pol)
{
   T u = (p > q) ? 0.5f - q / 2 : p / 2;
   T v = 1 - u; // u < 0.5 so no cancellation error
   T df = a * 2;
   T t = boost::math::detail::inverse_students_t(df, u, v, pol);
   T x = df / (df + t * t);
   *py = t * t / (df + t * t);
   return x;
}

template <class T, class Policy>
inline T fast_students_t_quantile_imp(T df, T p, const Policy& pol, const mpl::false_*)
{
   BOOST_MATH_STD_USING
   //
   // Need to use inverse incomplete beta to get
   // required precision so not so fast:
   //
   T probability = (p > 0.5) ? 1 - p : p;
   T t, x, y;
   x = ibeta_inv(df / 2, T(0.5), 2 * probability, &y, pol);
   if(df * y > tools::max_value<T>() * x)
      t = policies::raise_overflow_error<T>("boost::math::students_t_quantile<%1%>(%1%,%1%)", 0, pol);
   else
      t = sqrt(df * y / x);
   //
   // Figure out sign based on the size of p:
   //
   if(p < 0.5)
      t = -t;
   return t;
}

template <class T, class Policy>
T fast_students_t_quantile_imp(T df, T p, const Policy& pol, const mpl::true_*)
{
   BOOST_MATH_STD_USING
   bool invert = false;
   if((df < 2) && (floor(df) != df))
      return boost::math::detail::fast_students_t_quantile_imp(df, p, pol, static_cast<mpl::false_*>(0));
   if(p > 0.5)
   {
      p = 1 - p;
      invert = true;
   }
   //
   // Get an estimate of the result:
   //
   bool exact;
   T t = inverse_students_t(df, p, 1-p, pol, &exact);
   if((t == 0) || exact)
      return invert ? -t : t; // can't do better!
   //
   // Change variables to inverse incomplete beta:
   //
   T t2 = t * t;
   T xb = df / (df + t2);
   T y = t2 / (df + t2);
   T a = df / 2;
   //
   // t can be so large that x underflows,
   // just return our estimate in that case:
   //
   if(xb == 0)
      return t;
   //
   // Get incomplete beta and it's derivative:
   //
   T f1;
   T f0 = xb < y ? ibeta_imp(a, constants::half<T>(), xb, pol, false, true, &f1)
      : ibeta_imp(constants::half<T>(), a, y, pol, true, true, &f1);

   // Get cdf from incomplete beta result:
   T p0 = f0 / 2  - p;
   // Get pdf from derivative:
   T p1 = f1 * sqrt(y * xb * xb * xb / df);
   //
   // Second derivative divided by p1:
   //
   // yacas gives:
   //
   // In> PrettyForm(Simplify(D(t) (1 + t^2/v) ^ (-(v+1)/2)))
   //
   //  |                        | v + 1     |     |
   //  |                       -| ----- + 1 |     |
   //  |                        |   2       |     |
   // -|             |  2     |                   |
   //  |             | t      |                   |
   //  |             | -- + 1 |                   |
   //  | ( v + 1 ) * | v      |               * t |
   // ---------------------------------------------
   //                       v
   //
   // Which after some manipulation is:
   //
   // -p1 * t * (df + 1) / (t^2 + df)
   //
   T p2 = t * (df + 1) / (t * t + df);
   // Halley step:
   t = fabs(t);
   t += p0 / (p1 + p0 * p2 / 2);
   return !invert ? -t : t;
}

template <class T, class Policy>
inline T fast_students_t_quantile(T df, T p, const Policy& pol)
{
   typedef typename policies::evaluation<T, 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;

   typedef mpl::bool_<
      (std::numeric_limits<T>::digits <= 53)
       &&
      (std::numeric_limits<T>::is_specialized)> tag_type;
   return policies::checked_narrowing_cast<T, forwarding_policy>(fast_students_t_quantile_imp(static_cast<value_type>(df), static_cast<value_type>(p), pol, static_cast<tag_type*>(0)), "boost::math::students_t_quantile<%1%>(%1%,%1%,%1%)");
}

}}} // namespaces

#endif // BOOST_MATH_SF_DETAIL_INV_T_HPP