diff --git a/include/boost/math/distributions/detail/common_error_handling.hpp b/include/boost/math/distributions/detail/common_error_handling.hpp index 46ded33cd..931747485 100644 --- a/include/boost/math/distributions/detail/common_error_handling.hpp +++ b/include/boost/math/distributions/detail/common_error_handling.hpp @@ -13,6 +13,7 @@ #include // using boost::math::isfinite; + namespace boost{ namespace math{ namespace detail { diff --git a/include/boost/math/distributions/students_t.hpp b/include/boost/math/distributions/students_t.hpp index 5dfce88c9..5ed4c9877 100644 --- a/include/boost/math/distributions/students_t.hpp +++ b/include/boost/math/distributions/students_t.hpp @@ -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 +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( + function, + "Degrees of freedom argument is %1%, but must be > 0 !", df, pol); + return false; + } + return true; +} // check_df + template > 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 students_t; +typedef students_t_distribution students_t; // Convenience typedef for double version. template inline const std::pair range(const students_t_distribution& /*dist*/) { // Range of permissible values for random variable x. - // No including infinity. + // NOT including infinity. using boost::math::tools::max_value; return std::pair(-max_value(), max_value()); } @@ -81,81 +95,89 @@ template inline RealType pdf(const students_t_distribution& 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 n(0, 1); + RealType result = pdf(n, x); + return result; + } RealType limit = policies::get_epsilon(); // Use policies so that if policy requests lower precision, // then get the normal distribution approximation earlier. limit = static_cast(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 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 -inline RealType cdf(const students_t_distribution& dist, const RealType& t) +inline RealType cdf(const students_t_distribution& 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(0.5); + } + if (boost::math::isinf(x)) + { // +infinity. + normal_distribution n(0, 1); + RealType result = cdf(n, x); + return result; } - RealType limit = policies::get_epsilon(); // Use policies so that if policy requests lower precision, // then get the normal distribution approximation earlier. limit = static_cast(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 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& 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(0.5), degrees_of_freedom / 2, z, Policy()) / 2; + RealType z = x2 / (df + x2); + probability = ibetac(static_cast(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(0.5), z, Policy()) / 2; + RealType z = df / (df + x2); + probability = ibeta(df / 2, static_cast(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& 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(function, 0, Policy()); if (probability == 1) return policies::raise_overflow_error(function, 0, Policy()); if (probability == static_cast(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& 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 @@ -276,7 +295,9 @@ struct sample_size_func RealType operator()(const RealType& df) { if(df <= tools::min_value()) + { // return 1; + } students_t_distribution t(df); RealType qa = quantile(complement(t, alpha)); RealType qb = quantile(complement(t, beta)); @@ -329,13 +350,15 @@ RealType students_t_distribution::find_degrees_of_freedom( template inline RealType mode(const students_t_distribution& /*dist*/) { - return 0; + // Assume no checks on degrees of freedom are useful (unlike mean). + return 0; // Always zero by definition. } template inline RealType median(const students_t_distribution& /*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 inline RealType mean(const students_t_distribution& 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( "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::quiet_NaN(); } return 0; -} +} // mean template inline RealType variance(const students_t_distribution& 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( "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::quiet_NaN(); // Undefined. } - return df / (df - 2); + if (boost::math::isinf(df)) + { // +infinity. + return 1; + } + RealType limit = policies::get_epsilon(); + // Use policies so that if policy requests lower precision, + // then get the normal distribution approximation earlier. + limit = static_cast(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 inline RealType skewness(const students_t_distribution& 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( "boost::math::skewness(students_t_distribution<%1%> const&, %1%)", @@ -382,14 +421,14 @@ inline RealType skewness(const students_t_distribution& dist) dist.degrees_of_freedom(), Policy()); return std::numeric_limits::quiet_NaN(); } - return 0; -} + return 0; // For all valid df, including infinity. +} // skewness template inline RealType kurtosis(const students_t_distribution& 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( "boost::math::kurtosis(students_t_distribution<%1%> const&, %1%)", @@ -397,9 +436,25 @@ inline RealType kurtosis(const students_t_distribution& dist) df, Policy()); return std::numeric_limits::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(); + // Use policies so that if policy requests lower precision, + // then get the normal distribution approximation earlier. + limit = static_cast(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 inline RealType kurtosis_excess(const students_t_distribution& dist) @@ -407,7 +462,7 @@ inline RealType kurtosis_excess(const students_t_distribution& // 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( "boost::math::kurtosis_excess(students_t_distribution<%1%> const&, %1%)", @@ -415,7 +470,23 @@ inline RealType kurtosis_excess(const students_t_distribution& df, Policy()); return std::numeric_limits::quiet_NaN(); // Undefined. } - return 6 / (df - 4); + if (boost::math::isinf(df)) + { // +infinity. + return 0; + } + RealType limit = policies::get_epsilon(); + // Use policies so that if policy requests lower precision, + // then get the normal distribution approximation earlier. + limit = static_cast(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