From: john_at_[hidden]
Date: 2008-01-30 07:42:27
Author: johnmaddock
Date: 2008-01-30 07:42:26 EST (Wed, 30 Jan 2008)
New Revision: 43020
URL: http://svn.boost.org/trac/boost/changeset/43020
Log:
More or less finished off the non central beta.
Text files modified:
sandbox/math_toolkit/boost/math/distributions/non_central_beta.hpp | 136 ++++++++++++++++++++++++--------
sandbox/math_toolkit/libs/math/test/Jamfile.v2 | 25 +++++
sandbox/math_toolkit/libs/math/test/test_nc_beta.cpp | 169 +++++++++++++++++++++++----------------
3 files changed, 227 insertions(+), 103 deletions(-)
Modified: sandbox/math_toolkit/boost/math/distributions/non_central_beta.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/non_central_beta.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions/non_central_beta.hpp 2008-01-30 07:42:26 EST (Wed, 30 Jan 2008)
@@ -14,11 +14,11 @@
#include <boost/math/special_functions/beta.hpp> // for incomplete gamma. gamma_q
#include <boost/math/distributions/complement.hpp> // complements
#include <boost/math/distributions/beta.hpp> // central distribution
-#include <boost/math/distributions/normal.hpp> // central distribution
+//#include <boost/math/distributions/normal.hpp> // central distribution
#include <boost/math/distributions/detail/common_error_handling.hpp> // error checks
#include <boost/math/special_functions/fpclassify.hpp> // isnan.
#include <boost/math/tools/roots.hpp> // for root finding.
-#include <boost/math/tools/minima.hpp> // function minimization for mode
+//#include <boost/math/tools/minima.hpp> // function minimization for mode
namespace boost
{
@@ -73,7 +73,7 @@
xterm *= (a + i - 1) / (x * (a + b + i - 2));
last_term = term;
}
- for(int i = k + 1; i - k < max_iter; ++i)
+ for(int i = k + 1; static_cast<boost::uintmax_t>(i - k) < max_iter; ++i)
{
poisf *= l2 / i;
xtermf *= (x * (a + b + i - 2)) / (a + i - 1);
@@ -113,25 +113,13 @@
T xterm = ibeta_derivative(a + k, b, x, pol) * (1 - x) / (a + b + k - 1);
T poisf(pois), betaf(beta), xtermf(xterm);
T sum = init_val;
-
//
- // Backwards recursion first, this is the unstable
- // direction for recursion:
+ // Forwards recursion first, this is the stable
+ // direction for recursion, and the location
+ // of the bulk of the sum:
//
- for(int i = k; i >= 0; --i)
- {
- T term = beta * pois;
- sum += term;
- if(fabs(term/sum) < errtol)
- {
- break;
- }
- pois *= i / l2;
- beta -= xterm;
- xterm *= (a + i - 1) / (x * (a + b + i - 2));
- }
T last_term = 0;
- for(int i = k + 1; i - k < max_iter; ++i)
+ for(int i = k + 1; static_cast<boost::uintmax_t>(i - k) < max_iter; ++i)
{
poisf *= l2 / i;
xtermf *= (x * (a + b + i - 2)) / (a + i - 1);
@@ -145,6 +133,18 @@
}
last_term = term;
}
+ for(int i = k; i >= 0; --i)
+ {
+ T term = beta * pois;
+ sum += term;
+ if(fabs(term/sum) < errtol)
+ {
+ break;
+ }
+ pois *= i / l2;
+ beta -= xterm;
+ xterm *= (a + i - 1) / (x * (a + b + i - 2));
+ }
return sum;
}
@@ -162,9 +162,9 @@
BOOST_MATH_STD_USING
if(x == 0)
- return invert ? 1 : 0;
+ return invert ? 1.0f : 0.0f;
if(x == 1)
- return invert ? 0 : 1;
+ return invert ? 0.0f : 1.0f;
value_type result;
value_type c = a + b + l / 2;
value_type cross = 1 - (b / c) * (1 + l / (2 * c * c));
@@ -218,6 +218,11 @@
bool comp;
};
+ //
+ // This is more or less a copy of bracket_and_solve_root, but
+ // modified to search only the interval [0,1] using similar
+ // heuristics.
+ //
template <class F, class T, class Tol, class Policy>
std::pair<T, T> bracket_and_solve_root_01(F f, const T& guess, T factor, bool rising, Tol tol, boost::uintmax_t& max_iter, const Policy& pol)
{
@@ -252,7 +257,8 @@
if((max_iter - count) % 20 == 0)
factor *= 2;
//
- // Now go ahead and move are guess by "factor":
+ // Now go ahead and move are guess by "factor",
+ // we do this by reducing 1-guess by factor:
//
a = b;
fa = fb;
@@ -353,12 +359,12 @@
//
if(p == 0)
return comp
- ? 1
- : 0;
+ ? 1.0f
+ : 0.0f;
if(p == 1)
return !comp
- ? 1
- : 0;
+ ? 1.0f
+ : 0.0f;
value_type c = a + b + l / 2;
value_type mean = 1 - (b / c) * (1 + l / (2 * c * c));
@@ -431,6 +437,59 @@
function);
}
+ template <class T, class Policy>
+ T non_central_beta_pdf(T a, T b, T lam, T x, const Policy& pol)
+ {
+ BOOST_MATH_STD_USING
+ using namespace boost::math;
+ //
+ // Variables come first:
+ //
+ boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
+ T errtol = ldexp(1.0, -boost::math::policies::digits<T, Policy>());
+ T l2 = lam / 2;
+ //
+ // k is the starting point for iteration, and is the
+ // maximum of the poisson weighting term:
+ //
+ int k = itrunc(l2);
+ // Starting Poisson weight:
+ T pois = gamma_p_derivative(k+1, l2, pol);
+ // Starting beta term:
+ T beta = ibeta_derivative(a + k, b, x, pol);
+ T sum = 0;
+ T poisf(pois);
+ T betaf(beta);
+
+ //
+ // Stable backwards recursion first:
+ //
+ for(int i = k; i >= 0; --i)
+ {
+ T term = beta * pois;
+ sum += term;
+ if(fabs(term/sum) < errtol)
+ {
+ break;
+ }
+ pois *= i / l2;
+ beta *= (a + i - 1) / (x * (a + i + b - 1));
+ }
+ for(int i = k + 1; static_cast<boost::uintmax_t>(i - k) < max_iter; ++i)
+ {
+ poisf *= l2 / i;
+ betaf *= x * (a + b + i - 1) / (a + i - 1);
+
+ T term = poisf * betaf;
+ sum += term;
+ if(fabs(term/sum) < errtol)
+ {
+ break;
+ }
+ }
+ return sum;
+ }
+
template <class RealType, class Policy>
RealType nc_beta_pdf(const non_central_beta_distribution<RealType, Policy>& dist, const RealType& x)
{
@@ -462,20 +521,22 @@
&r,
Policy())
||
- !detail::check_probability(
+ !beta_detail::check_x(
function,
- static_cast<value_type>(p),
+ static_cast<value_type>(x),
&r,
Policy()))
return (RealType)r;
- BOOST_MATH_STD_USING
- if(l == 0)
- return pdf(boost::math::beta_distribution<RealType, forwarding_policy>(dist.degrees_of_freedom()), x);
-
+ BOOST_MATH_STD_USING
+ if(l == 0)
+ return pdf(boost::math::beta_distribution<RealType, Policy>(dist.alpha(), dist.beta()), x);
+ return policies::checked_narrowing_cast<RealType, forwarding_policy>(
+ non_central_beta_pdf(a, b, l, static_cast<value_type>(x), forwarding_policy()),
+ "function");
}
- }
+ } // namespace detail
template <class RealType = double, class Policy = policies::policy<> >
class non_central_beta_distribution
@@ -539,6 +600,13 @@
return std::pair<RealType, RealType>(0, 1);
}
+#if 0
+ //
+ // We don't have the necessary information to implement
+ // these at present. These are just disabled for now,
+ // prototypes retained so we can fill in the blanks
+ // later:
+ //
template <class RealType, class Policy>
inline RealType mean(const non_central_beta_distribution<RealType, Policy>& dist)
{
@@ -585,7 +653,7 @@
{
return kurtosis_excess(dist) + 3;
}
-
+#endif
template <class RealType, class Policy>
inline RealType pdf(const non_central_beta_distribution<RealType, Policy>& dist, const RealType& x)
{ // Probability Density/Mass Function.
Modified: sandbox/math_toolkit/libs/math/test/Jamfile.v2
==============================================================================
--- sandbox/math_toolkit/libs/math/test/Jamfile.v2 (original)
+++ sandbox/math_toolkit/libs/math/test/Jamfile.v2 2008-01-30 07:42:26 EST (Wed, 30 Jan 2008)
@@ -271,6 +271,30 @@
: # requirements
<define>TEST_REAL_CONCEPT
: test_nc_chi_squared_real_concept ;
+run test_nc_beta.cpp
+ : # command line
+ : # input files
+ : # requirements
+ <define>TEST_FLOAT
+ : test_nc_beta_float ;
+run test_nc_beta.cpp
+ : # command line
+ : # input files
+ : # requirements
+ <define>TEST_DOUBLE
+ : test_nc_beta_double ;
+run test_nc_beta.cpp
+ : # command line
+ : # input files
+ : # requirements
+ <define>TEST_LDOUBLE
+ : test_nc_beta_long_double ;
+run test_nc_beta.cpp
+ : # command line
+ : # input files
+ : # requirements
+ <define>TEST_REAL_CONCEPT
+ : test_nc_beta_real_concept ;
run test_normal.cpp ;
run test_pareto.cpp ;
run test_poisson.cpp
@@ -422,3 +446,4 @@
+
Modified: sandbox/math_toolkit/libs/math/test/test_nc_beta.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/test/test_nc_beta.cpp (original)
+++ sandbox/math_toolkit/libs/math/test/test_nc_beta.cpp 2008-01-30 07:42:26 EST (Wed, 30 Jan 2008)
@@ -20,6 +20,7 @@
#include <boost/math/concepts/real_concept.hpp> // for real_concept
#include <boost/math/distributions/non_central_beta.hpp> // for chi_squared_distribution
+#include <boost/math/distributions/poisson.hpp> // for chi_squared_distribution
#include <boost/test/included/test_exec_monitor.hpp> // for test_main
#include <boost/test/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE
@@ -78,6 +79,20 @@
largest_type = "(long\\s+)?double|real_concept";
#endif
+ if(boost::math::tools::digits<long double>() == 64)
+ {
+ //
+ // Allow a small amount of error leakage from long double to double:
+ //
+ add_expected_result(
+ "[^|]*", // compiler
+ "[^|]*", // stdlib
+ "[^|]*", // platform
+ "double", // test type(s)
+ "[^|]*large[^|]*", // test data group
+ "[^|]*", 5, 5); // test function
+ }
+
//
// Catch all cases come last:
//
@@ -92,6 +107,13 @@
"[^|]*", // compiler
"[^|]*", // stdlib
"[^|]*", // platform
+ "real_concept", // test type(s)
+ "[^|]*large[^|]*", // test data group
+ "[^|]*", 20000, 4000); // test function
+ add_expected_result(
+ "[^|]*", // compiler
+ "[^|]*", // stdlib
+ "[^|]*", // platform
largest_type, // test type(s)
"[^|]*large[^|]*", // test data group
"[^|]*", 9000, 2000); // test function
@@ -104,10 +126,23 @@
}
template <class RealType>
-RealType naive_pdf(RealType v, RealType lam, RealType x)
+RealType naive_pdf(RealType a, RealType b, RealType lam, RealType x)
{
- // TODO!
- return 0;
+ using namespace boost::math;
+
+ RealType term = pdf(poisson_distribution<RealType>(lam/2), 0)
+ * ibeta_derivative(a, b, x);
+ RealType sum = term;
+
+ int i = 1;
+ while(term / sum > tools::epsilon<RealType>())
+ {
+ term = pdf(poisson_distribution<RealType>(lam/2), i)
+ * ibeta_derivative(a + i, b, x);
+ ++i;
+ sum += term;
+ }
+ return sum;
}
template <class RealType>
@@ -118,17 +153,24 @@
RealType cs, // Chi Square statistic
RealType P, // CDF
RealType Q, // Complement of CDF
+ RealType D, // PDF
RealType tol) // Test tolerance
{
boost::math::non_central_beta_distribution<RealType> dist(a, b, ncp);
BOOST_CHECK_CLOSE(
cdf(dist, cs), P, tol);
- try{/*
+ //
+ // Sanity checking using the naive PDF calculation above fails at
+ // float precision:
+ //
+ if(!boost::is_same<float, RealType>::value)
+ {
BOOST_CHECK_CLOSE(
- pdf(dist, cs), naive_pdf(dist.degrees_of_freedom(), ncp, cs), tol * 50);*/
+ pdf(dist, cs), naive_pdf(dist.alpha(), dist.beta(), ncp, cs), tol);
}
- catch(const std::overflow_error&)
- {}
+ BOOST_CHECK_CLOSE(
+ pdf(dist, cs), D, tol);
+
if((P < 0.99) && (Q < 0.99))
{
//
@@ -148,20 +190,52 @@
void test_spots(RealType)
{
RealType tolerance = (std::max)(
- boost::math::tools::epsilon<RealType>(),
- (RealType)boost::math::tools::epsilon<double>() * 5) * 100;
- //
- // At float precision we need to up the tolerance, since
- // the input values are rounded off to inexact quantities
- // the results get thrown off by a noticeable amount.
- //
- if(boost::math::tools::digits<RealType>() < 50)
- tolerance *= 50;
- if(boost::is_floating_point<RealType>::value != 1)
- tolerance *= 20; // real_concept special functions are less accurate
+ boost::math::tools::epsilon<RealType>() * 100,
+ (RealType)1e-6) * 100;
cout << "Tolerance = " << tolerance << "%." << endl;
+ //
+ // Spot tests use values computed by the R statistical
+ // package and the pbeta and dbeta functions:
+ //
+ test_spot(
+ RealType(2), // alpha
+ RealType(5), // beta
+ RealType(1), // non-centrality param
+ RealType(0.25), // Chi Square statistic
+ RealType(0.3658349), // CDF
+ RealType(1-0.3658349), // Complement of CDF
+ RealType(2.184465), // PDF
+ RealType(tolerance));
+ test_spot(
+ RealType(20), // alpha
+ RealType(15), // beta
+ RealType(35), // non-centrality param
+ RealType(0.75), // Chi Square statistic
+ RealType(0.6994175), // CDF
+ RealType(1-0.6994175), // Complement of CDF
+ RealType(5.576146), // PDF
+ RealType(tolerance));
+ test_spot(
+ RealType(100), // alpha
+ RealType(3), // beta
+ RealType(63), // non-centrality param
+ RealType(0.95), // Chi Square statistic
+ RealType(0.03529306), // CDF
+ RealType(1-0.03529306), // Complement of CDF
+ RealType(3.637894), // PDF
+ RealType(tolerance));
+ test_spot(
+ RealType(0.25), // alpha
+ RealType(0.75), // beta
+ RealType(150), // non-centrality param
+ RealType(0.975), // Chi Square statistic
+ RealType(0.09752216), // CDF
+ RealType(1-0.09752216), // Complement of CDF
+ RealType(8.020935), // PDF
+ RealType(tolerance));
+
} // template <class RealType>void test_spots(RealType)
template <class T>
@@ -238,6 +312,13 @@
for(unsigned i = 0; i < data.size(); ++i)
{
+ //
+ // Test case 493 fails at float precision: not enough bits to get
+ // us back where we started:
+ //
+ if((i == 493) && boost::is_same<float, value_type>::value)
+ continue;
+
if(data[i][4] == 0)
{
BOOST_CHECK(0 == quantile(boost::math::non_central_beta_distribution<value_type>(data[i][0], data[i][1], data[i][2]), data[i][4]));
@@ -258,57 +339,6 @@
value_type pt = data[i][3];
BOOST_CHECK_CLOSE_EX(pt, p, precision, i);
}
- /*
- 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
- // the mode is at *best* the square root of the accuracy of the PDF:
- //
- try
- {
- value_type m = mode(boost::math::non_central_beta_distribution<value_type>(data[i][0], data[i][1]));
- value_type p = pdf(boost::math::non_central_beta_distribution<value_type>(data[i][0], data[i][1]), m);
- BOOST_CHECK_EX(pdf(boost::math::non_central_beta_distribution<value_type>(data[i][0], data[i][1]), m * (1 + sqrt(precision) * 10)) <= p, i);
- BOOST_CHECK_EX(pdf(boost::math::non_central_beta_distribution<value_type>(data[i][0], data[i][1]), m * (1 - sqrt(precision)) * 10) <= p, i);
- }
- catch(const std::overflow_error&)
- {
- }
- //
- // Sanity check degrees-of-freedom finder, don't bother at float
- // precision though as there's not enough data in the probability
- // values to get back to the correct degrees of freedom or
- // non-cenrality parameter:
- //
- try{
- if((data[i][3] < 0.99) && (data[i][3] != 0))
- {
- BOOST_CHECK_CLOSE_EX(
- boost::math::non_central_beta_distribution<value_type>::find_degrees_of_freedom(data[i][1], data[i][2], data[i][3]),
- data[i][0], precision, i);
- BOOST_CHECK_CLOSE_EX(
- boost::math::non_central_beta_distribution<value_type>::find_non_centrality(data[i][0], data[i][2], data[i][3]),
- data[i][1], precision, i);
- }
- if((data[i][4] < 0.99) && (data[i][4] != 0))
- {
- BOOST_CHECK_CLOSE_EX(
- boost::math::non_central_beta_distribution<value_type>::find_degrees_of_freedom(boost::math::complement(data[i][1], data[i][2], data[i][4])),
- data[i][0], precision, i);
- BOOST_CHECK_CLOSE_EX(
- boost::math::non_central_beta_distribution<value_type>::find_non_centrality(boost::math::complement(data[i][0], data[i][2], data[i][4])),
- data[i][1], precision, i);
- }
- }
- catch(const std::exception& e)
- {
- BOOST_ERROR(e.what());
- }
- }
- */
}
}
@@ -321,6 +351,7 @@
#include "ncbeta_big.ipp"
do_test_nc_chi_squared(ncbeta_big, type_name, "Non Central Beta, large parameters");
+ // Takes too long to run:
// quantile_sanity_check(ncbeta_big, type_name, "Non Central Beta, large parameters");
}