Skip to content

Add parameter finders to the non-central-T. #1355

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
jzmaddock wants to merge 1 commit into
base: develop
Choose a base branch
from
+60 −50
Open

Add parameter finders to the non-central-T. #1355

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
75 changes: 35 additions & 40 deletions include/boost/math/distributions/non_central_t.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
{
Expand All @@ -749,15 +744,21 @@ 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
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);
tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
Expand All @@ -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>())
{
Expand Down Expand Up @@ -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;
Expand All @@ -832,7 +833,6 @@ namespace boost
}
return result;
}
#endif
} // namespace detail ======================================================================

template <class RealType = double, class Policy = policies::policy<> >
Expand Down Expand Up @@ -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,
Expand All @@ -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,
Expand All @@ -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,
Expand All @@ -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
Expand Down
35 changes: 25 additions & 10 deletions test/test_nc_t.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -494,38 +494,53 @@ 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
// values to get back to the correct degrees of freedom or
// 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
}
}
#endif
Expand Down
Loading