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

Added finite sum to incomplete beta, plus tests and docs.

Browse Source 
Linked in lognormal docs.


[SVN r3310]
This commit is contained in:
John Maddock
2006-10-25 09:03:04 +00:00
parent 76f2122bad
commit 2f964f6ef6
11 changed files with 1188 additions and 37 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
doc/dist_reference.qbk
View File

@@ -36,6 +36,7 @@
[include distributions/exponential.qbk]
[include distributions/extreme_value.qbk]
[include distributions/fisher.qbk]
[include distributions/lognormal.qbk]
[include distributions/normal.qbk]
[include distributions/students_t.qbk]

76
doc/equations/ibeta.xml
View File

@@ -1405,5 +1405,81 @@
</msup>
</mrow>
</math>
<math xmlns="http://www.w3.org/1998/Math/MathML" display="block">
<mrow>
<msub>
<mi>I</mi>
<mi>x</mi>
</msub>
<mfenced>
<mrow>
<mi>a</mi>
<mo>,</mo>
<mi>b</mi>
</mrow>
</mfenced>
<mspace width="1em"/>
<mo>=</mo>
<mspace width="1em"/>
<munderover>
<mo>&Sum;</mo>
<mrow>
<mi>i</mi>
<mo>=</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
<mi>N</mi>
</munderover>
<mfrac>
<mrow>
<mi>N</mi>
<mo>!</mo>
</mrow>
<mrow>
<mi>i</mi>
<mo>!</mo>
<mfenced>
<mrow>
<mi>N</mi>
<mo>&minus;</mo>
<mi>i</mi>
</mrow>
</mfenced>
<mo>!</mo>
</mrow>
</mfrac>
<msup>
<mi>x</mi>
<mi>i</mi>
</msup>
<msup>
<mi>y</mi>
<mfenced>
<mrow>
<mi>N</mi>
<mo>&minus;</mo>
<mi>i</mi>
</mrow>
</mfenced>
</msup>
<mspace width="1em"/>
<mo>;</mo>
<mspace width="1em"/>
<mi>k</mi>
<mo>=</mo>
<mi>a</mi>
<mo>&minus;</mo>
<mn>1,</mn>
<mi>N</mi>
<mo>=</mo>
<mi>a</mi>
<mo>+</mo>
<mi>b</mi>
<mo>&minus;</mo>
<mn>1</mn>
</mrow>
</math>
</body>
</html>

BIN
doc/equations/ibeta12.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 1.4 KiB

5
doc/ibeta.qbk
View File

@@ -161,6 +161,11 @@ Note that this method relies
on keeping a table of all the p[sub n ] previously computed, which does limit
the precision of the method, depending upon the size of the table used.
When /a/ and /b/ are both small integers, then we can relate the incomplete
beta to the binomial distribution and use the following finite sum:
[$../equations/ibeta12.png]
Finally we can sidestep difficult areas, or move to an area with a more
efficient means of computation, by using the duplication formulae:

41
include/boost/math/distributions/binomial.hpp
View File

@@ -368,8 +368,8 @@ namespace boost
// -- ( i )
// i=0
// The terms are not summed directly (at least for larger k)
// instead the incomplete beta integral is employed,
// The terms are not summed directly instead
// the incomplete beta integral is employed,
// according to the formula:
// P = I[1-p]( n-k, k+1).
// = 1 - I[p](k + 1, n - k)
@@ -413,22 +413,16 @@ namespace boost
// and that case has been handled above already.
return 0;
}
if((k < 20) // Small (perhaps up to 34, the largest factorial for float).
&& (floor(k) == k)) // and integral.
{
// For small k, use a finite sum of pdfs:
// it's cheaper than the incomplete beta:
RealType result = 0;
for(unsigned i = 0; i <= k; ++i)
result += pdf(dist, static_cast<RealType>(i));
return result;
}
// else for larger and non-integral k,
// calculate cdf binomial using the incomplete beta function.
//
// P = I[1-p](n - k, k + 1)
// = 1 - I[p](k + 1, n - k)
// Use of ibetac here prevents cancellation errors in calculating
// 1-p if p is very small, perhaps smaller than machine epsilon.
//
// Note that we do not use a finite sum here, since the incomplete
// beta uses a finite sum internally for integer arguments, so
// we'll just let it take care of the necessary logic.
//
return ibetac(k + 1, n - k, p);
} // binomial cdf
@@ -447,8 +441,8 @@ namespace boost
// -- ( i )
// i=k+1
// The terms are not summed directly (at least for larger k)
// instead the incomplete beta integral is employed,
// The terms are not summed directly instead
// the incomplete beta integral is employed,
// according to the formula:
// Q = 1 -I[1-p]( n-k, k+1).
// = I[p](k + 1, n - k)
@@ -496,20 +490,17 @@ namespace boost
// The k == n case has already been handled above.
return 1;
}
if((n - k < 20) && (floor(k) == k) && (floor(n) == n))
{
// For small n-k use a finite sum, it's cheaper
// than the incomplete beta:
RealType result = 0;
for(RealType i = n; i > k; i -= 1)
result += pdf(dist, i);
return result;
}
//
// Calculate cdf binomial using the incomplete beta function.
// Q = 1 -I[1-p](n - k, k + 1)
// = I[p](k + 1, n - k)
// Use of ibeta here prevents cancellation errors in calculating
// 1-p if p is very small, perhaps smaller than machine epsilon.
//
// Note that we do not use a finite sum here, since the incomplete
// beta uses a finite sum internally for integer arguments, so
// we'll just let it take care of the necessary logic.
//
return ibeta(k + 1, n - k, p);
} // binomial cdf

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

@@ -773,6 +773,24 @@ T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const L& l, bool
return sum;
} // template <class T, class L>T beta_small_b_large_a_series(T a, T b, T x, T y, T s0, T mult, const L& l, bool normalised)
//
// For integer arguments we can relate the incomplete beta to the
// complement of the binomial distribution cdf and use this finite sum.
//
template <class T>
T binomial_ccdf(T n, T k, T x, T y)
{
T result = pow(x, n);
T term = result;
for(unsigned i = tools::real_cast<unsigned>(n - 1); i > k; --i)
{
term *= ((i + 1) * y) / ((n - i) * x) ;
result += term;
}
return result;
}
//
// The incomplete beta function implementation:
@@ -948,10 +966,19 @@ T ibeta_imp(T a, T b, T x, const L& l, bool inv, bool normalised)
std::swap(x, y);
invert = !invert;
}
if(b < 40)
{
if(b * x <= 0.7)
if((floor(a) == a) && (floor(b) == b))
{
// relate to the binomial distribution and use a finite sum:
T k = a - 1;
T n = b + k;
fract = binomial_ccdf(n, k, x, y);
if(!normalised)
fract *= boost::math::beta(a, b);
}
else if(b * x <= 0.7)
{
if(!invert)
fract = ibeta_series(a, b, x, T(0), l, normalised);

3
include/boost/math/special_functions/gamma.hpp
View File

@@ -1252,6 +1252,7 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert, const L& l)
using namespace std;
T result;
/*
if((floor(a) == a) && (a < 30) && (a <= x + 1) && (x > 0.6))
{
// calculate Q via finite sum:
@@ -1268,7 +1269,7 @@ T gamma_incomplete_imp(T a, T x, bool normalised, bool invert, const L& l)
if(normalised == false)
result *= gamma_imp(a, l);
}
else if(x < 0.5)
else*/ if(x < 0.5)
{
//
// Changeover criterion chosen to give a changeover at Q ~ 0.33

8
include/boost/math/tools/promotion.hpp
View File

@@ -66,15 +66,15 @@ namespace boost
typedef typename promote_arg<T2>::type T2P; // T2 perhaps promoted.
typedef typename mpl::if_c<
is_floating_point<T1P>::value && is_floating_point<T2P>::value, // both T1P and T2P are floating-point?
typename mpl::if_c<is_same<long double, T1P>::value || is_same<long double, T2P>::value, // either long double?
::boost::is_floating_point<T1P>::value && ::boost::is_floating_point<T2P>::value, // both T1P and T2P are floating-point?
typename mpl::if_c< ::boost::is_same<long double, T1P>::value || ::boost::is_same<long double, T2P>::value, // either long double?
long double, // then result type is long double.
typename mpl::if_c<is_same<double, T1P>::value || is_same<double, T2P>::value, // either double?
typename mpl::if_c< ::boost::is_same<double, T1P>::value || ::boost::is_same<double, T2P>::value, // either double?
double, // result type is double.
float // else result type is float.
>::type
>::type,
typename mpl::if_<is_convertible<T1P, T2P>, T2P, T1P>::type>::type type;
typename mpl::if_< ::boost::is_convertible<T1P, T2P>, T2P, T1P>::type>::type type;
}; // promote_arg2
template <class T1, class T2, class T3>

1005
test/ibeta_int_data.ipp Normal file
View File

File diff suppressed because it is too large Load Diff

30
test/test_ibeta.cpp
View File

@@ -277,6 +277,10 @@ void test_beta(T, const char* name)
# include "ibeta_large_data.ipp"
do_test_beta(ibeta_large_data, name, "Incomplete Beta Function: Large and Diverse Values");
# include "ibeta_int_data.ipp"
do_test_beta(ibeta_int_data, name, "Incomplete Beta Function: Small Integer Values");
}
template <class T>
@@ -389,6 +393,32 @@ void test_spots(T)
static_cast<T>(1),
static_cast<T>(32)/100),
static_cast<T>(0.94856839398626914764591440181367780660208493234722L), tolerance);
// try with some integer arguments:
BOOST_CHECK_CLOSE(
::boost::math::ibeta(
static_cast<T>(3),
static_cast<T>(8),
static_cast<T>(0.25)),
static_cast<T>(0.474407196044921875000000000000000000000000000000000000000000L), tolerance);
BOOST_CHECK_CLOSE(
::boost::math::ibeta(
static_cast<T>(6),
static_cast<T>(8),
static_cast<T>(0.25)),
static_cast<T>(0.0802125930786132812500000000000000000000000000000000000000000L), tolerance);
BOOST_CHECK_CLOSE(
::boost::math::ibeta(
static_cast<T>(12),
static_cast<T>(1),
static_cast<T>(0.25)),
static_cast<T>(5.96046447753906250000000000000000000000000000000000000000000e-8L), tolerance);
BOOST_CHECK_CLOSE(
::boost::math::ibeta(
static_cast<T>(1),
static_cast<T>(8),
static_cast<T>(0.25)),
static_cast<T>(0.899887084960937500000000000000000000000000000000000000000000L), tolerance);
// very naive check on derivative:
using namespace std; // For ADL of std functions

25
tools/ibeta_data.cpp
View File

@@ -4,6 +4,7 @@
// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <boost/math/tools/ntl.hpp>
#include <boost/test/included/test_exec_monitor.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/special_functions/beta.hpp>
#include <boost/math/constants/constants.hpp>
@@ -229,13 +230,24 @@ struct beta_data_generator_small
}
};
struct beta_data_generator_int
{
std::tr1::tuple<NTL::RR, NTL::RR, NTL::RR, NTL::RR, NTL::RR, NTL::RR, NTL::RR> operator()(NTL::RR a, NTL::RR b, NTL::RR x_)
{
float x = truncate_to_float(real_cast<float>(x_));
std::cout << a << " " << b << " " << x << std::endl;
std::pair<NTL::RR, NTL::RR> ib_full = ibeta_fraction1(a, b, NTL::RR(x));
std::pair<NTL::RR, NTL::RR> ib_reg = ibeta_fraction1_regular(a, b, NTL::RR(x));
return std::tr1::make_tuple(a, b, x, ib_full.first, ib_full.second, ib_reg.first, ib_reg.second);
}
};
int main(int argc, char* argv[])
int test_main(int, char* [])
{
NTL::RR::SetPrecision(1000);
NTL::RR::SetOutputPrecision(40);
parameter_info<NTL::RR> arg1, arg2, arg3, arg4;
parameter_info<NTL::RR> arg1, arg2, arg3, arg4, arg5;
test_data<NTL::RR> data;
std::cout << "Welcome.\n"
@@ -247,17 +259,20 @@ int main(int argc, char* argv[])
arg2 = make_periodic_param(NTL::RR(-5), NTL::RR(6), 11);
arg3 = make_random_param(NTL::RR(0.0001), NTL::RR(1), 10);
arg4 = make_random_param(NTL::RR(0.0001), NTL::RR(1), 100 /*500*/);
arg5 = make_periodic_param(NTL::RR(1), NTL::RR(41), 10);
arg1.type |= dummy_param;
arg2.type |= dummy_param;
arg3.type |= dummy_param;
arg4.type |= dummy_param;
arg5.type |= dummy_param;
// comment out all but one of the following when running
// or this program will take forever to complete!
data.insert(beta_data_generator(), arg1, arg2, arg3);
data.insert(beta_data_generator_medium(), arg4);
data.insert(beta_data_generator_small(), arg4);
//data.insert(beta_data_generator(), arg1, arg2, arg3);
//data.insert(beta_data_generator_medium(), arg4);
//data.insert(beta_data_generator_small(), arg4);
data.insert(beta_data_generator_int(), arg5, arg5, arg3);
test_data<NTL::RR>::const_iterator i, j;
i = data.begin();