|
|
Boost-Commit : |
From: john_at_[hidden]
Date: 2008-04-04 11:56:11
Author: johnmaddock
Date: 2008-04-04 11:56:10 EDT (Fri, 04 Apr 2008)
New Revision: 44035
URL: http://svn.boost.org/trac/boost/changeset/44035
Log:
Fix a few more interval comparisons.
Get rounded_arith_opp finally passing the tests.
Added non-central data generators.
Added:
sandbox/interval_math_toolkit/libs/math/test/mpfr_interval_concept_check.cpp (contents, props changed)
sandbox/interval_math_toolkit/libs/math/tools/non_central_beta_data.cpp (contents, props changed)
sandbox/interval_math_toolkit/libs/math/tools/non_central_chi_sq_data.cpp (contents, props changed)
sandbox/interval_math_toolkit/libs/math/tools/non_central_t_data.cpp (contents, props changed)
Text files modified:
sandbox/interval_math_toolkit/boost/math/special_functions/beta.hpp | 2 +-
sandbox/interval_math_toolkit/boost/math/tools/config.hpp | 3 ++-
sandbox/interval_math_toolkit/boost/math/tools/test_data.hpp | 27 ++++++++++++++++++++++++---
sandbox/interval_math_toolkit/boost/numeric/interval/rounded_arith.hpp | 21 +++++----------------
4 files changed, 32 insertions(+), 21 deletions(-)
Modified: sandbox/interval_math_toolkit/boost/math/special_functions/beta.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/special_functions/beta.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/special_functions/beta.hpp 2008-04-04 11:56:10 EDT (Fri, 04 Apr 2008)
@@ -1033,7 +1033,7 @@
{
lambda = (a + b) * y - b;
}
- if(lambda < 0)
+ if(tools::maybe_less(lambda, T(0)))
{
std::swap(a, b);
std::swap(x, y);
Modified: sandbox/interval_math_toolkit/boost/math/tools/config.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/tools/config.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/tools/config.hpp 2008-04-04 11:56:10 EDT (Fri, 04 Apr 2008)
@@ -211,7 +211,8 @@
template <class T, class Policy>
std::ostream& operator << (std::ostream& os, const interval<T, Policy>& val)
{
- os << "{ " << val.lower() << ", " << val.upper() << " }";
+ T eps = width(val) / median(val) / boost::math::tools::epsilon<T>();
+ os << "{ " << val.lower() << ", " << val.upper() << " (=" << eps << "eps) }";
return os;
}
#endif
Modified: sandbox/interval_math_toolkit/boost/math/tools/test_data.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/math/tools/test_data.hpp (original)
+++ sandbox/interval_math_toolkit/boost/math/tools/test_data.hpp 2008-04-04 11:56:10 EDT (Fri, 04 Apr 2008)
@@ -200,7 +200,14 @@
boost::math::tools::detail::unpack_and_append(row, func(*a));
m_data.insert(row);
}
- catch(const std::domain_error&){}
+ catch(const std::domain_error& e)
+ {
+ std::cerr << "Ignoring domain error from test point, error was: " << e.what() << std::endl;
+ }
+ catch(const std::runtime_error& e)
+ {
+ std::cerr << "Ignoring runtime error from test point, error was: " << e.what() << std::endl;
+ }
row.clear();
++a;
}
@@ -236,7 +243,14 @@
detail::unpack_and_append(row, func(*a, *c));
m_data.insert(row);
}
- catch(const std::domain_error&){}
+ catch(const std::domain_error& e)
+ {
+ std::cerr << "Ignoring domain error from test point, error was: " << e.what() << std::endl;
+ }
+ catch(const std::runtime_error& e)
+ {
+ std::cerr << "Ignoring runtime error from test point, error was: " << e.what() << std::endl;
+ }
row.clear();
++c;
}
@@ -281,7 +295,14 @@
detail::unpack_and_append(row, func(*a, *c, *e));
m_data.insert(row);
}
- catch(const std::domain_error&){}
+ catch(const std::domain_error& e)
+ {
+ std::cerr << "Ignoring domain error from test point, error was: " << e.what() << std::endl;
+ }
+ catch(const std::runtime_error& e)
+ {
+ std::cerr << "Ignoring runtime error from test point, error was: " << e.what() << std::endl;
+ }
row.clear();
++e;
}
Modified: sandbox/interval_math_toolkit/boost/numeric/interval/rounded_arith.hpp
==============================================================================
--- sandbox/interval_math_toolkit/boost/numeric/interval/rounded_arith.hpp (original)
+++ sandbox/interval_math_toolkit/boost/numeric/interval/rounded_arith.hpp 2008-04-04 11:56:10 EDT (Fri, 04 Apr 2008)
@@ -97,23 +97,13 @@
# define BOOST_UP(EXPR) return this->force_rounding(EXPR)
# define BOOST_UP_NEG(EXPR) return -this->force_rounding(EXPR)
private:
- // Downward conversion has three cases:
- template<class U> T conv_down(U const &v, const mpl::true_&, const mpl::true_&)
- {
- // U is an integer as wide as boost::intmax_t:
- if((v > static_cast<U>(std::numeric_limits<boost::intmax_t>::max()))
- || (v < static_cast<U>(std::numeric_limits<boost::intmax_t>::min())))
- {
- BOOST_DN(static_cast<T>(v));
- }
- BOOST_UP_NEG(static_cast<T>(-static_cast<boost::intmax_t>(v)));
- }
- template<class U> T conv_down(U const &v, const mpl::true_&, const mpl::false_&)
+ // Downward conversion has two cases:
+ template<class U> T conv_down(U const &v, const mpl::true_&)
{
// U is an integer smaller than boost::intmax_t:
- BOOST_UP_NEG(static_cast<T>(-static_cast<boost::intmax_t>(v)));
+ BOOST_DN(static_cast<T>(v));
}
- template<class U> T conv_down(U const &v, const mpl::false_&, const mpl::false_&)
+ template<class U> T conv_down(U const &v, const mpl::false_&)
{
// U is not an integer, so it had better be a signed type(!):
BOOST_UP_NEG(-v);
@@ -122,8 +112,7 @@
template<class U> T conv_down(U const &v)
{
typedef typename boost::is_integral<U>::type t1;
- typedef typename mpl::and_<t1, boost::is_unsigned<U>, mpl::bool_<sizeof(U) == sizeof(boost::intmax_t)> >::type t2;
- return conv_down(v, t1(), t2());
+ return conv_down(v, t1());
}
template<class U> T conv_up (U const &v) { BOOST_UP(static_cast<T>(v)); }
T add_down(const T& x, const T& y) { BOOST_UP_NEG((-x) - y); }
Added: sandbox/interval_math_toolkit/libs/math/test/mpfr_interval_concept_check.cpp
==============================================================================
--- (empty file)
+++ sandbox/interval_math_toolkit/libs/math/test/mpfr_interval_concept_check.cpp 2008-04-04 11:56:10 EDT (Fri, 04 Apr 2008)
@@ -0,0 +1,37 @@
+// Copyright John Maddock 2007-8.
+// 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)
+
+//
+// This tests two things: that boost::interval meets our
+// conceptual requirements, and that we can instantiate
+// all our distributions and special functions on this type.
+//
+#define TEST_MPFR
+#define BOOST_MATH_ASSERT_UNDEFINED_POLICY false
+
+#pragma warning(disable:4800)
+
+#include <boost/math/bindings/interval.hpp>
+#include <boost/math/bindings/mpfr_interval.hpp>
+#include <boost/math/concepts/real_type_concept.hpp>
+#include "compile_test/instantiate.hpp"
+
+
+void foo()
+{
+ instantiate(boost::math::mpfr_interval());
+}
+
+int main()
+{
+ BOOST_CONCEPT_ASSERT((boost::math::concepts::RealTypeConcept<double>));
+ //BOOST_CONCEPT_ASSERT((boost::math::concepts::RealTypeConcept<boost::numeric::interval<double> >));
+ BOOST_CONCEPT_ASSERT((boost::math::concepts::RealTypeConcept<boost::math::mpfr_interval>));
+
+ boost::math::mpfr_interval v;
+ v = 0.5f;
+ v -= 0.5f;
+}
+
Added: sandbox/interval_math_toolkit/libs/math/tools/non_central_beta_data.cpp
==============================================================================
--- (empty file)
+++ sandbox/interval_math_toolkit/libs/math/tools/non_central_beta_data.cpp 2008-04-04 11:56:10 EDT (Fri, 04 Apr 2008)
@@ -0,0 +1,281 @@
+// (C) Copyright John Maddock 2006.
+// 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)
+
+#include <boost/math/bindings/mpfr_interval.hpp>
+#include <boost/test/included/test_exec_monitor.hpp>
+#include <boost/math/special_functions.hpp>
+#include <boost/math/distributions.hpp>
+#include <boost/math/tools/test.hpp>
+#include <boost/lexical_cast.hpp>
+#include <fstream>
+#include <map>
+
+#include <boost/math/tools/test_data.hpp>
+
+using namespace boost::math::tools;
+using namespace boost::math;
+using namespace std;
+
+template <class T, class Policy>
+T non_central_beta_p(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 = boost::math::tools::real_cast<int>(l2);
+ // Starting Poisson weight:
+ T pois = gamma_p_derivative(k+1, l2, pol);
+ // Starting beta term:
+ T beta = ibeta(a + k, b, x, pol);
+ // recurance term:
+ T xterm = ibeta_derivative(a + k, b, x, pol) * (1 - x) / (a + b + k - 1);
+ T poisf(pois), betaf(beta), xtermf(xterm);
+ T sum = 0;
+
+ //
+ // Backwards recursion first, this is the stable
+ // direction for recursion:
+ //
+ for(int i = k; i >= 0; --i)
+ {
+ //T beta_c = ibeta(a + i, b, x);
+ //T pois_c = pdf(poisson_distribution<T>(l2), i);
+ T term = beta * pois;
+ //std::cout << i << " " << term << std::endl;
+ sum += term;
+ if(fabs(term/sum) < errtol)
+ break;
+ pois *= i / l2;
+ beta += xterm;
+ xterm *= (a + i - 1) / (x * (a + b + i - 2));
+ }
+ for(int i = k + 1; i - k < max_iter; ++i)
+ {
+ poisf *= l2 / i;
+ xtermf *= (x * (a + b + i - 2)) / (a + i - 1);
+ betaf -= xtermf;
+ //T beta_c = ibeta(a + i, b, x);
+ //T pois_c = pdf(poisson_distribution<T>(l2), i);
+
+ T term = poisf * betaf;
+ //std::cout << i << " " << term << std::endl;
+ sum += term;
+ if(fabs(term/sum) < errtol)
+ break;
+ }
+ return sum;
+}
+
+template <class T, class Policy>
+T non_central_beta_q(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 = boost::math::tools::real_cast<int>(l2);
+ // Starting Poisson weight:
+ T pois = gamma_p_derivative(k+1, l2, pol);
+ // Starting beta term:
+ T beta = ibetac(a + k, b, x, pol);
+ // recurance term:
+ T xterm = ibeta_derivative(a + k, b, x, pol) * (1 - x) / (a + b + k - 1);
+ T poisf(pois), betaf(beta), xtermf(xterm);
+ T sum = 0;
+
+ //
+ // Backwards recursion first, this is the unstable
+ // direction for recursion:
+ //
+ for(int i = k; i >= 0; --i)
+ {
+ T beta_c = ibetac(a + i, b, x);
+ T pois_c = pdf(poisson_distribution<T>(l2), i);
+ T term = beta * pois;
+ //std::cout << i << " " << term << std::endl;
+ sum += term;
+ if(fabs(term/sum) < errtol)
+ break;
+ pois *= i / l2;
+ beta -= xterm;
+ xterm *= (a + i - 1) / (x * (a + b + i - 2));
+ }
+ for(int i = k + 1; i - k < max_iter; ++i)
+ {
+ poisf *= l2 / i;
+ xtermf *= (x * (a + b + i - 2)) / (a + i - 1);
+ betaf += xtermf;
+ T beta_c = ibeta(a + i, b, x);
+ T pois_c = pdf(poisson_distribution<T>(l2), i);
+
+ T term = poisf * betaf;
+ //std::cout << i << " " << term << std::endl;
+ sum += term;
+ if(fabs(term/sum) < errtol)
+ break;
+ }
+ return sum;
+}
+
+
+int i = 0;
+int j = 4;
+float factors[] =
+{
+ 1/10.0F,
+ 0.25f,
+ 0.5f,
+ 0.75f,
+ 0.9f,
+ 0.99f,
+ 0.999f,
+ 1.0f,
+ 1.0001f,
+ 1.001f,
+ 1.01f,
+ 1.1f,
+ 1.2f,
+ 1.3f,
+ 1.5f,
+ 2,
+ 3,
+ 4,
+ 5
+};
+
+volatile float external_f;
+float force_truncate(const float* f)
+{
+ external_f = *f;
+ return external_f;
+}
+
+struct non_central_beta_cdf
+{
+ std::tr1::tuple<mpfr_class, mpfr_class, mpfr_class, mpfr_class, mpfr_class, mpfr_class>
+ operator()(
+ mpfr_class alpha,
+ mpfr_class ncp)
+ {
+ float factor = factors[j % (sizeof(factors)/sizeof(factors[0]))];
+ ++j;
+ float beta = factor * boost::math::tools::real_cast<float>(alpha);
+ factor = factors[i % (sizeof(factors)/sizeof(factors[0]))];
+ ++i;
+ mpfr_class c = alpha + beta + ncp / 2;
+ mpfr_class cross = 1 - (beta / c) * (1 + ncp / (2 * c * c));
+ float x;
+ if(factor < 1)
+ x = factor * boost::math::tools::real_cast<float>(cross);
+ else
+ x = 1 - (1 - boost::math::tools::real_cast<float>(cross)) / factor;
+ mpfr_class p, q;
+ //
+ // Use a Policy that ups the max number of iterations, so that
+ // we should get convergence even in extreme cases!
+ //
+ boost::math::policies::policy<boost::math::policies::max_series_iterations<BOOST_MATH_MAX_SERIES_ITERATION_POLICY*10> > pol;
+ std::cout << alpha << " " << beta << " " << ncp << " " << x << "(x" << factor << ") ";
+ if(x < cross)
+ {
+ mpfr_interval pi = non_central_beta_p(
+ mpfr_interval(alpha), mpfr_interval(beta),
+ mpfr_interval(ncp), mpfr_interval(x), pol);
+
+ std::cout << median(pi) << "("
+ << width(pi)/median(pi) << ") ";
+ std::cout << median(1-pi) << "("
+ << width(1-pi)/median(1-pi) << ")" << std::endl;
+
+ if(width(pi)/median(pi) > 1e-50)
+ throw std::domain_error("Error too large");
+ if(width(1-pi)/median(1-pi) > 1e-50)
+ throw std::domain_error("Error too large");
+ p = median(pi);
+ q = median(1-pi);
+ }
+ else
+ {
+ mpfr_interval qi = non_central_beta_q(
+ mpfr_interval(alpha), mpfr_interval(beta),
+ mpfr_interval(ncp), mpfr_interval(x), pol);
+
+ std::cout << median(1-qi) << "("
+ << width(1-qi)/median(1-qi) << ") ";
+ std::cout << median(qi) << "("
+ << width(qi)/median(qi) << ")" << std::endl;
+
+ if(width(qi)/median(qi) > 1e-50)
+ throw std::domain_error("Error too large");
+ if(width(1-qi)/median(1-qi) > 1e-50)
+ throw std::domain_error("Error too large");
+ p = median(1-qi);
+ q = median(qi);
+ }
+ return std::tr1::make_tuple(alpha, beta, ncp, x, p, q);
+ }
+};
+
+int test_main(int argc, char*argv [])
+{
+ using namespace boost::math::tools;
+
+ mpfr_class::set_dprec(300);
+
+ parameter_info<mpfr_class> arg1, arg2;
+ test_data<mpfr_class> data;
+ std::string line;
+ bool cont = false;
+
+ std::cout << "Welcome.\n"
+ "This program will generate spot tests for the non central beta CDF:\n\n";
+ do{
+ //get_user_parameter_info(arg1, "alpha");
+ //get_user_parameter_info(arg2, "non centrality");
+ arg1 = make_power_param(mpfr_class(0), 7, 15);
+ arg2 = make_power_param(mpfr_class(0), 7, 15);
+ arg1.type |=dummy_param;
+ arg2.type |=dummy_param;
+ data.insert(non_central_beta_cdf(), arg1, arg2);
+ /*
+ std::cout << "Any more data [y/n]?";
+ std::getline(std::cin, line);
+ boost::algorithm::trim(line);
+ */
+ cont = false; //(line == "y");
+ }while(cont);
+
+ /*
+ std::cout << "Enter name of test data file [default=ncbeta.ipp]";
+ std::getline(std::cin, line);
+ boost::algorithm::trim(line);
+ if(line == "")
+ line = "ncbeta.ipp";
+ */
+ std::ofstream ofs("ncbeta_big.ipp");
+ ofs << std::setprecision(40);
+ write_code(ofs, data, "ncbeta_big");
+
+ return 0;
+}
+
+
Added: sandbox/interval_math_toolkit/libs/math/tools/non_central_chi_sq_data.cpp
==============================================================================
--- (empty file)
+++ sandbox/interval_math_toolkit/libs/math/tools/non_central_chi_sq_data.cpp 2008-04-04 11:56:10 EDT (Fri, 04 Apr 2008)
@@ -0,0 +1,267 @@
+// (C) Copyright John Maddock 2006.
+// 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)
+
+#include <boost/math/bindings/mpfr_interval.hpp>
+#include <boost/test/included/test_exec_monitor.hpp>
+#include <boost/math/special_functions/gamma.hpp>
+#include <boost/math/tools/test.hpp>
+#include <boost/lexical_cast.hpp>
+#include <fstream>
+#include <map>
+
+#include <boost/math/tools/test_data.hpp>
+
+using namespace boost::math::tools;
+using namespace boost::math;
+using namespace std;
+
+template <class T, class Policy>
+T non_central_chi_square_q(T x, T f, T theta, const Policy& pol)
+{
+ //
+ // Computes the complement of the Non-Central Chi-Square
+ // Distribution CDF by summing a weighted sum of complements
+ // of the central-distributions. The weighting factor is
+ // a Poisson Distribution.
+ //
+ BOOST_MATH_STD_USING
+ using namespace boost::math;
+
+ T lambda = theta / 2;
+ T del = f / 2;
+ T y = x / 2;
+ boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
+ T errtol = ldexp(1.0, -boost::math::policies::digits<T, Policy>());
+ T sum = 0;
+ int k = tools::real_cast<int>(lambda);
+ T poisf = boost::math::gamma_p_derivative(1 + k, lambda);
+ T poisb = poisf * k / lambda;
+ T gamf = boost::math::gamma_q(del + k, y);
+ T xtermf = boost::math::gamma_p_derivative(del + 1 + k, y);
+ T xtermb = xtermf * (del + k) / y;
+ T gamb = gamf - xtermb;
+
+
+ for(int i = k; i < max_iter; ++i)
+ {
+ T term = poisf * gamf;
+ sum += term;
+ poisf *= lambda / (i + 1);
+ gamf += xtermf;
+ xtermf *= y / (del + i + 1);
+ if(tools::maybe_less(term / sum, errtol))
+ break;
+ }
+ for(int i = k - 1; i >= 0; --i)
+ {
+ T term = poisb * gamb;
+ sum += term;
+ poisb *= i / lambda;
+ xtermb *= (del + i) / y;
+ gamb -= xtermb;
+ if(tools::maybe_less(term / sum, errtol))
+ break;
+ }
+
+ return sum;
+}
+
+template <class T, class Policy>
+T non_central_chi_square_p(T y, T n, T lambda, const Policy& pol)
+{
+ using namespace boost::math;
+ 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 errorf, errorb;
+
+ T x = y / 2;
+ T del = lambda / 2;
+
+ int k = //tools::real_cast<int>(detail::inverse_poisson_cornish_fisher(del, 1 - errtol * 100, errtol * 100, pol));
+ tools::real_cast<int>(del);
+ T a = n / 2 + k;
+ T gamkf = boost::math::gamma_p(a, x, pol);
+ BOOST_MATH_INSTRUMENT_VARIABLE(gamkf);
+
+ if(lambda == 0)
+ return gamkf;
+
+ T gamkb = gamkf;
+ T poiskf = gamma_p_derivative(k+1, del, pol);
+ BOOST_MATH_INSTRUMENT_VARIABLE(poiskf)
+
+ T poiskb = poiskf;
+ T xtermf = boost::math::gamma_p_derivative(a, x, pol);
+ BOOST_MATH_INSTRUMENT_VARIABLE(xtermf)
+
+ T xtermb = xtermf * x / a;
+ T sum = poiskf * gamkf;
+ int i = 1;
+ //bool recurse_fwd(true), recurse_bwd(true);
+
+ while(i <= k)
+ {
+ xtermb *= (a - i + 1) / x;
+ gamkb += xtermb;
+ poiskb = poiskb * (k - i + 1) / del;
+ errorb = gamkb * poiskb;
+ sum += errorb;
+ if(tools::maybe_less(errorb / sum, errtol))
+ break;
+ ++i;
+ }
+ i = 1;
+ do
+ {
+ xtermf = xtermf * x / (a + i - 1);
+ gamkf = gamkf - xtermf;
+ poiskf = poiskf * del / (k + i);
+ errorf = poiskf * gamkf;
+ sum += errorf;
+ ++i;
+ }while(!tools::maybe_less(errorf / sum, errtol));
+
+ return sum;
+}
+
+int i = 0;
+float factors[] =
+{
+ 1/10.0F,
+ 0.25f,
+ 0.5f,
+ 0.75f,
+ 0.9f,
+ 0.99f,
+ 0.999f,
+ 1.0f,
+ 1.0001f,
+ 1.001f,
+ 1.01f,
+ 1.1f,
+ 1.2f,
+ 1.3f,
+ 1.5f,
+ 2,
+ 3,
+ 4,
+ 5
+};
+
+volatile float external_f;
+float force_truncate(const float* f)
+{
+ external_f = *f;
+ return external_f;
+}
+
+struct non_central_chi_square_cdf
+{
+ std::tr1::tuple<mpfr_class, mpfr_class, mpfr_class>
+ operator()(
+ mpfr_class df,
+ mpfr_class ncp)
+ {
+ float factor = factors[i % (sizeof(factors)/sizeof(factors[0]))];
+ ++i;
+ float x = factor * boost::math::tools::real_cast<float>(mpfr_class(df + ncp));
+ mpfr_class p, q;
+ boost::math::policies::policy<> pol;
+ std::cout << df << " " << ncp << " " << x << "(x" << factor << ") ";
+ if(x < df + ncp)
+ {
+ mpfr_interval pi = non_central_chi_square_p(
+ mpfr_interval(x), mpfr_interval(df), mpfr_interval(ncp), pol);
+
+ std::cout << median(pi) << "("
+ << width(pi)/median(pi) << ") ";
+ std::cout << median(1-pi) << "("
+ << width(1-pi)/median(1-pi) << ")" << std::endl;
+
+ if(width(pi)/median(pi) > 1e-50)
+ throw std::domain_error("Error too large");
+ if(width(1-pi)/median(1-pi) > 1e-50)
+ throw std::domain_error("Error too large");
+ p = median(pi);
+ q = median(1-pi);
+ }
+ else
+ {
+ mpfr_interval qi = non_central_chi_square_q(
+ mpfr_interval(x), mpfr_interval(df), mpfr_interval(ncp), pol);
+
+ std::cout << median(1-qi) << "("
+ << width(1-qi)/median(1-qi) << ") ";
+ std::cout << median(qi) << "("
+ << width(qi)/median(qi) << ")" << std::endl;
+
+ if(width(qi)/median(qi) > 1e-50)
+ throw std::domain_error("Error too large");
+ if(width(1-qi)/median(1-qi) > 1e-50)
+ throw std::domain_error("Error too large");
+ p = median(1-qi);
+ q = median(qi);
+ }
+ return std::tr1::make_tuple(x, p, q);
+ }
+};
+
+int test_main(int argc, char*argv [])
+{
+ using namespace boost::math::tools;
+
+ mpfr_class::set_dprec(300);
+
+ parameter_info<mpfr_class> arg1, arg2;
+ test_data<mpfr_class> data;
+ std::string line;
+ bool cont = false;
+
+ if((argc >= 2) && (std::strcmp(argv[1], "-inverse") == 0))
+ {
+ std::cout << "Welcome.\n"
+ "This program will generate spot tests for the non central chi square CDF:\n";
+
+ do{
+ get_user_parameter_info(arg1, "degrees of freedom");
+ get_user_parameter_info(arg2, "non centrality");
+ data.insert(non_central_chi_square_cdf(), arg1, arg2);
+
+ std::cout << "Any more data [y/n]?";
+ std::getline(std::cin, line);
+ boost::algorithm::trim(line);
+ cont = (line == "y");
+ }while(cont);
+ }
+ else
+ {
+ std::cout << "Welcome.\n"
+ "This program will generate spot tests for the non central chi square CDF:\n\n";
+ do{
+ get_user_parameter_info(arg1, "degrees of freedom");
+ get_user_parameter_info(arg2, "non centrality");
+ data.insert(non_central_chi_square_cdf(), arg1, arg2);
+
+ std::cout << "Any more data [y/n]?";
+ std::getline(std::cin, line);
+ boost::algorithm::trim(line);
+ cont = (line == "y");
+ }while(cont);
+ }
+
+ std::cout << "Enter name of test data file [default=nccs.ipp]";
+ std::getline(std::cin, line);
+ boost::algorithm::trim(line);
+ if(line == "")
+ line = "nccs.ipp";
+ std::ofstream ofs(line.c_str());
+ ofs << std::setprecision(40);
+ write_code(ofs, data, "nccs");
+
+ return 0;
+}
+
+
Added: sandbox/interval_math_toolkit/libs/math/tools/non_central_t_data.cpp
==============================================================================
--- (empty file)
+++ sandbox/interval_math_toolkit/libs/math/tools/non_central_t_data.cpp 2008-04-04 11:56:10 EDT (Fri, 04 Apr 2008)
@@ -0,0 +1,467 @@
+// (C) Copyright John Maddock 2006.
+// 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)
+
+#include <boost/math/bindings/mpfr_interval.hpp>
+#include <boost/test/included/test_exec_monitor.hpp>
+#include <boost/math/special_functions.hpp>
+#include <boost/math/distributions.hpp>
+#include <boost/math/tools/test.hpp>
+#include <boost/lexical_cast.hpp>
+#include <fstream>
+#include <map>
+
+#include <boost/math/tools/test_data.hpp>
+
+using namespace boost::math::tools;
+using namespace boost::math;
+using namespace std;
+
+template <class T, class Policy>
+T non_central_beta_p(T a, T b, T lam, T x, T y, 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 = tools::real_cast<int>(l2);
+ // Starting Poisson weight:
+ T pois = gamma_p_derivative(k+1, l2, pol);
+ // Starting beta term:
+ T beta = x < y
+ ? ibeta(a + k, b, x, pol)
+ : ibetac(b, a + k, y, pol);
+ // recurance term:
+ T xterm = x < y
+ ? ibeta_derivative(a + k, b, x, pol)
+ : ibeta_derivative(b, a + k, y, pol);
+ xterm *= (1 - x) / (a + b + k - 1);
+ T poisf(pois), betaf(beta), xtermf(xterm);
+ T sum = 0;
+
+ //
+ // Backwards recursion first, this is the stable
+ // direction for recursion:
+ //
+ int count = 0;
+ for(int i = k; i >= 0; --i)
+ {
+ T term = beta * pois;
+ sum += term;
+ if(tools::maybe_less(fabs(term/sum), errtol))
+ break;
+ pois *= i / l2;
+ beta += xterm;
+ xterm *= (a + i - 1) / (x * (a + b + i - 2));
+ ++count;
+ }
+ for(int i = k + 1; i - k < max_iter; ++i)
+ {
+ poisf *= l2 / i;
+ xtermf *= (x * (a + b + i - 2)) / (a + i - 1);
+ betaf -= xtermf;
+ T term = poisf * betaf;
+ sum += term;
+ if(tools::maybe_less(fabs(term/sum), errtol))
+ break;
+ ++count;
+ }
+ return sum;
+}
+
+template <class T, class Policy>
+T non_central_beta_q(T a, T b, T lam, T x, T y, 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 = tools::real_cast<int>(l2);
+ // Starting Poisson weight:
+ T pois = gamma_p_derivative(k+1, l2, pol);
+ // Starting beta term:
+ T beta = x < y
+ ? ibetac(a + k, b, x, pol)
+ : ibeta(b, a + k, y, pol);
+ // Recurance term:
+ T xterm = x < y
+ ? ibeta_derivative(a + k, b, x, pol)
+ : ibeta_derivative(b, a + k, y, pol);
+ xterm *= (1 - x) / (a + b + k - 1);
+ T poisf(pois), betaf(beta), xtermf(xterm);
+ T sum = 0;
+
+ //
+ // Backwards recursion first, this is the unstable
+ // direction for recursion:
+ //
+ for(int i = k; i >= 0; --i)
+ {
+ T term = beta * pois;
+ sum += term;
+ if(tools::maybe_less(fabs(term/sum), errtol))
+ break;
+ pois *= i / l2;
+ beta -= xterm;
+ xterm *= (a + i - 1) / (x * (a + b + i - 2));
+ }
+ for(int i = k + 1; i - k < max_iter; ++i)
+ {
+ poisf *= l2 / i;
+ xtermf *= (x * (a + b + i - 2)) / (a + i - 1);
+ betaf += xtermf;
+ T term = poisf * betaf;
+ sum += term;
+ if(tools::maybe_less(fabs(term/sum), errtol))
+ break;
+ }
+ return sum;
+}
+
+template <class T, class Policy>
+T non_central_t2_p(T n, T delta, T x, T y, const Policy& pol, T init_val)
+{
+ 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 d2 = delta * delta / 2;
+ //
+ // k is the starting point for iteration, and is the
+ // maximum of the poisson weighting term:
+ //
+ int k = tools::real_cast<int>(d2);
+ // Starting Poisson weight:
+ T pois = gamma_p_derivative(k+1, d2, pol)
+ * tgamma_delta_ratio(k + 1, T(0.5f))
+ * delta / constants::root_two<T>();
+ // Starting beta term:
+ T beta = x < y
+ ? ibeta(k + 1, n / 2, x, pol)
+ : ibetac(n / 2, k + 1, y, pol);
+ // Recurance term:
+ T xterm = x < y
+ ? ibeta_derivative(k + 1, n / 2, x, pol)
+ : ibeta_derivative(n / 2, k + 1, y, pol);
+ xterm *= (1 - x) / (n / 2 + k);
+ T poisf(pois), betaf(beta), xtermf(xterm);
+ T sum = init_val;
+
+ //
+ // Backwards recursion first, this is the stable
+ // direction for recursion:
+ //
+ int count = 0;
+ for(int i = k; i >= 0; --i)
+ {
+ T term = beta * pois;
+ sum += term;
+ if(tools::maybe_less(fabs(term/sum), errtol))
+ break;
+ pois *= (i + 0.5f) / d2;
+ beta += xterm;
+ xterm *= (i) / (x * (n / 2 + i - 1));
+ ++count;
+ }
+ for(int i = k + 1; i - k < max_iter; ++i)
+ {
+ poisf *= d2 / (i + 0.5f);
+ xtermf *= (x * (n / 2 + i - 1)) / (i);
+ betaf -= xtermf;
+ T term = poisf * betaf;
+ sum += term;
+ if(tools::maybe_less(fabs(term/sum), errtol))
+ break;
+ ++count;
+ }
+ return sum;
+}
+
+template <class T, class Policy>
+T non_central_t2_q(T n, T delta, T x, T y, const Policy& pol, T init_val)
+{
+ 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 d2 = delta * delta / 2;
+ //
+ // k is the starting point for iteration, and is the
+ // maximum of the poisson weighting term:
+ //
+ int k = tools::real_cast<int>(d2);
+ // Starting Poisson weight:
+ T pois = gamma_p_derivative(k+1, d2, pol)
+ * tgamma_delta_ratio(k + 1, T(0.5f))
+ * delta / constants::root_two<T>();
+ // Starting beta term:
+ T beta = x < y
+ ? ibetac(k + 1, n / 2, x, pol)
+ : ibeta(n / 2, k + 1, y, pol);
+ // Recurance term:
+ T xterm = x < y
+ ? ibeta_derivative(k + 1, n / 2, x, pol)
+ : ibeta_derivative(n / 2, k + 1, y, pol);
+ xterm *= (1 - x) / (n / 2 + k);
+ T poisf(pois), betaf(beta), xtermf(xterm);
+ T sum = init_val;
+
+ //
+ // Forward recursion first, this is the stable direction:
+ //
+ int count = 0;
+ for(int i = k + 1; i - k < max_iter; ++i)
+ {
+ poisf *= d2 / (i + 0.5f);
+ xtermf *= (x * (n / 2 + i - 1)) / (i);
+ betaf += xtermf;
+
+ T term = poisf * betaf;
+ sum += term;
+ if(tools::maybe_less(fabs(term/sum), errtol))
+ break;
+ ++count;
+ }
+ //
+ // Backwards recursion:
+ //
+ for(int i = k; i >= 0; --i)
+ {
+ T term = beta * pois;
+ sum += term;
+ if(tools::maybe_less(fabs(term/sum), errtol))
+ break;
+ pois *= (i + 0.5f) / d2;
+ beta -= xterm;
+ xterm *= (i) / (x * (n / 2 + i - 1));
+ ++count;
+ }
+ return sum;
+}
+
+template <class T, class Policy>
+T get_non_central_t_cdf(T n, T delta, T t, bool invert, const Policy& pol)
+{
+ if(t < T(0))
+ {
+ t = -t;
+ delta = -delta;
+ invert = !invert;
+ }
+ T x = t * t / (n + t * t);
+ T y = n / (n + t * t);
+ T d2 = delta * delta;
+ T a = 0.5f;
+ T b = n / 2;
+ T c = a + b + d2 / 2;
+ T cross = 1 - (b / c) * (1 + d2 / (2 * c * c));
+ T result;
+ if(x < cross)
+ {
+ //
+ // Calculate p:
+ //
+ result = non_central_beta_p(a, b, d2, x, y, pol);
+ result = non_central_t2_p(n, delta, x, y, pol, result);
+ result /= 2;
+ result += cdf(boost::math::normal_distribution<T, Policy>(), -delta);
+ }
+ else
+ {
+ //
+ // Calculate q:
+ //
+ invert = !invert;
+ result = non_central_beta_q(a, b, d2, x, y, pol);
+ result = non_central_t2_q(n, delta, x, y, pol, result);
+ result /= 2;
+ }
+ if(invert)
+ result = 1 - result;
+ return result;
+}
+
+template <class T, class Policy>
+T mean(T v, T delta, const Policy& pol)
+{
+ return v <= 1 ? 0 : delta * sqrt(v / 2) * tgamma_delta_ratio((v - 1) * 0.5f, T(0.5f), pol);
+}
+
+template <class T, class Policy>
+T variance(T v, T delta, const Policy& pol)
+{
+ if(v <= 2)
+ return 0.1;
+ T result = ((delta * delta + 1) * v) / (v - 2);
+ T m = mean(v, delta, pol);
+ result -= m * m;
+ return result;
+}
+
+int i = 0;
+int j = 4;
+float factors[] =
+{
+ 1/10.0F,
+ 0.25f,
+ 0.5f,
+ 0.75f,
+ 0.9f,
+ 0.99f,
+ 0.999f,
+ 1.0f,
+ 1.0001f,
+ 1.001f,
+ 1.01f,
+ 1.1f,
+ 1.2f,
+ 1.3f,
+ 1.5f,
+ 2,
+ 3,
+ 4,
+ 5
+};
+
+int signs[] = {
+1,1,1,-1,1,-1,-1,-1,1,1,-1,1,-1,-1,-1,-1,1,1,1,-1,-1,1,-1,1,-1,1,1,-1,1,-1,-1,-1,-1,-1,1,1
+};
+
+volatile float external_f;
+float force_truncate(const float* f)
+{
+ external_f = *f;
+ return external_f;
+}
+
+struct non_central_t_cdf
+{
+ std::tr1::tuple<mpfr_class, mpfr_class, mpfr_class, mpfr_class, mpfr_class>
+ operator()(
+ mpfr_class v)
+ {
+ //
+ // Use a Policy that ups the max number of iterations, so that
+ // we should get convergence even in extreme cases!
+ //
+ boost::math::policies::policy<boost::math::policies::max_series_iterations<BOOST_MATH_MAX_SERIES_ITERATION_POLICY*10> > pol;
+
+ float factor = factors[j % (sizeof(factors)/sizeof(factors[0]))];
+ ++j;
+ float delta = signs[j % (sizeof(signs)/sizeof(signs[0]))] * factor * boost::math::tools::real_cast<float>(v);
+ factor = factors[i % (sizeof(factors)/sizeof(factors[0]))];
+ ++i;
+ float m = mean(boost::math::tools::real_cast<float>(v), delta, pol);
+ float sd = sqrt(variance(boost::math::tools::real_cast<float>(v), delta, pol));
+ float x = m + sd * factor * signs[j % (sizeof(signs)/sizeof(signs[0]))];
+
+ mpfr_class p, q;
+ std::cout << v << " " << delta << " " << x << "(x" << factor << ") ";
+ if(x < m)
+ {
+ mpfr_interval pi = get_non_central_t_cdf(
+ mpfr_interval(v), mpfr_interval(delta),
+ mpfr_interval(x), false, pol);
+
+ std::cout << median(pi) << "("
+ << width(pi)/median(pi) << ") ";
+ std::cout << median(1-pi) << "("
+ << width(1-pi)/median(1-pi) << ")" << std::endl;
+
+ if(width(pi)/median(pi) > 1e-50)
+ throw std::domain_error("Error too large");
+ if(width(1-pi)/median(1-pi) > 1e-50)
+ throw std::domain_error("Error too large");
+ p = median(pi);
+ q = median(1-pi);
+ }
+ else
+ {
+ mpfr_interval qi = get_non_central_t_cdf(
+ mpfr_interval(v), mpfr_interval(delta),
+ mpfr_interval(x), true, pol);
+
+ std::cout << median(1-qi) << "("
+ << width(1-qi)/median(1-qi) << ") ";
+ std::cout << median(qi) << "("
+ << width(qi)/median(qi) << ")" << std::endl;
+
+ if(width(qi)/median(qi) > 1e-50)
+ throw std::domain_error("Error too large");
+ if(width(1-qi)/median(1-qi) > 1e-50)
+ throw std::domain_error("Error too large");
+ p = median(1-qi);
+ q = median(qi);
+ }
+ return std::tr1::make_tuple(v, delta, x, p, q);
+ }
+};
+
+int test_main(int argc, char*argv [])
+{
+ using namespace boost::math::tools;
+
+ mpfr_class::set_dprec(300);
+
+ parameter_info<mpfr_class> arg1, arg2;
+ test_data<mpfr_class> data;
+ std::string line;
+ bool cont = false;
+
+ std::cout << "Welcome.\n"
+ "This program will generate spot tests for the non central T CDF:\n\n";
+ do{
+ //get_user_parameter_info(arg1, "alpha");
+ //get_user_parameter_info(arg2, "non centrality");
+ arg1 = make_power_param(mpfr_class(0), 7, 15);
+ arg2 = make_random_param(mpfr_class(0.0001), mpfr_class(200), 200);
+ arg1.type |=dummy_param;
+ arg2.type |=dummy_param;
+ data.insert(non_central_t_cdf(), arg2);
+ data.insert(non_central_t_cdf(), arg1);
+ /*
+ std::cout << "Any more data [y/n]?";
+ std::getline(std::cin, line);
+ boost::algorithm::trim(line);
+ */
+ cont = false; //(line == "y");
+ }while(cont);
+
+ /*
+ std::cout << "Enter name of test data file [default=ncbeta.ipp]";
+ std::getline(std::cin, line);
+ boost::algorithm::trim(line);
+ if(line == "")
+ line = "ncbeta.ipp";
+ */
+ std::ofstream ofs("nct_big.ipp");
+ ofs << std::setprecision(40);
+ write_code(ofs, data, "nct_big");
+
+ return 0;
+}
+
+
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