Chris@16: //  Copyright John Maddock 2006.
Chris@16: //  Copyright Paul A. Bristow 2006, 2012.
Chris@16: //  Copyright Thomas Mang 2012.
Chris@16: 
Chris@16: //  Use, modification and distribution are subject to the
Chris@16: //  Boost Software License, Version 1.0. (See accompanying file
Chris@16: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_STATS_STUDENTS_T_HPP
Chris@16: #define BOOST_STATS_STUDENTS_T_HPP
Chris@16: 
Chris@16: // http://en.wikipedia.org/wiki/Student%27s_t_distribution
Chris@16: // http://www.itl.nist.gov/div898/handbook/eda/section3/eda3664.htm
Chris@16: 
Chris@16: #include <boost/math/distributions/fwd.hpp>
Chris@16: #include <boost/math/special_functions/beta.hpp> // for ibeta(a, b, x).
Chris@16: #include <boost/math/distributions/complement.hpp>
Chris@16: #include <boost/math/distributions/detail/common_error_handling.hpp>
Chris@16: #include <boost/math/distributions/normal.hpp> 
Chris@16: 
Chris@16: #include <utility>
Chris@16: 
Chris@16: #ifdef BOOST_MSVC
Chris@16: # pragma warning(push)
Chris@16: # pragma warning(disable: 4702) // unreachable code (return after domain_error throw).
Chris@16: #endif
Chris@16: 
Chris@16: namespace boost{ namespace math{
Chris@16: 
Chris@16: template <class RealType = double, class Policy = policies::policy<> >
Chris@16: class students_t_distribution
Chris@16: {
Chris@16: public:
Chris@16:    typedef RealType value_type;
Chris@16:    typedef Policy policy_type;
Chris@16: 
Chris@16:    students_t_distribution(RealType df) : df_(df)
Chris@16:    { // Constructor.
Chris@16:       RealType result;
Chris@16:       detail::check_df_gt0_to_inf( // Checks that df > 0 or df == inf.
Chris@16:          "boost::math::students_t_distribution<%1%>::students_t_distribution", df_, &result, Policy());
Chris@16:    } // students_t_distribution
Chris@16: 
Chris@16:    RealType degrees_of_freedom()const
Chris@16:    {
Chris@16:       return df_;
Chris@16:    }
Chris@16: 
Chris@16:    // Parameter estimation:
Chris@16:    static RealType find_degrees_of_freedom(
Chris@16:       RealType difference_from_mean,
Chris@16:       RealType alpha,
Chris@16:       RealType beta,
Chris@16:       RealType sd,
Chris@16:       RealType hint = 100);
Chris@16: 
Chris@16: private:
Chris@16:    // Data member:
Chris@16:    RealType df_;  // degrees of freedom is a real number or +infinity.
Chris@16: };
Chris@16: 
Chris@16: typedef students_t_distribution<double> students_t; // Convenience typedef for double version.
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline const std::pair<RealType, RealType> range(const students_t_distribution<RealType, Policy>& /*dist*/)
Chris@16: { // Range of permissible values for random variable x.
Chris@16:   // NOT including infinity.
Chris@16:    using boost::math::tools::max_value;
Chris@16:    return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>());
Chris@16: }
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline const std::pair<RealType, RealType> support(const students_t_distribution<RealType, Policy>& /*dist*/)
Chris@16: { // Range of supported values for random variable x.
Chris@16:    // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
Chris@16:    using boost::math::tools::max_value;
Chris@16:    return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>());
Chris@16: }
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType pdf(const students_t_distribution<RealType, Policy>& dist, const RealType& x)
Chris@16: {
Chris@16:    BOOST_FPU_EXCEPTION_GUARD
Chris@16:    BOOST_MATH_STD_USING  // for ADL of std functions.
Chris@16: 
Chris@16:    RealType error_result;
Chris@16:    if(false == detail::check_x(
Chris@16:       "boost::math::pdf(const students_t_distribution<%1%>&, %1%)", x, &error_result, Policy()))
Chris@16:       return error_result;
Chris@16:    RealType df = dist.degrees_of_freedom();
Chris@16:    if(false == detail::check_df_gt0_to_inf( // Check that df > 0 or == +infinity.
Chris@16:       "boost::math::pdf(const students_t_distribution<%1%>&, %1%)", df, &error_result, Policy()))
Chris@16:       return error_result;
Chris@16: 
Chris@16:    RealType result;
Chris@16:    if ((boost::math::isinf)(x))
Chris@16:    { // +infinity.
Chris@16:      normal_distribution<RealType, Policy> n(0, 1); 
Chris@16:      result = pdf(n, x);
Chris@16:      return result;
Chris@16:    }
Chris@16:    RealType limit = policies::get_epsilon<RealType, Policy>();
Chris@16:    // Use policies so that if policy requests lower precision, 
Chris@16:    // then get the normal distribution approximation earlier.
Chris@16:    limit = static_cast<RealType>(1) / limit; // 1/eps
Chris@16:    // for 64-bit double 1/eps = 4503599627370496
Chris@16:    if (df > limit)
Chris@16:    { // Special case for really big degrees_of_freedom > 1 / eps 
Chris@16:      // - use normal distribution which is much faster and more accurate.
Chris@16:      normal_distribution<RealType, Policy> n(0, 1); 
Chris@16:      result = pdf(n, x);
Chris@16:    }
Chris@16:    else
Chris@16:    { // 
Chris@16:      RealType basem1 = x * x / df;
Chris@16:      if(basem1 < 0.125)
Chris@16:      {
Chris@16:         result = exp(-boost::math::log1p(basem1, Policy()) * (1+df) / 2);
Chris@16:      }
Chris@16:      else
Chris@16:      {
Chris@16:         result = pow(1 / (1 + basem1), (df + 1) / 2);
Chris@16:      }
Chris@16:      result /= sqrt(df) * boost::math::beta(df / 2, RealType(0.5f), Policy());
Chris@16:    }
Chris@16:    return result;
Chris@16: } // pdf
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType cdf(const students_t_distribution<RealType, Policy>& dist, const RealType& x)
Chris@16: {
Chris@16:    RealType error_result;
Chris@16:    if(false == detail::check_x(
Chris@16:       "boost::math::pdf(const students_t_distribution<%1%>&, %1%)", x, &error_result, Policy()))
Chris@16:       return error_result;
Chris@16:    RealType df = dist.degrees_of_freedom();
Chris@16:    // Error check:
Chris@16: 
Chris@16:    if(false == detail::check_df_gt0_to_inf(  // Check that df > 0 or == +infinity.
Chris@16:       "boost::math::cdf(const students_t_distribution<%1%>&, %1%)", df, &error_result, Policy()))
Chris@16:       return error_result;
Chris@16: 
Chris@16:    if (x == 0)
Chris@16:    { // Special case with exact result.
Chris@16:      return static_cast<RealType>(0.5);
Chris@16:    }
Chris@16:    if ((boost::math::isinf)(x))
Chris@16:    { // +infinity.
Chris@16:      normal_distribution<RealType, Policy> n(0, 1); 
Chris@16:      RealType result = cdf(n, x);
Chris@16:      return result;
Chris@16:    }
Chris@16:    RealType limit = policies::get_epsilon<RealType, Policy>();
Chris@16:    // Use policies so that if policy requests lower precision, 
Chris@16:    // then get the normal distribution approximation earlier.
Chris@16:    limit = static_cast<RealType>(1) / limit; // 1/eps
Chris@16:    // for 64-bit double 1/eps = 4503599627370496
Chris@16:    if (df > limit)
Chris@16:    { // Special case for really big degrees_of_freedom > 1 / eps (perhaps infinite?)
Chris@16:      // - use normal distribution which is much faster and more accurate.
Chris@16:      normal_distribution<RealType, Policy> n(0, 1); 
Chris@16:      RealType result = cdf(n, x);
Chris@16:      return result;
Chris@16:    }
Chris@16:    else
Chris@16:    { // normal df case.
Chris@16:      //
Chris@16:      // Calculate probability of Student's t using the incomplete beta function.
Chris@16:      // probability = ibeta(degrees_of_freedom / 2, 1/2, degrees_of_freedom / (degrees_of_freedom + t*t))
Chris@16:      //
Chris@16:      // However when t is small compared to the degrees of freedom, that formula
Chris@16:      // suffers from rounding error, use the identity formula to work around
Chris@16:      // the problem:
Chris@16:      //
Chris@16:      // I[x](a,b) = 1 - I[1-x](b,a)
Chris@16:      //
Chris@16:      // and:
Chris@16:      //
Chris@16:      //     x = df / (df + t^2)
Chris@16:      //
Chris@16:      // so:
Chris@16:      //
Chris@16:      // 1 - x = t^2 / (df + t^2)
Chris@16:      //
Chris@16:      RealType x2 = x * x;
Chris@16:      RealType probability;
Chris@16:      if(df > 2 * x2)
Chris@16:      {
Chris@16:         RealType z = x2 / (df + x2);
Chris@16:         probability = ibetac(static_cast<RealType>(0.5), df / 2, z, Policy()) / 2;
Chris@16:      }
Chris@16:      else
Chris@16:      {
Chris@16:         RealType z = df / (df + x2);
Chris@16:         probability = ibeta(df / 2, static_cast<RealType>(0.5), z, Policy()) / 2;
Chris@16:      }
Chris@16:      return (x > 0 ? 1   - probability : probability);
Chris@16:   }
Chris@16: } // cdf
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType quantile(const students_t_distribution<RealType, Policy>& dist, const RealType& p)
Chris@16: {
Chris@16:    BOOST_MATH_STD_USING // for ADL of std functions
Chris@16:    //
Chris@16:    // Obtain parameters:
Chris@16:    RealType probability = p;
Chris@16:  
Chris@16:    // Check for domain errors:
Chris@16:    RealType df = dist.degrees_of_freedom();
Chris@16:    static const char* function = "boost::math::quantile(const students_t_distribution<%1%>&, %1%)";
Chris@16:    RealType error_result;
Chris@16:    if(false == (detail::check_df_gt0_to_inf( // Check that df > 0 or == +infinity.
Chris@16:       function, df, &error_result, Policy())
Chris@16:          && detail::check_probability(function, probability, &error_result, Policy())))
Chris@16:       return error_result;
Chris@16:    // Special cases, regardless of degrees_of_freedom.
Chris@16:    if (probability == 0)
Chris@16:       return -policies::raise_overflow_error<RealType>(function, 0, Policy());
Chris@16:    if (probability == 1)
Chris@16:      return policies::raise_overflow_error<RealType>(function, 0, Policy());
Chris@16:    if (probability == static_cast<RealType>(0.5))
Chris@16:      return 0;  //
Chris@16:    //
Chris@16: #if 0
Chris@16:    // This next block is disabled in favour of a faster method than
Chris@16:    // incomplete beta inverse, but code retained for future reference:
Chris@16:    //
Chris@16:    // Calculate quantile of Student's t using the incomplete beta function inverse:
Chris@16:    //
Chris@16:    probability = (probability > 0.5) ? 1 - probability : probability;
Chris@16:    RealType t, x, y;
Chris@16:    x = ibeta_inv(degrees_of_freedom / 2, RealType(0.5), 2 * probability, &y);
Chris@16:    if(degrees_of_freedom * y > tools::max_value<RealType>() * x)
Chris@16:       t = tools::overflow_error<RealType>(function);
Chris@16:    else
Chris@16:       t = sqrt(degrees_of_freedom * y / x);
Chris@16:    //
Chris@16:    // Figure out sign based on the size of p:
Chris@16:    //
Chris@16:    if(p < 0.5)
Chris@16:       t = -t;
Chris@16: 
Chris@16:    return t;
Chris@16: #endif
Chris@16:    //
Chris@16:    // Depending on how many digits RealType has, this may forward
Chris@16:    // to the incomplete beta inverse as above.  Otherwise uses a
Chris@16:    // faster method that is accurate to ~15 digits everywhere
Chris@16:    // and a couple of epsilon at double precision and in the central 
Chris@16:    // region where most use cases will occur...
Chris@16:    //
Chris@16:    return boost::math::detail::fast_students_t_quantile(df, probability, Policy());
Chris@16: } // quantile
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType cdf(const complemented2_type<students_t_distribution<RealType, Policy>, RealType>& c)
Chris@16: {
Chris@16:    return cdf(c.dist, -c.param);
Chris@16: }
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType quantile(const complemented2_type<students_t_distribution<RealType, Policy>, RealType>& c)
Chris@16: {
Chris@16:    return -quantile(c.dist, c.param);
Chris@16: }
Chris@16: 
Chris@16: //
Chris@16: // Parameter estimation follows:
Chris@16: //
Chris@16: namespace detail{
Chris@16: //
Chris@16: // Functors for finding degrees of freedom:
Chris@16: //
Chris@16: template <class RealType, class Policy>
Chris@16: struct sample_size_func
Chris@16: {
Chris@16:    sample_size_func(RealType a, RealType b, RealType s, RealType d)
Chris@16:       : alpha(a), beta(b), ratio(s*s/(d*d)) {}
Chris@16: 
Chris@16:    RealType operator()(const RealType& df)
Chris@16:    {
Chris@16:       if(df <= tools::min_value<RealType>())
Chris@16:       { // 
Chris@16:          return 1;
Chris@16:       }
Chris@16:       students_t_distribution<RealType, Policy> t(df);
Chris@16:       RealType qa = quantile(complement(t, alpha));
Chris@16:       RealType qb = quantile(complement(t, beta));
Chris@16:       qa += qb;
Chris@16:       qa *= qa;
Chris@16:       qa *= ratio;
Chris@16:       qa -= (df + 1);
Chris@16:       return qa;
Chris@16:    }
Chris@16:    RealType alpha, beta, ratio;
Chris@16: };
Chris@16: 
Chris@16: }  // namespace detail
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: RealType students_t_distribution<RealType, Policy>::find_degrees_of_freedom(
Chris@16:       RealType difference_from_mean,
Chris@16:       RealType alpha,
Chris@16:       RealType beta,
Chris@16:       RealType sd,
Chris@16:       RealType hint)
Chris@16: {
Chris@16:    static const char* function = "boost::math::students_t_distribution<%1%>::find_degrees_of_freedom";
Chris@16:    //
Chris@16:    // Check for domain errors:
Chris@16:    //
Chris@16:    RealType error_result;
Chris@16:    if(false == detail::check_probability(
Chris@16:       function, alpha, &error_result, Policy())
Chris@16:          && detail::check_probability(function, beta, &error_result, Policy()))
Chris@16:       return error_result;
Chris@16: 
Chris@16:    if(hint <= 0)
Chris@16:       hint = 1;
Chris@16: 
Chris@16:    detail::sample_size_func<RealType, Policy> f(alpha, beta, sd, difference_from_mean);
Chris@16:    tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
Chris@16:    boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
Chris@16:    std::pair<RealType, RealType> r = tools::bracket_and_solve_root(f, hint, RealType(2), false, tol, max_iter, Policy());
Chris@16:    RealType result = r.first + (r.second - r.first) / 2;
Chris@16:    if(max_iter >= policies::get_max_root_iterations<Policy>())
Chris@16:    {
Chris@101:       return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
Chris@16:          " either there is no answer to how many degrees of freedom are required"
Chris@16:          " or the answer is infinite.  Current best guess is %1%", result, Policy());
Chris@16:    }
Chris@16:    return result;
Chris@16: }
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType mode(const students_t_distribution<RealType, Policy>& /*dist*/)
Chris@16: {
Chris@16:   // Assume no checks on degrees of freedom are useful (unlike mean).
Chris@16:    return 0; // Always zero by definition.
Chris@16: }
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType median(const students_t_distribution<RealType, Policy>& /*dist*/)
Chris@16: {
Chris@16:    // Assume no checks on degrees of freedom are useful (unlike mean).
Chris@16:    return 0; // Always zero by definition.
Chris@16: }
Chris@16: 
Chris@16: // See section 5.1 on moments at  http://en.wikipedia.org/wiki/Student%27s_t-distribution
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType mean(const students_t_distribution<RealType, Policy>& dist)
Chris@16: {  // Revised for https://svn.boost.org/trac/boost/ticket/7177
Chris@16:    RealType df = dist.degrees_of_freedom();
Chris@16:    if(((boost::math::isnan)(df)) || (df <= 1) ) 
Chris@16:    { // mean is undefined for moment <= 1!
Chris@101:       return policies::raise_domain_error<RealType>(
Chris@16:       "boost::math::mean(students_t_distribution<%1%> const&, %1%)",
Chris@16:       "Mean is undefined for degrees of freedom < 1 but got %1%.", df, Policy());
Chris@16:       return std::numeric_limits<RealType>::quiet_NaN();
Chris@16:    }
Chris@16:    return 0;
Chris@16: } // mean
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType variance(const students_t_distribution<RealType, Policy>& dist)
Chris@16: { // http://en.wikipedia.org/wiki/Student%27s_t-distribution
Chris@16:   // Revised for https://svn.boost.org/trac/boost/ticket/7177
Chris@16:   RealType df = dist.degrees_of_freedom();
Chris@16:   if ((boost::math::isnan)(df) || (df <= 2))
Chris@16:   { // NaN or undefined for <= 2.
Chris@101:      return policies::raise_domain_error<RealType>(
Chris@16:       "boost::math::variance(students_t_distribution<%1%> const&, %1%)",
Chris@16:       "variance is undefined for degrees of freedom <= 2, but got %1%.",
Chris@16:       df, Policy());
Chris@16:     return std::numeric_limits<RealType>::quiet_NaN(); // Undefined.
Chris@16:   }
Chris@16:   if ((boost::math::isinf)(df))
Chris@16:   { // +infinity.
Chris@16:     return 1;
Chris@16:   }
Chris@16:   RealType limit = policies::get_epsilon<RealType, Policy>();
Chris@16:   // Use policies so that if policy requests lower precision, 
Chris@16:   // then get the normal distribution approximation earlier.
Chris@16:   limit = static_cast<RealType>(1) / limit; // 1/eps
Chris@16:   // for 64-bit double 1/eps = 4503599627370496
Chris@16:   if (df > limit)
Chris@16:   { // Special case for really big degrees_of_freedom > 1 / eps.
Chris@16:     return 1;
Chris@16:   }
Chris@16:   else
Chris@16:   {
Chris@16:     return df / (df - 2);
Chris@16:   }
Chris@16: } // variance
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType skewness(const students_t_distribution<RealType, Policy>& dist)
Chris@16: {
Chris@16:     RealType df = dist.degrees_of_freedom();
Chris@16:    if( ((boost::math::isnan)(df)) || (dist.degrees_of_freedom() <= 3))
Chris@16:    { // Undefined for moment k = 3.
Chris@101:       return policies::raise_domain_error<RealType>(
Chris@16:          "boost::math::skewness(students_t_distribution<%1%> const&, %1%)",
Chris@16:          "Skewness is undefined for degrees of freedom <= 3, but got %1%.",
Chris@16:          dist.degrees_of_freedom(), Policy());
Chris@16:       return std::numeric_limits<RealType>::quiet_NaN();
Chris@16:    }
Chris@16:    return 0; // For all valid df, including infinity.
Chris@16: } // skewness
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType kurtosis(const students_t_distribution<RealType, Policy>& dist)
Chris@16: {
Chris@16:    RealType df = dist.degrees_of_freedom();
Chris@16:    if(((boost::math::isnan)(df)) || (df <= 4))
Chris@16:    { // Undefined or infinity for moment k = 4.
Chris@101:       return policies::raise_domain_error<RealType>(
Chris@16:        "boost::math::kurtosis(students_t_distribution<%1%> const&, %1%)",
Chris@16:        "Kurtosis is undefined for degrees of freedom <= 4, but got %1%.",
Chris@16:         df, Policy());
Chris@16:         return std::numeric_limits<RealType>::quiet_NaN(); // Undefined.
Chris@16:    }
Chris@16:    if ((boost::math::isinf)(df))
Chris@16:    { // +infinity.
Chris@16:      return 3;
Chris@16:    }
Chris@16:    RealType limit = policies::get_epsilon<RealType, Policy>();
Chris@16:    // Use policies so that if policy requests lower precision, 
Chris@16:    // then get the normal distribution approximation earlier.
Chris@16:    limit = static_cast<RealType>(1) / limit; // 1/eps
Chris@16:    // for 64-bit double 1/eps = 4503599627370496
Chris@16:    if (df > limit)
Chris@16:    { // Special case for really big degrees_of_freedom > 1 / eps.
Chris@16:      return 3;
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:      //return 3 * (df - 2) / (df - 4); re-arranged to
Chris@16:      return 6 / (df - 4) + 3;
Chris@16:    }
Chris@16: } // kurtosis
Chris@16: 
Chris@16: template <class RealType, class Policy>
Chris@16: inline RealType kurtosis_excess(const students_t_distribution<RealType, Policy>& dist)
Chris@16: {
Chris@16:    // see http://mathworld.wolfram.com/Kurtosis.html
Chris@16: 
Chris@16:    RealType df = dist.degrees_of_freedom();
Chris@16:    if(((boost::math::isnan)(df)) || (df <= 4))
Chris@16:    { // Undefined or infinity for moment k = 4.
Chris@101:      return policies::raise_domain_error<RealType>(
Chris@16:        "boost::math::kurtosis_excess(students_t_distribution<%1%> const&, %1%)",
Chris@16:        "Kurtosis_excess is undefined for degrees of freedom <= 4, but got %1%.",
Chris@16:       df, Policy());
Chris@16:      return std::numeric_limits<RealType>::quiet_NaN(); // Undefined.
Chris@16:    }
Chris@16:    if ((boost::math::isinf)(df))
Chris@16:    { // +infinity.
Chris@16:      return 0;
Chris@16:    }
Chris@16:    RealType limit = policies::get_epsilon<RealType, Policy>();
Chris@16:    // Use policies so that if policy requests lower precision, 
Chris@16:    // then get the normal distribution approximation earlier.
Chris@16:    limit = static_cast<RealType>(1) / limit; // 1/eps
Chris@16:    // for 64-bit double 1/eps = 4503599627370496
Chris@16:    if (df > limit)
Chris@16:    { // Special case for really big degrees_of_freedom > 1 / eps.
Chris@16:      return 0;
Chris@16:    }
Chris@16:    else
Chris@16:    {
Chris@16:      return 6 / (df - 4);
Chris@16:    }
Chris@16: }
Chris@16: 
Chris@16: } // namespace math
Chris@16: } // namespace boost
Chris@16: 
Chris@16: #ifdef BOOST_MSVC
Chris@16: # pragma warning(pop)
Chris@16: #endif
Chris@16: 
Chris@16: // This include must be at the end, *after* the accessors
Chris@16: // for this distribution have been defined, in order to
Chris@16: // keep compilers that support two-phase lookup happy.
Chris@16: #include <boost/math/distributions/detail/derived_accessors.hpp>
Chris@16: 
Chris@16: #endif // BOOST_STATS_STUDENTS_T_HPP