|
|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r85572 - in trunk: boost/math/special_functions boost/math/special_functions/detail libs/math/test libs/math/tools
From: john_at_[hidden]
Date: 2013-09-05 11:56:17
Author: johnmaddock
Date: 2013-09-05 11:56:17 EDT (Thu, 05 Sep 2013)
New Revision: 85572
URL: http://svn.boost.org/trac/boost/changeset/85572
Log:
Fix bug in incomplete beta inverse estimation routine (when estimating from student's t).
Add special cases to incomplete beta and inverse for a=b=0.5 and b=1.
Added tool for generating high precision gamma function test values.
Added:
trunk/libs/math/tools/tgamma_large_data.cpp (contents, props changed)
Text files modified:
trunk/boost/math/special_functions/beta.hpp | 43 ++++++++++++++++++++
trunk/boost/math/special_functions/detail/ibeta_inverse.hpp | 71 ++++++++++++++++++++++++++++-----
trunk/boost/math/special_functions/detail/t_distribution_inv.hpp | 9 +--
trunk/libs/math/test/test_ibeta.hpp | 67 +++++++++++++++++++++++++++++++
trunk/libs/math/test/test_ibeta_inv.hpp | 84 ++++++++++++++++++++++++++++++++++++++++
trunk/libs/math/tools/tgamma_large_data.cpp | 75 +++++++++++++++++++++++++++++++++++
6 files changed, 332 insertions(+), 17 deletions(-)
Modified: trunk/boost/math/special_functions/beta.hpp
==============================================================================
--- trunk/boost/math/special_functions/beta.hpp Thu Sep 5 08:57:15 2013 (r85571)
+++ trunk/boost/math/special_functions/beta.hpp 2013-09-05 11:56:17 EDT (Thu, 05 Sep 2013) (r85572)
@@ -935,6 +935,49 @@
}
return (invert == 0 ? (normalised ? 1 : boost::math::beta(a, b, pol)) : 0);
}
+ if((a == 0.5f) && (b == 0.5f))
+ {
+ // We have an arcsine distribution:
+ if(p_derivative)
+ {
+ *p_derivative = (invert ? -1 : 1) / constants::pi<T>() * sqrt(y * x);
+ }
+ T p = invert ? asin(sqrt(y)) / constants::half_pi<T>() : asin(sqrt(x)) / constants::half_pi<T>();
+ if(!normalised)
+ p *= constants::pi<T>();
+ return p;
+ }
+ if(a == 1)
+ {
+ std::swap(a, b);
+ std::swap(x, y);
+ invert = !invert;
+ }
+ if(b == 1)
+ {
+ //
+ // Special case see: http://functions.wolfram.com/GammaBetaErf/BetaRegularized/03/01/01/
+ //
+ if(a == 1)
+ {
+ if(p_derivative)
+ *p_derivative = invert ? -1 : 1;
+ return invert ? y : x;
+ }
+
+ if(p_derivative)
+ {
+ *p_derivative = (invert ? -1 : 1) * a * pow(x, a - 1);
+ }
+ T p;
+ if(y < 0.5)
+ p = invert ? T(-expm1(a * log1p(-y))) : T(exp(a * log1p(-y)));
+ else
+ p = invert ? T(-powm1(x, a)) : T(pow(x, a));
+ if(!normalised)
+ p /= a;
+ return p;
+ }
if((std::min)(a, b) <= 1)
{
Modified: trunk/boost/math/special_functions/detail/ibeta_inverse.hpp
==============================================================================
--- trunk/boost/math/special_functions/detail/ibeta_inverse.hpp Thu Sep 5 08:57:15 2013 (r85571)
+++ trunk/boost/math/special_functions/detail/ibeta_inverse.hpp 2013-09-05 11:56:17 EDT (Thu, 05 Sep 2013) (r85572)
@@ -455,6 +455,11 @@
BOOST_MATH_STD_USING // For ADL of math functions.
//
+ // The flag invert is set to true if we swap a for b and p for q,
+ // in which case the result has to be subtracted from 1:
+ //
+ bool invert = false;
+ //
// Handle trivial cases first:
//
if(q == 0)
@@ -467,17 +472,19 @@
if(py) *py = 1;
return 0;
}
- else if((a == 1) && (b == 1))
+ else if(a == 1)
{
- if(py) *py = 1 - p;
- return p;
+ if(b == 1)
+ {
+ if(py) *py = 1 - p;
+ return p;
+ }
+ // Change things around so we can handle as b == 1 special case below:
+ std::swap(a, b);
+ std::swap(p, q);
+ invert = true;
}
//
- // The flag invert is set to true if we swap a for b and p for q,
- // in which case the result has to be subtracted from 1:
- //
- bool invert = false;
- //
// Depending upon which approximation method we use, we may end up
// calculating either x or y initially (where y = 1-x):
//
@@ -495,11 +502,25 @@
// Student's T with b = 0.5 gets handled as a special case, swap
// around if the arguments are in the "wrong" order:
//
- if((a == 0.5f) && (b >= 0.5f))
+ if(a == 0.5f)
{
- std::swap(a, b);
- std::swap(p, q);
- invert = !invert;
+ if(b == 0.5f)
+ {
+ x = sin(p * constants::half_pi<T>());
+ x *= x;
+ if(py)
+ {
+ *py = sin(q * constants::half_pi<T>());
+ *py *= *py;
+ }
+ return x;
+ }
+ else if(b > 0.5f)
+ {
+ std::swap(a, b);
+ std::swap(p, q);
+ invert = !invert;
+ }
}
//
// Select calculation method for the initial estimate:
@@ -510,6 +531,32 @@
// We have a Student's T distribution:
x = find_ibeta_inv_from_t_dist(a, p, q, &y, pol);
}
+ else if(b == 1)
+ {
+ if(p < q)
+ {
+ if(a > 1)
+ {
+ x = pow(p, 1 / a);
+ y = -expm1(log(p) / a);
+ }
+ else
+ {
+ x = pow(p, 1 / a);
+ y = 1 - x;
+ }
+ }
+ else
+ {
+ x = exp(log1p(-q) / a);
+ y = -expm1(log1p(-q) / a);
+ }
+ if(invert)
+ std::swap(x, y);
+ if(py)
+ *py = y;
+ return x;
+ }
else if(a + b > 5)
{
//
Modified: trunk/boost/math/special_functions/detail/t_distribution_inv.hpp
==============================================================================
--- trunk/boost/math/special_functions/detail/t_distribution_inv.hpp Thu Sep 5 08:57:15 2013 (r85571)
+++ trunk/boost/math/special_functions/detail/t_distribution_inv.hpp 2013-09-05 11:56:17 EDT (Thu, 05 Sep 2013) (r85572)
@@ -416,15 +416,14 @@
}
template <class T, class Policy>
-inline T find_ibeta_inv_from_t_dist(T a, T p, T q, T* py, const Policy& pol)
+inline T find_ibeta_inv_from_t_dist(T a, T p, T /*q*/, T* py, const Policy& pol)
{
- T u = (p > q) ? T(0.5f - q) / T(2) : T(p / 2);
- T v = 1 - u; // u < 0.5 so no cancellation error
+ T u = p / 2;
+ T v = 1 - u;
T df = a * 2;
T t = boost::math::detail::inverse_students_t(df, u, v, pol);
- T x = df / (df + t * t);
*py = t * t / (df + t * t);
- return x;
+ return df / (df + t * t);
}
template <class T, class Policy>
Modified: trunk/libs/math/test/test_ibeta.hpp
==============================================================================
--- trunk/libs/math/test/test_ibeta.hpp Thu Sep 5 08:57:15 2013 (r85571)
+++ trunk/libs/math/test/test_ibeta.hpp 2013-09-05 11:56:17 EDT (Thu, 05 Sep 2013) (r85572)
@@ -295,5 +295,72 @@
BOOST_CHECK_THROW(::boost::math::ibetac(static_cast<T>(2), static_cast<T>(-2), static_cast<T>(0.5)), std::domain_error);
BOOST_CHECK_THROW(::boost::math::ibetac(static_cast<T>(2), static_cast<T>(2), static_cast<T>(-0.5)), std::domain_error);
BOOST_CHECK_THROW(::boost::math::ibetac(static_cast<T>(2), static_cast<T>(2), static_cast<T>(1.5)), std::domain_error);
+
+ //
+ // a = b = 0.5 is a special case:
+ //
+ BOOST_CHECK_CLOSE(
+ ::boost::math::ibeta(
+ static_cast<T>(0.5f),
+ static_cast<T>(0.5f),
+ static_cast<T>(0.25)),
+ static_cast<T>(1) / 3, tolerance);
+ BOOST_CHECK_CLOSE(
+ ::boost::math::ibetac(
+ static_cast<T>(0.5f),
+ static_cast<T>(0.5f),
+ static_cast<T>(0.25)),
+ static_cast<T>(2) / 3, tolerance);
+ BOOST_CHECK_CLOSE(
+ ::boost::math::ibeta(
+ static_cast<T>(0.5f),
+ static_cast<T>(0.5f),
+ static_cast<T>(0.125)),
+ static_cast<T>(0.230053456162615885213780567705142893009911395270714102055874L), tolerance);
+ BOOST_CHECK_CLOSE(
+ ::boost::math::ibetac(
+ static_cast<T>(0.5f),
+ static_cast<T>(0.5f),
+ static_cast<T>(0.125)),
+ static_cast<T>(0.769946543837384114786219432294857106990088604729285897944125L), tolerance);
+ BOOST_CHECK_CLOSE(
+ ::boost::math::ibeta(
+ static_cast<T>(0.5f),
+ static_cast<T>(0.5f),
+ static_cast<T>(0.825)),
+ static_cast<T>(0.725231121519469565327291851560156562956885802608457839260161L), tolerance);
+ BOOST_CHECK_CLOSE(
+ ::boost::math::ibetac(
+ static_cast<T>(0.5f),
+ static_cast<T>(0.5f),
+ static_cast<T>(0.825)),
+ static_cast<T>(0.274768878480530434672708148439843437043114197391542160739838L), tolerance);
+ //
+ // Second argument is 1 is a special case, see http://functions.wolfram.com/GammaBetaErf/BetaRegularized/03/01/01/
+ //
+ BOOST_CHECK_CLOSE(
+ ::boost::math::ibeta(
+ static_cast<T>(0.5f),
+ static_cast<T>(1),
+ static_cast<T>(0.825)),
+ static_cast<T>(0.908295106229247499626759842915458109758420750043003849691665L), tolerance);
+ BOOST_CHECK_CLOSE(
+ ::boost::math::ibetac(
+ static_cast<T>(0.5f),
+ static_cast<T>(1),
+ static_cast<T>(0.825)),
+ static_cast<T>(0.091704893770752500373240157084541890241579249956996150308334L), tolerance);
+ BOOST_CHECK_CLOSE(
+ ::boost::math::ibeta(
+ static_cast<T>(30),
+ static_cast<T>(1),
+ static_cast<T>(0.825)),
+ static_cast<T>(0.003116150729395132012981654047222541793435357905008020740211L), tolerance);
+ BOOST_CHECK_CLOSE(
+ ::boost::math::ibetac(
+ static_cast<T>(30),
+ static_cast<T>(1),
+ static_cast<T>(0.825)),
+ static_cast<T>(0.996883849270604867987018345952777458206564642094991979259788L), tolerance);
}
Modified: trunk/libs/math/test/test_ibeta_inv.hpp
==============================================================================
--- trunk/libs/math/test/test_ibeta_inv.hpp Thu Sep 5 08:57:15 2013 (r85571)
+++ trunk/libs/math/test/test_ibeta_inv.hpp 2013-09-05 11:56:17 EDT (Thu, 05 Sep 2013) (r85572)
@@ -177,5 +177,89 @@
static_cast<T>(80),
static_cast<T>(0.5)),
static_cast<T>(0.33240456430025026300937492802591128972548660643778L), tolerance);
+ BOOST_CHECK_CLOSE(
+ ::boost::math::ibeta_inv(
+ static_cast<T>(40),
+ static_cast<T>(0.5),
+ ldexp(T(1), -30)),
+ static_cast<T>(0.624305407878048788716096298053941618358257550305573588792717L), tolerance);
+ BOOST_CHECK_CLOSE(
+ ::boost::math::ibeta_inv(
+ static_cast<T>(40),
+ static_cast<T>(0.5),
+ 1 - ldexp(T(1), -30)),
+ static_cast<T>(0.99999999999999999998286262026583217516676792408012252456039L), tolerance);
+ BOOST_CHECK_CLOSE(
+ ::boost::math::ibeta_inv(
+ static_cast<T>(0.5),
+ static_cast<T>(40),
+ ldexp(T(1), -30)),
+ static_cast<T>(1.713737973416782483323207591987747543960774485649459249e-20L), tolerance);
+ BOOST_CHECK_CLOSE(
+ ::boost::math::ibeta_inv(
+ static_cast<T>(0.5),
+ static_cast<T>(0.75),
+ ldexp(T(1), -30)),
+ static_cast<T>(1.245132488513853853809715434621955746959615015005382639e-18L), tolerance);
+ BOOST_CHECK_CLOSE(
+ ::boost::math::ibeta_inv(
+ static_cast<T>(0.5),
+ static_cast<T>(0.5),
+ static_cast<T>(0.25)),
+ static_cast<T>(0.1464466094067262377995778189475754803575820311557629L), tolerance);
+ BOOST_CHECK_CLOSE(
+ ::boost::math::ibeta_inv(
+ static_cast<T>(0.5),
+ static_cast<T>(0.5),
+ static_cast<T>(0.75)),
+ static_cast<T>(0.853553390593273762200422181052424519642417968844237018294169L), tolerance);
+ BOOST_CHECK_CLOSE(
+ ::boost::math::ibeta_inv(
+ static_cast<T>(1),
+ static_cast<T>(5),
+ static_cast<T>(0.125)),
+ static_cast<T>(0.026352819384831863473794894078665766580641189002729204514544L), tolerance);
+ BOOST_CHECK_CLOSE(
+ ::boost::math::ibeta_inv(
+ static_cast<T>(5),
+ static_cast<T>(1),
+ static_cast<T>(0.125)),
+ static_cast<T>(0.659753955386447129687000985614820066516734506596709340752903L), tolerance);
+ BOOST_CHECK_CLOSE(
+ ::boost::math::ibeta_inv(
+ static_cast<T>(1),
+ static_cast<T>(0.125),
+ static_cast<T>(0.125)),
+ static_cast<T>(0.656391084194183349609374999999999999999999999999999999999999L), tolerance);
+ BOOST_CHECK_CLOSE(
+ ::boost::math::ibeta_inv(
+ static_cast<T>(0.125),
+ static_cast<T>(1),
+ static_cast<T>(0.125)),
+ static_cast<T>(5.960464477539062500000e-8), tolerance);
+ BOOST_CHECK_CLOSE(
+ ::boost::math::ibetac_inv(
+ static_cast<T>(5),
+ static_cast<T>(1),
+ static_cast<T>(0.125)),
+ static_cast<T>(0.973647180615168136526205105921334233419358810997270795485455L), tolerance);
+ BOOST_CHECK_CLOSE(
+ ::boost::math::ibetac_inv(
+ static_cast<T>(1),
+ static_cast<T>(5),
+ static_cast<T>(0.125)),
+ static_cast<T>(0.340246044613552870312999014385179933483265493403290659247096L), tolerance);
+ BOOST_CHECK_CLOSE(
+ ::boost::math::ibetac_inv(
+ static_cast<T>(0.125),
+ static_cast<T>(1),
+ static_cast<T>(0.125)),
+ static_cast<T>(0.343608915805816650390625000000000000000000000000000000000000L), tolerance);
+ BOOST_CHECK_CLOSE(
+ ::boost::math::ibetac_inv(
+ static_cast<T>(1),
+ static_cast<T>(0.125),
+ static_cast<T>(0.125)),
+ static_cast<T>(0.99999994039535522460937500000000000000000000000L), tolerance);
}
Added: trunk/libs/math/tools/tgamma_large_data.cpp
==============================================================================
--- /dev/null 00:00:00 1970 (empty, because file is newly added)
+++ trunk/libs/math/tools/tgamma_large_data.cpp 2013-09-05 11:56:17 EDT (Thu, 05 Sep 2013) (r85572)
@@ -0,0 +1,75 @@
+// Copyright John Maddock 2013.
+// 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)
+
+//[special_data_example
+
+#include <boost/multiprecision/mpfr.hpp>
+#include <boost/math/tools/test_data.hpp>
+#include <boost/test/included/prg_exec_monitor.hpp>
+#include <boost/math/tools/tuple.hpp>
+#include <fstream>
+
+using namespace boost::math::tools;
+using namespace boost::math;
+using namespace std;
+using namespace boost::multiprecision;
+
+typedef number<mpfr_float_backend<1000> > mp_type;
+
+
+boost::math::tuple<mp_type, mp_type, mp_type> generate(mp_type a)
+{
+ mp_type tg, lg;
+ mpfr_gamma(tg.backend().data(), a.backend().data(), GMP_RNDN);
+ mpfr_lngamma(lg.backend().data(), a.backend().data(), GMP_RNDN);
+ return boost::math::make_tuple(a, tg, lg);
+}
+
+int cpp_main(int argc, char*argv [])
+{
+ parameter_info<mp_type> arg1, arg2;
+ test_data<mp_type> data;
+
+ bool cont;
+ std::string line;
+
+ if(argc < 1)
+ return 1;
+
+ do{
+ //
+ // User interface which prompts for
+ // range of input parameters:
+ //
+ if(0 == get_user_parameter_info(arg1, "a"))
+ return 1;
+ arg1.type |= dummy_param;
+
+ data.insert(&generate, arg1);
+
+ std::cout << "Any more data [y/n]?";
+ std::getline(std::cin, line);
+ boost::algorithm::trim(line);
+ cont = (line == "y");
+ }while(cont);
+ //
+ // Just need to write the results to a file:
+ //
+ std::cout << "Enter name of test data file [default=gamma.ipp]";
+ std::getline(std::cin, line);
+ boost::algorithm::trim(line);
+ if(line == "")
+ line = "gamma.ipp";
+ std::ofstream ofs(line.c_str());
+ line.erase(line.find('.'));
+ ofs << std::scientific << std::setprecision(500);
+ write_code(ofs, data, line.c_str());
+
+ 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