[Boost-commit] svn:boost r79911 - in trunk/boost/math/distributions: . detail

Subject: [Boost-commit] svn:boost r79911 - in trunk/boost/math/distributions: . detail
From: pbristow_at_[hidden]
Date: 2012-08-07 11:53:55

Author: pbristow
Date: 2012-08-07 11:53:54 EDT (Tue, 07 Aug 2012)
New Revision: 79911
URL: http://svn.boost.org/trac/boost/changeset/79911

Log:
Major update to allow df == +infinity.
Text files modified: 
   trunk/boost/math/distributions/detail/common_error_handling.hpp |     1                                         
   trunk/boost/math/distributions/students_t.hpp                   |   215 ++++++++++++++++++++++++++------------- 
   2 files changed, 144 insertions(+), 72 deletions(-)

Modified: trunk/boost/math/distributions/detail/common_error_handling.hpp
==============================================================================
--- trunk/boost/math/distributions/detail/common_error_handling.hpp	(original)
+++ trunk/boost/math/distributions/detail/common_error_handling.hpp	2012-08-07 11:53:54 EDT (Tue, 07 Aug 2012)
@@ -13,6 +13,7 @@
 #include <boost/math/special_functions/fpclassify.hpp>
 // using boost::math::isfinite;
 
+
 namespace boost{ namespace math{ namespace detail
 {
 

Modified: trunk/boost/math/distributions/students_t.hpp
==============================================================================
--- trunk/boost/math/distributions/students_t.hpp	(original)
+++ trunk/boost/math/distributions/students_t.hpp	2012-08-07 11:53:54 EDT (Tue, 07 Aug 2012)
@@ -1,5 +1,7 @@
 //  Copyright John Maddock 2006.
-//  Copyright Paul A. Bristow 2006.
+//  Copyright Paul A. Bristow 2006, 2012.
+//  Copyright Thomas Mang 2012.
+
 //  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)
@@ -25,6 +27,20 @@
 
 namespace boost{ namespace math{
 
+template <class RealType, class Policy>
+inline bool check_df(const char* function, RealType const& df, RealType* result, const Policy& pol)
+{  // df > 0 or +infinity are allowed.
+   // if((df <= 0) || (boost::math::isnan)(df)) // but use signbit to ensure catch -inf and -NaN.
+   if(((boost::math::signbit)(df) != 0) || (boost::math::isnan)(df))
+   { // is bad df <= 0 or NaN or -infinity.
+      *result = policies::raise_domain_error<RealType>(
+         function,
+         "Degrees of freedom argument is %1%, but must be > 0 !", df, pol);
+      return false;
+   }
+   return true;
+} // check_df
+
 template <class RealType = double, class Policy = policies::policy<> >
 class students_t_distribution
 {
@@ -32,16 +48,16 @@
    typedef RealType value_type;
    typedef Policy policy_type;
 
-   students_t_distribution(RealType i) : m_df(i)
+   students_t_distribution(RealType df) : df_(df)
    { // Constructor.
       RealType result;
-      detail::check_df(
-         "boost::math::students_t_distribution<%1%>::students_t_distribution", m_df, &result, Policy());
+      check_df( // Checks that df > 0 or df == inf.
+         "boost::math::students_t_distribution<%1%>::students_t_distribution", df_, &result, Policy());
    } // students_t_distribution
 
    RealType degrees_of_freedom()const
    {
-      return m_df;
+      return df_;
    }
 
    // Parameter estimation:
@@ -53,18 +69,16 @@
       RealType hint = 100);
 
 private:
-   //
    // Data member:
-   //
-   RealType m_df;  // degrees of freedom is a real number.
+   RealType df_;  // degrees of freedom is a real number or +infinity.
 };
 
-typedef students_t_distribution<double> students_t;
+typedef students_t_distribution<double> students_t; // Convenience typedef for double version.
 
 template <class RealType, class Policy>
 inline const std::pair<RealType, RealType> range(const students_t_distribution<RealType, Policy>& /*dist*/)
 { // Range of permissible values for random variable x.
-  // No including infinity.
+  // NOT including infinity.
    using boost::math::tools::max_value;
    return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>());
 }
@@ -81,81 +95,89 @@
 inline RealType pdf(const students_t_distribution<RealType, Policy>& dist, const RealType& x)
 {
    BOOST_FPU_EXCEPTION_GUARD
-   BOOST_MATH_STD_USING  // for ADL of std functions
+   BOOST_MATH_STD_USING  // for ADL of std functions.
 
-   RealType degrees_of_freedom = dist.degrees_of_freedom();
-   // common handling Error check is for finite and > 0:
-   // Might conceivably permit df = +infinity too?
    RealType error_result;
-   if(false == detail::check_df(
-      "boost::math::pdf(const students_t_distribution<%1%>&, %1%)", degrees_of_freedom, &error_result, Policy()))
-      return error_result;
    if(false == detail::check_x(
       "boost::math::pdf(const students_t_distribution<%1%>&, %1%)", x, &error_result, Policy()))
       return error_result;
+   RealType df = dist.degrees_of_freedom();
+   if(false == check_df( // Check that df > 0 or == +infinity.
+      "boost::math::pdf(const students_t_distribution<%1%>&, %1%)", df, &error_result, Policy()))
+      return error_result;
 
    RealType result;
-
+   if (boost::math::isinf(x))
+   { // +infinity.
+     normal_distribution<RealType, Policy> n(0, 1); 
+     RealType result = pdf(n, x);
+     return result;
+   }
    RealType limit = policies::get_epsilon<RealType, Policy>();
    // Use policies so that if policy requests lower precision, 
    // then get the normal distribution approximation earlier.
    limit = static_cast<RealType>(1) / limit; // 1/eps
    // for 64-bit double 1/eps = 4503599627370496
-   if (degrees_of_freedom > limit)
-   { // Special case for really big degrees_of_freedom > 1 / eps (perhaps infinite?)
+   if (df > limit)
+   { // Special case for really big degrees_of_freedom > 1 / eps 
      // - use normal distribution which is much faster and more accurate.
      normal_distribution<RealType, Policy> n(0, 1); 
      result = pdf(n, x);
    }
    else
    { // 
-     RealType basem1 = x * x / degrees_of_freedom;
+     RealType basem1 = x * x / df;
      if(basem1 < 0.125)
      {
-        result = exp(-boost::math::log1p(basem1, Policy()) * (1+degrees_of_freedom) / 2);
+        result = exp(-boost::math::log1p(basem1, Policy()) * (1+df) / 2);
      }
      else
      {
-        result = pow(1 / (1 + basem1), (degrees_of_freedom + 1) / 2);
+        result = pow(1 / (1 + basem1), (df + 1) / 2);
      }
-     result /= sqrt(degrees_of_freedom) * boost::math::beta(degrees_of_freedom / 2, RealType(0.5f), Policy());
+     result /= sqrt(df) * boost::math::beta(df / 2, RealType(0.5f), Policy());
    }
    return result;
 } // pdf
 
 template <class RealType, class Policy>
-inline RealType cdf(const students_t_distribution<RealType, Policy>& dist, const RealType& t)
+inline RealType cdf(const students_t_distribution<RealType, Policy>& dist, const RealType& x)
 {
-   RealType degrees_of_freedom = dist.degrees_of_freedom();
-   // Error check:
    RealType error_result;
-   if(false == detail::check_df(
-      "boost::math::cdf(const students_t_distribution<%1%>&, %1%)", degrees_of_freedom, &error_result, Policy()))
-      return error_result;
    if(false == detail::check_x(
-      "boost::math::pdf(const students_t_distribution<%1%>&, %1%)", t, &error_result, Policy()))
+      "boost::math::pdf(const students_t_distribution<%1%>&, %1%)", x, &error_result, Policy()))
       return error_result;
+   RealType df = dist.degrees_of_freedom();
+   // Error check:
 
-   if (t == 0)
-   {
-     return 0.5;
-   }
+   if(false == check_df(  // Check that df > 0 or == +infinity.
+      "boost::math::cdf(const students_t_distribution<%1%>&, %1%)", df, &error_result, Policy()))
+      return error_result;
 
+   if (x == 0)
+   { // Special case with exact result.
+     return static_cast<RealType>(0.5);
+   }
+   if (boost::math::isinf(x))
+   { // +infinity.
+     normal_distribution<RealType, Policy> n(0, 1); 
+     RealType result = cdf(n, x);
+     return result;
+   }
    RealType limit = policies::get_epsilon<RealType, Policy>();
    // Use policies so that if policy requests lower precision, 
    // then get the normal distribution approximation earlier.
    limit = static_cast<RealType>(1) / limit; // 1/eps
    // for 64-bit double 1/eps = 4503599627370496
-   if (degrees_of_freedom > limit)
+   if (df > limit)
    { // Special case for really big degrees_of_freedom > 1 / eps (perhaps infinite?)
      // - use normal distribution which is much faster and more accurate.
      normal_distribution<RealType, Policy> n(0, 1); 
-     RealType result = cdf(n, t);
+     RealType result = cdf(n, x);
      return result;
    }
    else
-   { // normal case.
-
+   { // normal df case.
      //
      // Calculate probability of Student's t using the incomplete beta function.
      // probability = ibeta(degrees_of_freedom / 2, 1/2, degrees_of_freedom / (degrees_of_freedom + t*t))
@@ -174,19 +196,19 @@
      //
      // 1 - x = t^2 / (df + t^2)
      //
-     RealType t2 = t * t;
+     RealType x2 = x * x;
      RealType probability;
-     if(degrees_of_freedom > 2 * t2)
+     if(df > 2 * x2)
      {
-        RealType z = t2 / (degrees_of_freedom + t2);
-        probability = ibetac(static_cast<RealType>(0.5), degrees_of_freedom / 2, z, Policy()) / 2;
+        RealType z = x2 / (df + x2);
+        probability = ibetac(static_cast<RealType>(0.5), df / 2, z, Policy()) / 2;
      }
      else
      {
-        RealType z = degrees_of_freedom / (degrees_of_freedom + t2);
-        probability = ibeta(degrees_of_freedom / 2, static_cast<RealType>(0.5), z, Policy()) / 2;
+        RealType z = df / (df + x2);
+        probability = ibeta(df / 2, static_cast<RealType>(0.5), z, Policy()) / 2;
      }
-     return (t > 0 ? 1   - probability : probability);
+     return (x > 0 ? 1   - probability : probability);
   }
 } // cdf
 
@@ -196,26 +218,23 @@
    BOOST_MATH_STD_USING // for ADL of std functions
    //
    // Obtain parameters:
-   //
-   RealType degrees_of_freedom = dist.degrees_of_freedom();
    RealType probability = p;
-   //
+ 
    // Check for domain errors:
-   //
+   RealType df = dist.degrees_of_freedom();
    static const char* function = "boost::math::quantile(const students_t_distribution<%1%>&, %1%)";
    RealType error_result;
-   if(false == (detail::check_df(
-      function, degrees_of_freedom, &error_result, Policy())
+   if(false == (check_df( // Check that df > 0 or == +infinity.
+      function, df, &error_result, Policy())
          && detail::check_probability(function, probability, &error_result, Policy())))
       return error_result;
-
    // Special cases, regardless of degrees_of_freedom.
    if (probability == 0)
       return -policies::raise_overflow_error<RealType>(function, 0, Policy());
    if (probability == 1)
      return policies::raise_overflow_error<RealType>(function, 0, Policy());
    if (probability == static_cast<RealType>(0.5))
-     return 0;
+     return 0;  //
    //
 #if 0
    // This next block is disabled in favour of a faster method than
@@ -245,7 +264,7 @@
    // and a couple of epsilon at double precision and in the central 
    // region where most use cases will occur...
    //
-   return boost::math::detail::fast_students_t_quantile(degrees_of_freedom, probability, Policy());
+   return boost::math::detail::fast_students_t_quantile(df, probability, Policy());
 } // quantile
 
 template <class RealType, class Policy>
@@ -276,7 +295,9 @@
    RealType operator()(const RealType& df)
    {
       if(df <= tools::min_value<RealType>())
+      { // 
          return 1;
+      }
       students_t_distribution<RealType, Policy> t(df);
       RealType qa = quantile(complement(t, alpha));
       RealType qb = quantile(complement(t, beta));
@@ -329,13 +350,15 @@
 template <class RealType, class Policy>
 inline RealType mode(const students_t_distribution<RealType, Policy>& /*dist*/)
 {
-   return 0;
+  // Assume no checks on degrees of freedom are useful (unlike mean).
+   return 0; // Always zero by definition.
 }
 
 template <class RealType, class Policy>
 inline RealType median(const students_t_distribution<RealType, Policy>& /*dist*/)
 {
-   return 0;
+   // Assume no checks on degrees of freedom are useful (unlike mean).
+   return 0; // Always zero by definition.
 }
 
 // See section 5.1 on moments at  http://en.wikipedia.org/wiki/Student%27s_t-distribution
@@ -344,37 +367,53 @@
 inline RealType mean(const students_t_distribution<RealType, Policy>& dist)
 {  // Revised for https://svn.boost.org/trac/boost/ticket/7177
    RealType df = dist.degrees_of_freedom();
-   if((!(boost::math::isfinite)(df)) || (df <= 1) )// Undefined for moment <= 1
-   {
+   if(((boost::math::isnan)(df)) || (df <= 1) ) 
+   { // mean is undefined for moment <= 1!
       policies::raise_domain_error<RealType>(
       "boost::math::mean(students_t_distribution<%1%> const&, %1%)",
       "Mean is undefined for degrees of freedom < 1 but got %1%.", df, Policy());
       return std::numeric_limits<RealType>::quiet_NaN();
    }
    return 0;
-}
+} // mean
 
 template <class RealType, class Policy>
 inline RealType variance(const students_t_distribution<RealType, Policy>& dist)
 { // http://en.wikipedia.org/wiki/Student%27s_t-distribution
   // Revised for https://svn.boost.org/trac/boost/ticket/7177
   RealType df = dist.degrees_of_freedom();
-  if (!(boost::math::isfinite)(df) ||  (df <= 2))
-  { // Infinity or NaN.
+  if ((boost::math::isnan)(df) || (df <= 2))
+  { // NaN or undefined for <= 2.
      policies::raise_domain_error<RealType>(
       "boost::math::variance(students_t_distribution<%1%> const&, %1%)",
       "variance is undefined for degrees of freedom <= 2, but got %1%.",
       df, Policy());
     return std::numeric_limits<RealType>::quiet_NaN(); // Undefined.
   }
-  return df / (df - 2);
+  if (boost::math::isinf(df))
+  { // +infinity.
+    return 1;
+  }
+  RealType limit = policies::get_epsilon<RealType, Policy>();
+  // Use policies so that if policy requests lower precision, 
+  // then get the normal distribution approximation earlier.
+  limit = static_cast<RealType>(1) / limit; // 1/eps
+  // for 64-bit double 1/eps = 4503599627370496
+  if (df > limit)
+  { // Special case for really big degrees_of_freedom > 1 / eps.
+    return 1;
+  }
+  else
+  {
+    return df / (df - 2);
+  }
 } // variance
 
 template <class RealType, class Policy>
 inline RealType skewness(const students_t_distribution<RealType, Policy>& dist)
 {
     RealType df = dist.degrees_of_freedom();
-   if( (!(boost::math::isfinite)(df)) || (dist.degrees_of_freedom() <= 3))
+   if( ((boost::math::isnan)(df)) || (dist.degrees_of_freedom() <= 3))
    { // Undefined for moment k = 3.
       policies::raise_domain_error<RealType>(
          "boost::math::skewness(students_t_distribution<%1%> const&, %1%)",
@@ -382,14 +421,14 @@
          dist.degrees_of_freedom(), Policy());
       return std::numeric_limits<RealType>::quiet_NaN();
    }
-   return 0;
-}
+   return 0; // For all valid df, including infinity.
+} // skewness
 
 template <class RealType, class Policy>
 inline RealType kurtosis(const students_t_distribution<RealType, Policy>& dist)
 {
    RealType df = dist.degrees_of_freedom();
-   if((!(boost::math::isfinite)(df)) || (df <= 4))
+   if(((boost::math::isnan)(df)) || (df <= 4))
    { // Undefined or infinity for moment k = 4.
       policies::raise_domain_error<RealType>(
        "boost::math::kurtosis(students_t_distribution<%1%> const&, %1%)",
@@ -397,9 +436,25 @@
         df, Policy());
         return std::numeric_limits<RealType>::quiet_NaN(); // Undefined.
    }
-   //return 3 * (df - 2) / (df - 4);
-   return 6 / (df - 4) + 3;
-}
+   if (boost::math::isinf(df))
+   { // +infinity.
+     return 3;
+   }
+   RealType limit = policies::get_epsilon<RealType, Policy>();
+   // Use policies so that if policy requests lower precision, 
+   // then get the normal distribution approximation earlier.
+   limit = static_cast<RealType>(1) / limit; // 1/eps
+   // for 64-bit double 1/eps = 4503599627370496
+   if (df > limit)
+   { // Special case for really big degrees_of_freedom > 1 / eps.
+     return 3;
+   }
+   else
+   {
+     //return 3 * (df - 2) / (df - 4); re-arranged to
+     return 6 / (df - 4) + 3;
+   }
+} // kurtosis
 
 template <class RealType, class Policy>
 inline RealType kurtosis_excess(const students_t_distribution<RealType, Policy>& dist)
@@ -407,7 +462,7 @@
    // see http://mathworld.wolfram.com/Kurtosis.html
 
    RealType df = dist.degrees_of_freedom();
-   if((!(boost::math::isfinite)(df)) || (df <= 4))
+   if(((boost::math::isnan)(df)) || (df <= 4))
    { // Undefined or infinity for moment k = 4.
      policies::raise_domain_error<RealType>(
        "boost::math::kurtosis_excess(students_t_distribution<%1%> const&, %1%)",
@@ -415,7 +470,23 @@
       df, Policy());
      return std::numeric_limits<RealType>::quiet_NaN(); // Undefined.
    }
-   return 6 / (df - 4);
+   if (boost::math::isinf(df))
+   { // +infinity.
+     return 0;
+   }
+   RealType limit = policies::get_epsilon<RealType, Policy>();
+   // Use policies so that if policy requests lower precision, 
+   // then get the normal distribution approximation earlier.
+   limit = static_cast<RealType>(1) / limit; // 1/eps
+   // for 64-bit double 1/eps = 4503599627370496
+   if (df > limit)
+   { // Special case for really big degrees_of_freedom > 1 / eps.
+     return 0;
+   }
+   else
+   {
+     return 6 / (df - 4);
+   }
 }
 
 } // namespace math