|
|
Boost-Commit : |
From: john_at_[hidden]
Date: 2008-02-04 13:57:38
Author: johnmaddock
Date: 2008-02-04 13:57:37 EST (Mon, 04 Feb 2008)
New Revision: 43098
URL: http://svn.boost.org/trac/boost/changeset/43098
Log:
Made the non-central chi-squared PDF more robust.
Set generic mode calculation to throw if the original guess yields a zero PDF.
Updated tests to match.
Text files modified:
sandbox/math_toolkit/boost/math/distributions/detail/generic_mode.hpp | 10 +++++
sandbox/math_toolkit/boost/math/distributions/non_central_chi_squared.hpp | 68 ++++++++++++++++++++++++++++++++++-----
sandbox/math_toolkit/libs/math/test/test_nc_chi_squared.cpp | 15 ++-----
3 files changed, 73 insertions(+), 20 deletions(-)
Modified: sandbox/math_toolkit/boost/math/distributions/detail/generic_mode.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/detail/generic_mode.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions/detail/generic_mode.hpp 2008-02-04 13:57:37 EST (Mon, 04 Feb 2008)
@@ -41,6 +41,16 @@
value_type upper_bound = guess;
value_type lower_bound;
value_type v = pdf(dist, guess);
+ if(v == 0)
+ {
+ //
+ // Oops we don't know how to handle this, or even in which
+ // direction we should move in, treat as an evaluation error:
+ //
+ policies::raise_evaluation_error(
+ function,
+ "Could not locate a starting location for the search for the mode, original guess was %1%", guess, policy_type());
+ }
do
{
maxval = v;
Modified: sandbox/math_toolkit/boost/math/distributions/non_central_chi_squared.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/non_central_chi_squared.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions/non_central_chi_squared.hpp 2008-02-04 13:57:37 EST (Mon, 04 Feb 2008)
@@ -280,6 +280,45 @@
return sum;
}
+ template <class T, class Policy>
+ T non_central_chi_square_pdf(T x, T n, T lambda, const Policy& pol)
+ {
+ //
+ // As above but for the PDF:
+ //
+ BOOST_MATH_STD_USING
+ boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
+ T errtol = ldexp(1.0, -boost::math::policies::digits<T, Policy>());
+ T x2 = x / 2;
+ T n2 = n / 2;
+ T l2 = lambda / 2;
+ T sum = 0;
+ int k = itrunc(l2);
+ T pois = gamma_p_derivative(k + 1, l2, pol) * gamma_p_derivative(n2 + k, x2);
+ if(pois == 0)
+ return 0;
+ T poisb = pois;
+ for(int i = k; ; ++i)
+ {
+ sum += pois;
+ if(pois / sum < errtol)
+ break;
+ if(static_cast<boost::uintmax_t>(i - k) >= max_iter)
+ return policies::raise_evaluation_error(
+ "pdf(non_central_chi_squared_distribution<%1%>, %1%)",
+ "Series did not converge, closest value was %1%", sum, pol);
+ pois *= l2 * x2 / ((i + 1) * (n2 + i));
+ }
+ for(int i = k - 1; i >= 0; --i)
+ {
+ poisb *= (i + 1) * (n2 + i) / (l2 * x2);
+ sum += poisb;
+ if(poisb / sum < errtol)
+ break;
+ }
+ return sum / 2;
+ }
+
template <class RealType, class Policy>
inline RealType non_central_chi_squared_cdf(RealType x, RealType k, RealType l, bool invert, const Policy&)
{
@@ -470,16 +509,25 @@
// Special case:
if(x == 0)
return 0;
- r = log(x / l) * (k / 4 - 0.5f) - (x + l) / 2;
- if(r >= tools::log_max_value<RealType>())
- return policies::raise_overflow_error<RealType>(function, 0, forwarding_policy());
- if(r <= -tools::log_max_value<RealType>())
- return policies::raise_underflow_error<RealType>(function, 0, forwarding_policy());
-
- r = exp(r);
- r = 0.5f * r
- * boost::math::cyl_bessel_i(k/2 - 1, sqrt(l * x), forwarding_policy());
- return policies::checked_narrowing_cast<RealType, forwarding_policy>(
+ if(l > 50)
+ {
+ r = non_central_chi_square_pdf(static_cast<value_type>(x), k, l, forwarding_policy());
+ }
+ else
+ {
+ r = log(x / l) * (k / 4 - 0.5f) - (x + l) / 2;
+ if(fabs(r) >= tools::log_max_value<RealType>() / 4)
+ {
+ r = non_central_chi_square_pdf(static_cast<value_type>(x), k, l, forwarding_policy());
+ }
+ else
+ {
+ r = exp(r);
+ r = 0.5f * r
+ * boost::math::cyl_bessel_i(k/2 - 1, sqrt(l * x), forwarding_policy());
+ }
+ }
+ return policies::checked_narrowing_cast<RealType, forwarding_policy>(
r,
function);
}
Modified: sandbox/math_toolkit/libs/math/test/test_nc_chi_squared.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/test/test_nc_chi_squared.cpp (original)
+++ sandbox/math_toolkit/libs/math/test/test_nc_chi_squared.cpp 2008-02-04 13:57:37 EST (Mon, 04 Feb 2008)
@@ -441,21 +441,16 @@
if(boost::math::tools::digits<value_type>() > 50)
{
//
- // Sanity check mode, note this may well overflow
- // since we don't know how to compute PDF's (and hence the mode)
- // for large values of the parameters, in addition accuracy of
+ // Sanity check mode, the accuracy of
// the mode is at *best* the square root of the accuracy of the PDF:
//
- try
- {
+ try{
value_type m = mode(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]));
value_type p = pdf(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), m);
- BOOST_CHECK_EX(pdf(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), m * (1 + sqrt(precision) * 10)) <= p, i);
- BOOST_CHECK_EX(pdf(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), m * (1 - sqrt(precision)) * 10) <= p, i);
- }
- catch(const std::overflow_error&)
- {
+ BOOST_CHECK_EX(pdf(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), m * (1 + sqrt(precision) * 50)) <= p, i);
+ BOOST_CHECK_EX(pdf(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), m * (1 - sqrt(precision)) * 50) <= p, i);
}
+ catch(const boost::math::evaluation_error& ) {}
//
// Sanity check degrees-of-freedom finder, don't bother at float
// precision though as there's not enough data in the probability
Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk