boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-02-24 04:02:18 +00:00
Code Issues Projects Releases Wiki Activity

Merge pull request #1355 from boostorg/nc_t_find

Browse Source 
Add parameter finders to the non-central-T.
This commit is contained in:
jzmaddock
2026-02-19 17:35:55 +00:00
committed by GitHub
parent 0b5e38aa25 585d317947
commit a10f3c86fc
2 changed files with 68 additions and 51 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -719,11 +719,6 @@ namespace boost
return result;
}
#if 0
//
// This code is disabled, since there can be multiple answers to the
// question, and it's not clear how to find the "right" one.
//
template <class RealType, class Policy>
struct t_degrees_of_freedom_finder
{
@@ -749,13 +744,19 @@ namespace boost
inline RealType find_t_degrees_of_freedom(
RealType delta, RealType x, RealType p, RealType q, const Policy& pol)
{
using std::fabs;
const char* function = "non_central_t<%1%>::find_degrees_of_freedom";
if((p == 0) || (q == 0))
{
//
// Can't a thing if one of p and q is zero:
// Can't find a thing if one of p and q is zero:
//
return policies::raise_evaluation_error<RealType>(function, "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%", // LCOV_EXCL_LINE
return policies::raise_evaluation_error<RealType>(function, "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%",
RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy()); // LCOV_EXCL_LINE
}
if (fabs(x) < tools::epsilon<RealType>())
{
return policies::raise_evaluation_error<RealType>(function, "Can't find degrees of freedom when the abscissa value is very close to zero as all degrees of freedom generate the same CDF at x=0: try again further out in the tails!!",
RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy()); // LCOV_EXCL_LINE
}
t_degrees_of_freedom_finder<RealType, Policy> f(delta, x, p < q ? p : q, p < q ? false : true);
@@ -766,7 +767,7 @@ namespace boost
//
RealType guess = 200;
std::pair<RealType, RealType> ir = tools::bracket_and_solve_root(
f, guess, RealType(2), false, tol, max_iter, pol);
f, guess, RealType(2), x < 0 ? false : true, tol, max_iter, pol);
RealType result = ir.first + (ir.second - ir.first) / 2;
if(max_iter >= policies::get_max_root_iterations<Policy>())
{
@@ -819,9 +820,9 @@ namespace boost
//
RealType guess;
if(f(0) < 0)
guess = 1;
else
guess = -1;
else
guess = 1;
std::pair<RealType, RealType> ir = tools::bracket_and_solve_root(
f, guess, RealType(2), false, tol, max_iter, pol);
RealType result = ir.first + (ir.second - ir.first) / 2;
@@ -832,7 +833,6 @@ namespace boost
}
return result;
}
#endif
} // namespace detail ======================================================================
template <class RealType = double, class Policy = policies::policy<> >
@@ -864,26 +864,22 @@ namespace boost
{ // Private data getter function.
return ncp;
}
#if 0
//
// This code is disabled, since there can be multiple answers to the
// question, and it's not clear how to find the "right" one.
//
static RealType find_degrees_of_freedom(RealType delta, RealType x, RealType p)
{
const char* function = "non_central_t<%1%>::find_degrees_of_freedom";
typedef typename policies::evaluation<RealType, Policy>::type value_type;
typedef typename policies::evaluation<RealType, Policy>::type eval_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
value_type result = detail::find_t_degrees_of_freedom(
static_cast<value_type>(delta),
static_cast<value_type>(x),
static_cast<value_type>(p),
static_cast<value_type>(1-p),
eval_type result = detail::find_t_degrees_of_freedom(
static_cast<eval_type>(delta),
static_cast<eval_type>(x),
static_cast<eval_type>(p),
static_cast<eval_type>(1-p),
forwarding_policy());
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
result,
@@ -893,18 +889,18 @@ namespace boost
static RealType find_degrees_of_freedom(const complemented3_type<A,B,C>& c)
{
const char* function = "non_central_t<%1%>::find_degrees_of_freedom";
typedef typename policies::evaluation<RealType, Policy>::type value_type;
typedef typename policies::evaluation<RealType, Policy>::type eval_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
value_type result = detail::find_t_degrees_of_freedom(
static_cast<value_type>(c.dist),
static_cast<value_type>(c.param1),
static_cast<value_type>(1-c.param2),
static_cast<value_type>(c.param2),
eval_type result = detail::find_t_degrees_of_freedom(
static_cast<eval_type>(c.dist),
static_cast<eval_type>(c.param1),
static_cast<eval_type>(1-c.param2),
static_cast<eval_type>(c.param2),
forwarding_policy());
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
result,
@@ -913,18 +909,18 @@ namespace boost
static RealType find_non_centrality(RealType v, RealType x, RealType p)
{
const char* function = "non_central_t<%1%>::find_t_non_centrality";
typedef typename policies::evaluation<RealType, Policy>::type value_type;
typedef typename policies::evaluation<RealType, Policy>::type eval_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
value_type result = detail::find_t_non_centrality(
static_cast<value_type>(v),
static_cast<value_type>(x),
static_cast<value_type>(p),
static_cast<value_type>(1-p),
eval_type result = detail::find_t_non_centrality(
static_cast<eval_type>(v),
static_cast<eval_type>(x),
static_cast<eval_type>(p),
static_cast<eval_type>(1-p),
forwarding_policy());
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
result,
@@ -934,24 +930,23 @@ namespace boost
static RealType find_non_centrality(const complemented3_type<A,B,C>& c)
{
const char* function = "non_central_t<%1%>::find_t_non_centrality";
typedef typename policies::evaluation<RealType, Policy>::type value_type;
typedef typename policies::evaluation<RealType, Policy>::type eval_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
value_type result = detail::find_t_non_centrality(
static_cast<value_type>(c.dist),
static_cast<value_type>(c.param1),
static_cast<value_type>(1-c.param2),
static_cast<value_type>(c.param2),
eval_type result = detail::find_t_non_centrality(
static_cast<eval_type>(c.dist),
static_cast<eval_type>(c.param1),
static_cast<eval_type>(1-c.param2),
static_cast<eval_type>(c.param2),
forwarding_policy());
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
result,
function);
}
#endif
private:
// Data member, initialized by constructor.
RealType v; // degrees of freedom

42
test/test_nc_t.hpp
View File

@@ -494,7 +494,6 @@ void quantile_sanity_check(T& data, const char* type_name, const char* test)
}
catch(const boost::math::evaluation_error&) {}
#endif
#if 0
//
// Sanity check degrees-of-freedom finder, don't bother at float
// precision though as there's not enough data in the probability
@@ -502,30 +501,53 @@ void quantile_sanity_check(T& data, const char* type_name, const char* test)
// non-centrality parameter:
//
try{
Real non_centrality_precision_multiplier = 1;
Real df_precision_multiplier = 1;
if (data[i][0] > 10000)
{
non_centrality_precision_multiplier = 500;
df_precision_multiplier = 500;
}
if (data[i][0] > 1000000000)
{
df_precision_multiplier *= 50; // very little precision left at this point
}
if((data[i][3] < 0.99) && (data[i][3] != 0))
{
BOOST_CHECK_CLOSE_EX(
boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(data[i][1], data[i][2], data[i][3]),
data[i][0], precision, i);
if (data[i][0] < 1.0e12) // no precision above this
{
BOOST_CHECK_CLOSE_EX(
boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(data[i][1], data[i][2], data[i][3]),
data[i][0], precision * df_precision_multiplier, i);
}
BOOST_CHECK_CLOSE_EX(
boost::math::non_central_t_distribution<value_type>::find_non_centrality(data[i][0], data[i][2], data[i][3]),
data[i][1], precision, i);
data[i][1], precision * non_centrality_precision_multiplier, i);
}
if((data[i][4] < 0.99) && (data[i][4] != 0))
{
BOOST_CHECK_CLOSE_EX(
boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(boost::math::complement(data[i][1], data[i][2], data[i][4])),
data[i][0], precision, i);
if (data[i][0] < 1.0e12) // no precision above this
{
BOOST_CHECK_CLOSE_EX(
boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(boost::math::complement(data[i][1], data[i][2], data[i][4])),
data[i][0], precision * df_precision_multiplier, i);
}
BOOST_CHECK_CLOSE_EX(
boost::math::non_central_t_distribution<value_type>::find_non_centrality(boost::math::complement(data[i][0], data[i][2], data[i][4])),
data[i][1], precision, i);
data[i][1], precision * non_centrality_precision_multiplier, i);
}
}
catch(const std::exception& e)
{
BOOST_ERROR(e.what());
}
#endif
// Code coverage:
BOOST_CHECK_THROW(boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(data[i][1], boost::math::tools::epsilon<value_type>() / 2, data[i][3]), boost::math::evaluation_error);
BOOST_CHECK_THROW(boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(boost::math::complement(data[i][1], boost::math::tools::epsilon<value_type>() / 2, data[i][3])), boost::math::evaluation_error);
BOOST_CHECK_THROW(boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(data[i][1], data[i][2], value_type(0)), boost::math::evaluation_error);
BOOST_CHECK_THROW(boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(boost::math::complement(data[i][1], data[i][2], value_type(0))), boost::math::evaluation_error);
BOOST_CHECK_THROW(boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(data[i][1], data[i][2], value_type(1)), boost::math::evaluation_error);
BOOST_CHECK_THROW(boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(boost::math::complement(data[i][1], data[i][2], value_type(1))), boost::math::evaluation_error);
}
}
#endif