boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity

Fix bug in incomplete beta inverse estimation routine (when estimating from student's t).

Browse Source 
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.

[SVN r85572]
This commit is contained in:
John Maddock
2013-09-05 15:56:17 +00:00
parent 041dffcf2e
commit a3e3e86eaf
6 changed files with 332 additions and 17 deletions
Download Patch File Download Diff File Expand all files Collapse all files

43
include/boost/math/special_functions/beta.hpp
View File

@@ -935,6 +935,49 @@ T ibeta_imp(T a, T b, T x, const Policy& pol, bool inv, bool normalised, T* p_de
}
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)
{

71
include/boost/math/special_functions/detail/ibeta_inverse.hpp
View File

@@ -454,6 +454,11 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
{
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:
//
@@ -467,17 +472,19 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
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 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
// 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 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
// 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)
{
//

9
include/boost/math/special_functions/detail/t_distribution_inv.hpp
View File

@@ -416,15 +416,14 @@ calculate_real:
}
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>

67
test/test_ibeta.hpp
View File

@@ -295,5 +295,72 @@ void test_spots(T)
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);
}

84
test/test_ibeta_inv.hpp
View File

@@ -177,5 +177,89 @@ void test_spots(T)
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);
}

75
tools/tgamma_large_data.cpp Normal file
View File

@@ -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;
}
//]