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

Tightened up root finding tests.

Browse Source 
Added docs for roots without derivatives.
Added derivatives of incomplete gamma and beta.


[SVN r3112]
This commit is contained in:
John Maddock
2006-08-01 17:23:10 +00:00
parent cd90a9f0d4
commit 46d9b33645
12 changed files with 678 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files

44
doc/beta_gamma_derivatives.qbk Normal file
View File

@@ -0,0 +1,44 @@
[#beta_gamma_derivatives][section Derivatives of the Incomplete Beta and Gamma Functions]
[caution __caution ]
[h4 Synopsis]
``
#include <boost/math/tools/roots.hpp>
``
namespace boost{ namespace math{
template <class T>
T gamma_P_derivative(T a, T x);
template <class T>
T ibeta_derivative(T a, T b, T x);
}} // namespaces
[h4 Description]
These two functions find some uses in statistical distributions, they
implement the partial derivative with respect to /x/ of the incomplete
gamma and beta functions respectively.
[$../equations/derivative1.png]
[$../equations/derivative2.png]
[h4 Accuracy]
Almost identical to the incomplete gamma and beta functions, refer to
the documentation for those function for more information.
[h4 Implementation]
These two functions just expose some of the internals of the incomplete
gamma and beta functions, refer to the documentation for those functions
for more information.
[endsect]

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

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 1.1 KiB

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

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 1.2 KiB

145
doc/equations/derivatives.xml Normal file
View File

@@ -0,0 +1,145 @@
<?xml version='1.0'?>
<!DOCTYPE html PUBLIC '-//W3C//DTD XHTML 1.1 plus MathML 2.0//EN'
'http://www.w3.org/TR/MathML2/dtd/xhtml-math11-f.dtd'
[<!ENTITY mathml 'http://www.w3.org/1998/Math/MathML'>]>
<html xmlns='http://www.w3.org/1999/xhtml'>
<head>
<!-- MathML created with MathCast Equation Editor version 0.83 -->
</head>
<body>
<math xmlns="http://www.w3.org/1998/Math/MathML" display="block">
<mrow>
<mtext>gamma_P_derivative</mtext>
<mfenced>
<mrow>
<mi>a</mi>
<mo>,</mo>
<mi>x</mi>
</mrow>
</mfenced>
<mspace width="1em"/>
<mo>=</mo>
<mspace width="1em"/>
<mfrac>
<mo>&PartialD;</mo>
<mrow>
<mo>&PartialD;</mo>
<mi>x</mi>
</mrow>
</mfrac>
<mi>P</mi>
<mfenced>
<mrow>
<mi>a</mi>
<mo>,</mo>
<mi>x</mi>
</mrow>
</mfenced>
<mspace width="1em"/>
<mo>=</mo>
<mspace width="1em"/>
<mfrac>
<mrow>
<msup>
<mi>e</mi>
<mrow>
<mo>&minus;</mo>
<mi>x</mi>
</mrow>
</msup>
<msup>
<mi>x</mi>
<mrow>
<mi>a</mi>
<mo>&minus;</mo>
<mn>1</mn>
</mrow>
</msup>
</mrow>
<mrow>
<mi>&Gamma;</mi>
<mfenced>
<mrow>
<mi>a</mi>
</mrow>
</mfenced>
</mrow>
</mfrac>
</mrow>
</math>
<math xmlns="http://www.w3.org/1998/Math/MathML" display="block">
<mrow>
<mtext>ibeta_derivative</mtext>
<mfenced>
<mrow>
<mi>a</mi>
<mo>,</mo>
<mi>b</mi>
<mo>,</mo>
<mi>x</mi>
</mrow>
</mfenced>
<mspace width="1em"/>
<mo>=</mo>
<mspace width="1em"/>
<mfrac>
<mo>&PartialD;</mo>
<mrow>
<mo>&PartialD;</mo>
<mi>x</mi>
</mrow>
</mfrac>
<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"/>
<mfrac>
<mrow>
<msup>
<mfenced>
<mrow>
<mn>1</mn>
<mo>&minus;</mo>
<mi>x</mi>
</mrow>
</mfenced>
<mrow>
<mi>b</mi>
<mo>&minus;</mo>
<mn>1</mn>
</mrow>
</msup>
<msup>
<mi>x</mi>
<mrow>
<mi>a</mi>
<mo>&minus;</mo>
<mn>1</mn>
</mrow>
</msup>
</mrow>
<mrow>
<mi>B</mi>
<mfenced>
<mrow>
<mi>a</mi>
<mo>,</mo>
<mi>b</mi>
</mrow>
</mfenced>
</mrow>
</mfrac>
</mrow>
</math>
</body>
</html>

1
doc/math.qbk
View File

@@ -81,6 +81,7 @@ data and/or data for output to an external graphing application.
[include beta.qbk]
[include ibeta.qbk]
[include ibeta_inv.qbk]
[include beta_gamma_derivatives.qbk]
[include fpclassify.qbk]
[include error_handling.qbk]

310
doc/roots_without_derivatives.qbk Normal file
View File

@@ -0,0 +1,310 @@
[#roots2][section Root Finding Without Derivatives]
[caution __caution ]
[h4 Synopsis]
``
#include <boost/math/tools/roots.hpp>
``
namespace boost{ namespace math{ namespace tools{
template <class F, class T, class Tol>
std::pair<T, T>
bisect(
F f,
T min,
T max,
Tol tol,
boost::uintmax_t& max_iter);
template <class F, class T, class Tol>
std::pair<T, T>
bisect(
F f,
T min,
T max,
Tol tol);
template <class F, class T, class Tol>
std::pair<T, T>
bracket_and_solve_root(
F f,
const T& guess,
const T& factor,
bool rising,
Tol tol,
boost::uintmax_t& max_iter);
template <class F, class T, class Tol>
std::pair<T, T>
toms748_solve(
F f,
const T& a,
const T& b,
Tol tol,
boost::uintmax_t& max_iter);
template <class F, class T, class Tol>
std::pair<T, T>
toms748_solve(
F f,
const T& a,
const T& b,
const T& fa,
const T& fb,
Tol tol,
boost::uintmax_t& max_iter);
// Termination conditions:
template <class T>
struct eps_tolerance;
struct equal_floor;
struct equal_ceil;
}}} // namespaces
[h4 Description]
These functions solve the root of some function /f(x)/ without the
need for the derivatives of /f(x)/. The functions here that use TOMS
Algorithm 748 are asymptotically the most efficient known, and have
been shown to be optimal for a certain classes of smooth functions.
Alternatively there is a simple bisection routine which can be useful
in it's own right in some situations, or alternatively for narrowing
down the range containing the root prior to calling a more advanced
algorithm.
All the algorithms in this section reduce the diameter of the enclosing
interval with the same asymptotic efficiency with which they locate the
root. This is in contrast to the derivative based methods which may /never/
significantly reduce the enclosing interval, even though they rapidly approach
the root. This is also in contrast to some other derivate free methods
(for example the methods of Brent or Dekker) which only reduce the enclosing
interval on the final step. Therefore these methods return a std::pair
containing the enclosing
interval found, and accept a function object specifying the termination condition.
Three function objects are provided for ready-made termination conditions:
/eps_tolerance/ causes termination when the relative error in the enclosing
interval is below a certain threshold, while /equal_floor/ and /equal_ceil/ are
useful for certain statistical applications where the result is known to be
an integer. Other user defined termination conditions are likely to be used
only rairly, but may be useful in some specific circumstances.
template <class F, class T, class Tol>
std::pair<T, T>
bisect(
F f,
T min,
T max,
Tol tol,
boost::uintmax_t& max_iter);
template <class F, class T, class Tol>
std::pair<T, T>
bisect(
F f,
T min,
T max,
Tol tol);
These two functions locate the root using bisection, function arguments are:
[variablelist
[[f] [A unary functor which the function whose root is to be found.]]
[[min] [The left bracket of the interval known to contain the root.]]
[[max] [The right bracket of the interval known to contain the root.
It is a precondition that /min < max/ and /f(min)*f(max) <= 0/,
the function calls __logic_error if these preconditions are violated.]]
[[tol] [A binary functor that specifies the termination condition: the function
will return the current brackets enclosing the root when /tol(min,max)/ becomes true.]]
[[max_iter][The maximum number of invocations of /f(x)/ to make while searching for the root.]]
]
Returns: a pair of values /r/ that bracket the root so that:
f(r.first) * f(r.second) <= 0
and either
tol(r.first, r.second) == true
or
max_iter >= m
where /m/ is the initial value of /max_iter/ passed to the function.
In other words it's up to the caller to verify whether termination occured
as a result of exceeding /max_iter/ function invocations (easily done by
checking the value of /max_iter/), rather than because the termination
condition /tol/ was satisfied.
template <class F, class T, class Tol>
std::pair<T, T>
bracket_and_solve_root(
F f,
const T& guess,
const T& factor,
bool rising,
Tol tol,
boost::uintmax_t& max_iter);
This is a convenience function that calls /toms748_solve/ internally
to find the root of /f(x)/. It's usable only when /f(x)/ is a monotonic
function, and the location of the root is known approximately, and in
particular it is known whether the root is occurs for positive or negative
/x/. The parameters are:
[variablelist
[[f][A unary functor that is the function whose root is to be solved.
f(x) must be uniformly increasing or decreasing on /x/.]]
[[guess][An initial approximation to the root]]
[[factor][A scaling factor that is used to bracket the root: the value
/guess/ is mutiplied (or divided as appropriate) by /factor/
until two values are found that bracket the root. A value
such as 2 is a typical choice for /factor/.]]
[[rising][Set to /true/ if /f(x)/ is rising on /x/ and /false/ if /f(x)/
is falling on /x/. This value is used along with the result
of /f(guess)/ to determine if /guess/ is
above or below the root.]]
[[tol] [A binary functor that determines the termination condition for the search
for the root. /tol/ is passed the current brackets at each step,
when it returns true then the current brackets are returned as the result.]]
[[max_iter] [The maximum number of function invocations to perform in the search
for the root.]]
]
Returns: a pair of values /r/ that bracket the root so that:
f(r.first) * f(r.second) <= 0
and either
tol(r.first, r.second) == true
or
max_iter >= m
where /m/ is the initial value of /max_iter/ passed to the function.
In other words it's up to the caller to verify whether termination occured
as a result of exceeding /max_iter/ function invocations (easily done by
checking the value of /max_iter/), rather than because the termination
condition /tol/ was satisfied.
template <class F, class T, class Tol>
std::pair<T, T>
toms748_solve(
F f,
const T& a,
const T& b,
Tol tol,
boost::uintmax_t& max_iter);
template <class F, class T, class Tol>
std::pair<T, T>
toms748_solve(
F f,
const T& a,
const T& b,
const T& fa,
const T& fb,
Tol tol,
boost::uintmax_t& max_iter);
These two functions implement TOMS Algorithm 748: it uses a mixture of
cubic, quadratic and linear (secant) interpolation to locate the root of
/f(x)/. The two functions differ only by whether values for /f(a)/ and
/f(b)/ are already available. The parameters are:
[variablelist
[[f] [A unary functor that is the function whose root is to be solved.
f(x) need not be uniformly increasing or decreasing on /x/ and
may have mutilple roots.]]
[[a] [ The lower bound for the initial bracket of the root.]]
[[b] [The upper bound for the initial bracket of the root.
It is a precondition that /a < b/ and that /a/ and /b/
bracket the root to find so that /f(a)*f(b) < 0/.]]
[[fa] [Optional: the value of /f(a)/.]]
[[fb] [Optional: the value of /f(b)/.]]
[[tol] [A binary functor that determines the termination condition for the search
for the root. /tol/ is passed the current brackets at each step,
when it returns true then the current brackets are returned as the result.]]
[[max_iter] [The maximum number of function invocations to perform in the search
for the root. On exit /max_iter/ is set to actual number of function
invocations used.]]
]
Returns: a pair of values /r/ that bracket the root so that:
f(r.first) * f(r.second) <= 0
and either
tol(r.first, r.second) == true
or
max_iter >= m
where /m/ is the initial value of /max_iter/ passed to the function.
In other words it's up to the caller to verify whether termination occured
as a result of exceeding /max_iter/ function invocations (easily done by
checking the value of /max_iter/), rather than because the termination
condition /tol/ was satisfied.
template <class T>
struct eps_tolerance
{
eps_tolerance(int bits);
bool operator()(const T& a, const T& b)const;
};
This is the usual termination condition used with these root finding functions.
It's operator() will return true when the relative distance between /a/ and /b/
is less than twice the machine epsilon for T, or 2[super 1-bits], whichever is
the larger. In other words you set /bits/ to the number of bits of precision you
want in the result. The minimal tolerance of twice the machine epsilon of T is
required to ensure that we get back a bracketing interval: since this must clearly
be at least 1 epsilon in size.
struct equal_floor
{
equal_floor();
template <class T> bool operator()(const T& a, const T& b)const;
};
This termination condition is used when you want to find an integer result
that is the /floor/ of the true root. It will terminate as soon as both ends
of the interval have the same /floor/.
struct equal_ceil
{
equal_ceil();
template <class T> bool operator()(const T& a, const T& b)const;
};
This termination condition is used when you want to find an integer result
that is the /ceil/ of the true root. It will terminate as soon as both ends
of the interval have the same /ceil/.
[h4 Implementation]
The implementation of the bisection algorithm is extrememly straightforward
and not detailed here. TOMS algorithm 748 is described in detail in:
['Algorithm 748: Enclosing Zeros of Continuous Functions,
G. E. Alefeld, F. A. Potra and Yixun Shi,
ACM Transactions on Mathematica1 Software, Vol. 21. No. 3. September 1995.
Pages 327-344.]
The implementation here is a faithful translation of this paper into C++.
[endsect]

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

@@ -926,6 +926,51 @@ T ibeta_imp(T a, T b, T x, const L& l, bool inv, bool normalised)
return invert ? (normalised ? 1 : beta_imp(a, b, l)) - fract : fract;
} // template <class T, class L>T ibeta_imp(T a, T b, T x, const L& l, bool inv, bool normalised)
template <class T, class L>
T ibeta_derivative_imp(T a, T b, T x, const L& l)
{
//
// start with the usual error checks:
//
if(a <= 0)
tools::domain_error<T>(BOOST_CURRENT_FUNCTION, "The argument a to the incomplete beta function must be greater than zero (got a=%1%).", a);
if(b <= 0)
tools::domain_error<T>(BOOST_CURRENT_FUNCTION, "The argument b to the incomplete beta function must be greater than zero (got b=%1%).", b);
if((x < 0) || (x > 1))
tools::domain_error<T>(BOOST_CURRENT_FUNCTION, "Parameter x outside the range [0,1] in the incomplete beta function (got x=%1%).", x);
//
// Now the corner cases:
//
if(x == 0)
{
return (a > 1) ? 0 :
(a == 1) ? 1 : tools::overflow_error<T>(BOOST_CURRENT_FUNCTION, 0);
}
else if(x == 1)
{
return (b > 1) ? 0 :
(b == 1) ? 1 : tools::overflow_error<T>(BOOST_CURRENT_FUNCTION, 0);
}
//
// Now the regular cases:
//
T f1 = ibeta_power_terms(a, b, x, 1 - x, l, true);
T y = (1 - x) * x;
if(f1 == 0)
return 0;
if((tools::max_value<T>() * y < f1))
{
// overflow:
return tools::overflow_error<T>(BOOST_CURRENT_FUNCTION, 0);
}
f1 /= y;
return f1;
}
} // namespace detail
//
@@ -973,6 +1018,14 @@ T ibetac(T a, T b, T x)
return tools::checked_narrowing_cast<T>(detail::ibeta_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), evaluation_type(), true, true), BOOST_CURRENT_FUNCTION);
}
template <class T>
T ibeta_derivative(T a, T b, T x)
{
typedef typename lanczos::lanczos_traits<T>::value_type value_type;
typedef typename lanczos::lanczos_traits<T>::evaluation_type evaluation_type;
return tools::checked_narrowing_cast<T>(detail::ibeta_derivative_imp(static_cast<value_type>(a), static_cast<value_type>(b), static_cast<value_type>(x), evaluation_type()), BOOST_CURRENT_FUNCTION);
}
} // namespace math
} // namespace boost

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

@@ -873,6 +873,39 @@ T tgamma_delta_ratio_imp(T z, T delta, const lanczos::undefined_lanczos&)
return sum * prefix;
}
template <class T, class L>
T gamma_P_derivative_imp(T a, T x, const L& l)
{
//
// Usual error checks first:
//
if(a <= 0)
tools::domain_error<T>(BOOST_CURRENT_FUNCTION, "Argument a to the incomplete gamma function must be greater than zero (got a=%1%).", a);
if(x < 0)
tools::domain_error<T>(BOOST_CURRENT_FUNCTION, "Argument x to the incomplete gamma function must be >= 0 (got x=%1%).", x);
//
// Now special cases:
//
if(x == 0)
{
return (a > 1) ? 0 :
(a == 1) ? 1 : tools::overflow_error<T>(BOOST_CURRENT_FUNCTION, 0);
}
//
// Normal case:
//
T f1 = detail::regularised_gamma_prefix(a, x, l);
if((x < 1) && (tools::max_value<T>() * x < f1))
{
// overflow:
return tools::overflow_error<T>(BOOST_CURRENT_FUNCTION, 0);
}
f1 /= x;
return f1;
}
} // namespace detail
template <class T>
@@ -974,6 +1007,14 @@ inline T tgamma_ratio(T a, T b)
return tools::checked_narrowing_cast<T>(detail::tgamma_delta_ratio_imp(static_cast<value_type>(a), static_cast<value_type>(b - a), evaluation_type()), BOOST_CURRENT_FUNCTION);
}
template <class T>
T gamma_P_derivative(T a, T x)
{
typedef typename lanczos::lanczos_traits<T>::value_type value_type;
typedef typename lanczos::lanczos_traits<T>::evaluation_type evaluation_type;
return tools::checked_narrowing_cast<T>(detail::gamma_P_derivative_imp(static_cast<value_type>(a), static_cast<value_type>(x), evaluation_type()), BOOST_CURRENT_FUNCTION);
}
} // namespace math
} // namespace boost

10
test/test_ibeta.cpp
View File

@@ -319,6 +319,16 @@ void test_spots(T)
static_cast<T>(1),
static_cast<T>(0.32)),
static_cast<T>(0.948565954109602496577407403168592262389L), tolerance);
// very naive check on derivative:
using namespace std; // For ADL of std functions
tolerance = boost::math::tools::epsilon<T>() * 10000; // 100 eps
BOOST_CHECK_CLOSE(
::boost::math::ibeta_derivative(
static_cast<T>(2),
static_cast<T>(3),
static_cast<T>(0.5)),
pow(static_cast<T>(0.5), static_cast<T>(2)) * pow(static_cast<T>(0.5), static_cast<T>(1)) / boost::math::beta(static_cast<T>(2), static_cast<T>(3)), tolerance);
}
int test_main(int, char* [])

7
test/test_igamma.cpp
View File

@@ -275,6 +275,13 @@ void test_spots(T)
BOOST_CHECK_CLOSE(::boost::math::gamma_P(static_cast<T>(5), static_cast<T>(100)), static_cast<T>(0.9999999999999999999999999999999999998386069466302L), tolerance);
BOOST_CHECK_CLOSE(::boost::math::gamma_P(static_cast<T>(1.5), static_cast<T>(2)), static_cast<T>(0.73853587005088937779717792402407879809718939080921L), tolerance);
BOOST_CHECK_CLOSE(::boost::math::gamma_P(static_cast<T>(20.5), static_cast<T>(22)), static_cast<T>(0.65424667956532673185028409120341593367429721070928L), tolerance);
// naive check on derivative function:
using namespace std; // For ADL of std functions
tolerance = boost::math::tools::epsilon<T>() * 5000; // 50 eps
BOOST_CHECK_CLOSE(::boost::math::gamma_P_derivative(static_cast<T>(20.5), static_cast<T>(22)),
exp(static_cast<T>(-22)) * pow(static_cast<T>(22), static_cast<T>(19.5)) / boost::math::tgamma(static_cast<T>(20.5)), tolerance);
}
int test_main(int, char* [])

5
test/test_toms748_solve.cpp
View File

@@ -120,6 +120,7 @@ void run_test(T a, T b, int id)
b,
boost::math::tools::eps_tolerance<double>(std::numeric_limits<double>::digits),
c);
BOOST_CHECK_EQUAL(c, toms748tester<double>::total_calls());
total += c;
invocations += toms748tester<double>::total_calls();
std::cout << "Function " << id << "\nresult={" << r.first << ", " << r.second << "} total calls=" << toms748tester<double>::total_calls() << "\n\n";
@@ -135,6 +136,7 @@ void run_test(T a, T b, int id, int p)
b,
boost::math::tools::eps_tolerance<double>(std::numeric_limits<double>::digits),
c);
BOOST_CHECK_EQUAL(c, toms748tester<double>::total_calls());
total += c;
invocations += toms748tester<double>::total_calls();
std::cout << "Function " << id << "\nresult={" << r.first << ", " << r.second << "} total calls=" << toms748tester<double>::total_calls() << "\n\n";
@@ -150,6 +152,7 @@ void run_test(T a, T b, int id, T p1, T p2)
b,
boost::math::tools::eps_tolerance<double>(std::numeric_limits<double>::digits),
c);
BOOST_CHECK_EQUAL(c, toms748tester<double>::total_calls());
total += c;
invocations += toms748tester<double>::total_calls();
std::cout << "Function " << id << "\n Result={" << r.first << ", " << r.second << "} total calls=" << toms748tester<double>::total_calls() << "\n\n";
@@ -259,7 +262,7 @@ int test_main(int, char* [])
std::cout << std::setprecision(18);
std::cout << "Function " << 4 << "\n Result={" << r.first << ", " << r.second << "} total calls=" << toms748tester<double>::total_calls() << "\n\n";
toms748tester<double>::reset();
BOOST_CHECK(c < 5);
BOOST_CHECK(c < 20);
}
return 0;

63
tools/log1p_expm1_data.cpp Normal file
View File

@@ -0,0 +1,63 @@
// (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/tools/ntl.hpp>
#include <fstream>
#include <boost/math/tools/test_data.hpp>
#include <boost/math/special_functions/log1p.hpp>
#include <boost/math/special_functions/expm1.hpp>
using namespace boost::math::tools;
using namespace std;
struct data_generator
{
std::tr1::tuple<NTL::RR, NTL::RR> operator()(NTL::RR z)
{
return std::tr1::make_tuple(boost::math::log1p(z), boost::math::expm1(z));
}
};
int main(int argc, char* argv[])
{
NTL::RR::SetPrecision(1000);
NTL::RR::SetOutputPrecision(40);
parameter_info<NTL::RR> arg1;
test_data<NTL::RR> data;
std::cout << "Welcome.\n"
"This program will generate spot tests for the log1p and expm1 functions:\n\n";
bool cont;
std::string line;
do{
if(0 == get_user_parameter_info(arg1, "z"))
return 1;
data.insert(data_generator(), arg1);
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=log1p_expm1_data.ipp]";
std::getline(std::cin, line);
boost::algorithm::trim(line);
if(line == "")
line = "log1p_expm1_data.ipp";
std::ofstream ofs(line.c_str());
ofs << std::scientific;
write_code(ofs, data, "log1p_expm1_data");
return 0;
}