boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-29 07:42:11 +00:00
Code Issues Projects Releases Wiki Activity

Major update to allow df == +infinity.

Browse Source 
[SVN r79911]
This commit is contained in:
Paul A. Bristow
2012-08-07 15:53:54 +00:00
parent 5d887432cd
commit cf1d644c4e
2 changed files with 144 additions and 72 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
include/boost/math/distributions/detail/common_error_handling.hpp
View File

@@ -13,6 +13,7 @@
#include <boost/math/special_functions/fpclassify.hpp>
// using boost::math::isfinite;
namespace boost{ namespace math{ namespace detail
{

215
include/boost/math/distributions/students_t.hpp
View File

@@ -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 @@ public:
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 @@ public:
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 @@ template <class RealType, class Policy>
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(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 (t == 0)
{
return 0.5;
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 @@ inline RealType cdf(const students_t_distribution<RealType, Policy>& dist, const
//
// 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 @@ inline RealType quantile(const students_t_distribution<RealType, Policy>& dist,
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 @@ inline RealType quantile(const students_t_distribution<RealType, Policy>& dist,
// 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 @@ struct sample_size_func
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 @@ RealType students_t_distribution<RealType, Policy>::find_degrees_of_freedom(
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 @@ template <class RealType, class Policy>
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 @@ inline RealType skewness(const students_t_distribution<RealType, Policy>& dist)
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 @@ inline RealType kurtosis(const students_t_distribution<RealType, Policy>& dist)
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 @@ inline RealType kurtosis_excess(const students_t_distribution<RealType, Policy>&
// 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 @@ inline RealType kurtosis_excess(const students_t_distribution<RealType, Policy>&
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