boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-02-26 16:52:27 +00:00
Code Issues Projects Releases Wiki Activity

Merge branch 'develop'

Browse Source
This commit is contained in:
John Maddock
2023-03-07 17:59:58 +00:00
parent 310fc00f67 cf6cd75477
commit c56f334348
14 changed files with 236 additions and 67 deletions
Download Patch File Download Diff File Expand all files Collapse all files

7
include/boost/math/distributions/detail/common_error_handling.hpp
View File

@@ -185,11 +185,12 @@ inline bool check_non_centrality(
RealType* result,
const Policy& pol)
{
if((ncp < 0) || !(boost::math::isfinite)(ncp))
{ // Assume scale == 0 is NOT valid for any distribution.
static const RealType upper_limit = static_cast<RealType>((std::numeric_limits<long long>::max)()) - boost::math::policies::get_max_root_iterations<Policy>();
if((ncp < 0) || !(boost::math::isfinite)(ncp) || ncp > upper_limit)
{
*result = policies::raise_domain_error<RealType>(
function,
"Non centrality parameter is %1%, but must be > 0 !", ncp, pol);
"Non centrality parameter is %1%, but must be > 0, and a countable value such that x+1 != x", ncp, pol);
return false;
}
return true;

10
include/boost/math/distributions/non_central_beta.hpp
View File

@@ -430,7 +430,7 @@ namespace boost
static_cast<value_type>(p),
&r,
Policy()))
return (RealType)r;
return static_cast<RealType>(r);
//
// Special cases first:
//
@@ -624,7 +624,7 @@ namespace boost
static_cast<value_type>(x),
&r,
Policy()))
return (RealType)r;
return static_cast<RealType>(r);
if(l == 0)
return pdf(boost::math::beta_distribution<RealType, Policy>(dist.alpha(), dist.beta()), x);
@@ -761,7 +761,7 @@ namespace boost
l,
&r,
Policy()))
return (RealType)r;
return static_cast<RealType>(r);
RealType c = a + b + l / 2;
RealType mean = 1 - (b / c) * (1 + l / (2 * c * c));
return detail::generic_find_mode_01(
@@ -872,7 +872,7 @@ namespace boost
x,
&r,
Policy()))
return (RealType)r;
return static_cast<RealType>(r);
if(l == 0)
return cdf(beta_distribution<RealType, Policy>(a, b), x);
@@ -909,7 +909,7 @@ namespace boost
x,
&r,
Policy()))
return (RealType)r;
return static_cast<RealType>(r);
if(l == 0)
return cdf(complement(beta_distribution<RealType, Policy>(a, b), x));

24
include/boost/math/distributions/non_central_chi_squared.hpp
View File

@@ -88,7 +88,7 @@ namespace boost
// stable direction for the gamma function
// recurrences:
//
int i;
long long i;
for(i = k; static_cast<std::uintmax_t>(i-k) < max_iter; ++i)
{
T term = poisf * gamf;
@@ -299,7 +299,7 @@ namespace boost
if(pois == 0)
return 0;
T poisb = pois;
for(int i = k; ; ++i)
for(long long i = k; ; ++i)
{
sum += pois;
if(pois / sum < errtol)
@@ -310,7 +310,7 @@ namespace boost
"Series did not converge, closest value was %1%", sum, pol);
pois *= l2 * x2 / ((i + 1) * (n2 + i));
}
for(int i = k - 1; i >= 0; --i)
for(long long i = k - 1; i >= 0; --i)
{
poisb *= (i + 1) * (n2 + i) / (l2 * x2);
sum += poisb;
@@ -428,7 +428,7 @@ namespace boost
static_cast<value_type>(p),
&r,
Policy()))
return (RealType)r;
return static_cast<RealType>(r);
//
// Special cases get short-circuited first:
//
@@ -519,7 +519,7 @@ namespace boost
(value_type)x,
&r,
Policy()))
return (RealType)r;
return static_cast<RealType>(r);
if(l == 0)
return pdf(boost::math::chi_squared_distribution<RealType, forwarding_policy>(dist.degrees_of_freedom()), x);
@@ -821,7 +821,7 @@ namespace boost
l,
&r,
Policy()))
return r;
return static_cast<RealType>(r);
return k + l;
} // mean
@@ -842,7 +842,7 @@ namespace boost
l,
&r,
Policy()))
return (RealType)r;
return static_cast<RealType>(r);
bool asymptotic_mode = k < l/4;
RealType starting_point = asymptotic_mode ? k + l - RealType(3) : RealType(1) + k;
return detail::generic_find_mode(dist, starting_point, function);
@@ -864,7 +864,7 @@ namespace boost
l,
&r,
Policy()))
return r;
return static_cast<RealType>(r);
return 2 * (2 * l + k);
}
@@ -887,7 +887,7 @@ namespace boost
l,
&r,
Policy()))
return r;
return static_cast<RealType>(r);
BOOST_MATH_STD_USING
return pow(2 / (k + 2 * l), RealType(3)/2) * (k + 3 * l);
}
@@ -908,7 +908,7 @@ namespace boost
l,
&r,
Policy()))
return r;
return static_cast<RealType>(r);
return 12 * (k + 4 * l) / ((k + 2 * l) * (k + 2 * l));
} // kurtosis_excess
@@ -946,7 +946,7 @@ namespace boost
x,
&r,
Policy()))
return r;
return static_cast<RealType>(r);
return detail::non_central_chi_squared_cdf(x, k, l, false, Policy());
} // cdf
@@ -975,7 +975,7 @@ namespace boost
x,
&r,
Policy()))
return r;
return static_cast<RealType>(r);
return detail::non_central_chi_squared_cdf(x, k, l, true, Policy());
} // ccdf

14
include/boost/math/distributions/non_central_f.hpp
View File

@@ -106,7 +106,7 @@ namespace boost
l,
&r,
Policy()))
return r;
return r;
if(v2 <= 2)
return policies::raise_domain_error(
function,
@@ -137,7 +137,7 @@ namespace boost
l,
&r,
Policy()))
return r;
return r;
RealType guess = m > 2 ? RealType(m * (n + l) / (n * (m - 2))) : RealType(1);
return detail::generic_find_mode(
dist,
@@ -166,7 +166,7 @@ namespace boost
l,
&r,
Policy()))
return r;
return r;
if(m <= 4)
return policies::raise_domain_error(
function,
@@ -203,7 +203,7 @@ namespace boost
l,
&r,
Policy()))
return r;
return r;
if(m <= 6)
return policies::raise_domain_error(
function,
@@ -240,7 +240,7 @@ namespace boost
l,
&r,
Policy()))
return r;
return r;
if(m <= 8)
return policies::raise_domain_error(
function,
@@ -309,7 +309,7 @@ namespace boost
dist.non_centrality(),
&r,
Policy()))
return r;
return r;
if((x < 0) || !(boost::math::isfinite)(x))
{
@@ -350,7 +350,7 @@ namespace boost
c.dist.non_centrality(),
&r,
Policy()))
return r;
return r;
if((c.param < 0) || !(boost::math::isfinite)(c.param))
{

58
include/boost/math/distributions/non_central_t.hpp
View File

@@ -307,9 +307,9 @@ namespace boost
function,
v, &r, Policy())
||
!detail::check_finite(
!detail::check_non_centrality(
function,
delta,
static_cast<T>(delta * delta),
&r,
Policy())
||
@@ -730,9 +730,9 @@ namespace boost
detail::check_df_gt0_to_inf(
function,
v, &r, Policy());
detail::check_finite(
detail::check_non_centrality(
function,
lambda,
static_cast<value_type>(lambda * lambda),
&r,
Policy());
} // non_central_t_distribution constructor.
@@ -874,12 +874,12 @@ namespace boost
function,
v, &r, Policy())
||
!detail::check_finite(
!detail::check_non_centrality(
function,
l,
static_cast<RealType>(l * l),
&r,
Policy()))
return (RealType)r;
return static_cast<RealType>(r);
BOOST_MATH_STD_USING
@@ -912,12 +912,12 @@ namespace boost
function,
v, &r, Policy())
||
!detail::check_finite(
!detail::check_non_centrality(
function,
l,
static_cast<RealType>(l * l),
&r,
Policy()))
return (RealType)r;
return static_cast<RealType>(r);
if(v <= 1)
return policies::raise_domain_error<RealType>(
function,
@@ -947,12 +947,12 @@ namespace boost
function,
v, &r, Policy())
||
!detail::check_finite(
!detail::check_non_centrality(
function,
l,
static_cast<RealType>(l * l),
&r,
Policy()))
return (RealType)r;
return static_cast<RealType>(r);
if(v <= 2)
return policies::raise_domain_error<RealType>(
function,
@@ -982,12 +982,12 @@ namespace boost
function,
v, &r, Policy())
||
!detail::check_finite(
!detail::check_non_centrality(
function,
l,
static_cast<RealType>(l * l),
&r,
Policy()))
return (RealType)r;
return static_cast<RealType>(r);
if(v <= 3)
return policies::raise_domain_error<RealType>(
function,
@@ -1014,12 +1014,12 @@ namespace boost
function,
v, &r, Policy())
||
!detail::check_finite(
!detail::check_non_centrality(
function,
l,
static_cast<RealType>(l * l),
&r,
Policy()))
return (RealType)r;
return static_cast<RealType>(r);
if(v <= 4)
return policies::raise_domain_error<RealType>(
function,
@@ -1053,9 +1053,9 @@ namespace boost
function,
v, &r, Policy())
||
!detail::check_finite(
!detail::check_non_centrality(
function,
l,
static_cast<RealType>(l * l), // we need l^2 to be countable.
&r,
Policy())
||
@@ -1064,7 +1064,7 @@ namespace boost
t,
&r,
Policy()))
return (RealType)r;
return static_cast<RealType>(r);
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
detail::non_central_t_pdf(static_cast<value_type>(v),
static_cast<value_type>(l),
@@ -1093,9 +1093,9 @@ namespace boost
function,
v, &r, Policy())
||
!detail::check_finite(
!detail::check_non_centrality(
function,
l,
static_cast<RealType>(l * l),
&r,
Policy())
||
@@ -1104,8 +1104,8 @@ namespace boost
x,
&r,
Policy()))
return (RealType)r;
if ((boost::math::isinf)(v))
return static_cast<RealType>(r);
if ((boost::math::isinf)(v))
{ // Infinite degrees of freedom, so use normal distribution located at delta.
normal_distribution<RealType, Policy> n(l, 1);
cdf(n, x);
@@ -1147,9 +1147,9 @@ namespace boost
function,
v, &r, Policy())
||
!detail::check_finite(
!detail::check_non_centrality(
function,
l,
static_cast<RealType>(l * l),
&r,
Policy())
||
@@ -1158,7 +1158,7 @@ namespace boost
x,
&r,
Policy()))
return (RealType)r;
return static_cast<RealType>(r);
if ((boost::math::isinf)(v))
{ // Infinite degrees of freedom, so use normal distribution located at delta.

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

@@ -118,9 +118,17 @@ T temme_method_1_ibeta_inverse(T a, T b, T z, const Policy& pol)
x = 0.5;
else
x = (1 + eta * sqrt((1 + c) / eta_2)) / 2;
BOOST_MATH_ASSERT(x >= 0);
BOOST_MATH_ASSERT(x <= 1);
//
// These are post-conditions of the method, but the addition above
// may result in us being out by 1ulp either side of the boundary,
// so just check that we're in bounds and adjust as needed.
// See https://github.com/boostorg/math/issues/961
//
if (x < 0)
x = 0;
else if (x > 1)
x = 1;
BOOST_MATH_ASSERT(eta * (x - 0.5) >= 0);
#ifdef BOOST_INSTRUMENT
std::cout << "Estimating x with Temme method 1: " << x << std::endl;
@@ -901,6 +909,8 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
if(x < lower)
x = lower;
}
std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
std::uintmax_t max_iter_used = 0;
//
// Figure out how many digits to iterate towards:
//
@@ -923,10 +933,9 @@ T ibeta_inv_imp(T a, T b, T p, T q, const Policy& pol, T* py)
// Now iterate, we can use either p or q as the target here
// depending on which is smaller:
//
std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
x = boost::math::tools::halley_iterate(
boost::math::detail::ibeta_roots<T, Policy>(a, b, (p < q ? p : q), (p < q ? false : true)), x, lower, upper, digits, max_iter);
policies::check_root_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%)", max_iter, pol);
policies::check_root_iterations<T>("boost::math::ibeta<%1%>(%1%, %1%, %1%)", max_iter + max_iter_used, pol);
//
// We don't really want these asserts here, but they are useful for sanity
// checking that we have the limits right, uncomment if you suspect bugs *only*.

24
include/boost/math/tools/roots.hpp
View File

@@ -586,8 +586,8 @@ namespace detail {
// we can jump way out of bounds if we're not careful.
// See https://svn.boost.org/trac/boost/ticket/8314.
delta = f0 / f1;
if (fabs(delta) > 2 * fabs(guess))
delta = (delta < 0 ? -1 : 1) * 2 * fabs(guess);
if (fabs(delta) > 2 * fabs(result))
delta = (delta < 0 ? -1 : 1) * 2 * fabs(result);
}
}
else
@@ -600,9 +600,19 @@ namespace detail {
if ((convergence > 0.8) && (convergence < 2))
{
// last two steps haven't converged.
delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
if ((result != 0) && (fabs(delta) > result))
delta = sign(delta) * fabs(result) * 0.9f; // protect against huge jumps!
if (fabs(min) < 1 ? fabs(1000 * min) < fabs(max) : fabs(max / min) > 1000)
{
if(delta > 0)
delta = bracket_root_towards_min(f, result, f0, min, max, count);
else
delta = bracket_root_towards_max(f, result, f0, min, max, count);
}
else
{
delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
if ((result != 0) && (fabs(delta) > result))
delta = sign(delta) * fabs(result) * 0.9f; // protect against huge jumps!
}
// reset delta2 so that this branch will *not* be taken on the
// next iteration:
delta2 = delta * 3;
@@ -641,6 +651,8 @@ namespace detail {
result = guess - delta;
if (result <= min)
result = float_next(min);
if (result >= max)
result = float_prior(max);
guess = min;
continue;
}
@@ -669,6 +681,8 @@ namespace detail {
result = guess - delta;
if (result >= max)
result = float_prior(max);
if (result <= min)
result = float_next(min);
guess = min;
continue;
}