Skip to content

Commit a10f3c8

Browse files
jzmaddockjzmaddock
authored
Merge pull request #1355 from boostorg/nc_t_find
Add parameter finders to the non-central-T.
2 parents 0b5e38a + 585d317 commit a10f3c8

File tree

2 files changed

+68
-51
lines changed

2 files changed

+68
-51
lines changed

‎include/boost/math/distributions/non_central_t.hpp‎

Lines changed: 36 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -719,11 +719,6 @@ namespace boost
719719
return result;
720720
}
721721

722-
#if 0
723-
//
724-
// This code is disabled, since there can be multiple answers to the
725-
// question, and it's not clear how to find the "right" one.
726-
//
727722
template <class RealType, class Policy>
728723
struct t_degrees_of_freedom_finder
729724
{
@@ -749,13 +744,19 @@ namespace boost
749744
inline RealType find_t_degrees_of_freedom(
750745
RealType delta, RealType x, RealType p, RealType q, const Policy& pol)
751746
{
747+
using std::fabs;
752748
const char* function = "non_central_t<%1%>::find_degrees_of_freedom";
753749
if((p == 0) || (q == 0))
754750
{
755751
//
756-
// Can't a thing if one of p and q is zero:
752+
// Can't find a thing if one of p and q is zero:
757753
//
758-
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
754+
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%",
755+
RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy()); // LCOV_EXCL_LINE
756+
}
757+
if (fabs(x) < tools::epsilon<RealType>())
758+
{
759+
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!!",
759760
RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy()); // LCOV_EXCL_LINE
760761
}
761762
t_degrees_of_freedom_finder<RealType, Policy> f(delta, x, p < q ? p : q, p < q ? false : true);
@@ -766,7 +767,7 @@ namespace boost
766767
//
767768
RealType guess = 200;
768769
std::pair<RealType, RealType> ir = tools::bracket_and_solve_root(
769-
f, guess, RealType(2), false, tol, max_iter, pol);
770+
f, guess, RealType(2), x < 0 ? false : true, tol, max_iter, pol);
770771
RealType result = ir.first + (ir.second - ir.first) / 2;
771772
if(max_iter >= policies::get_max_root_iterations<Policy>())
772773
{
@@ -819,9 +820,9 @@ namespace boost
819820
//
820821
RealType guess;
821822
if(f(0) < 0)
822-
guess = 1;
823-
else
824823
guess = -1;
824+
else
825+
guess = 1;
825826
std::pair<RealType, RealType> ir = tools::bracket_and_solve_root(
826827
f, guess, RealType(2), false, tol, max_iter, pol);
827828
RealType result = ir.first + (ir.second - ir.first) / 2;
@@ -832,7 +833,6 @@ namespace boost
832833
}
833834
return result;
834835
}
835-
#endif
836836
} // namespace detail ======================================================================
837837

838838
template <class RealType = double, class Policy = policies::policy<> >
@@ -864,26 +864,22 @@ namespace boost
864864
{ // Private data getter function.
865865
return ncp;
866866
}
867-
#if 0
868-
//
869-
// This code is disabled, since there can be multiple answers to the
870-
// question, and it's not clear how to find the "right" one.
871-
//
867+
872868
static RealType find_degrees_of_freedom(RealType delta, RealType x, RealType p)
873869
{
874870
const char* function = "non_central_t<%1%>::find_degrees_of_freedom";
875-
typedef typename policies::evaluation<RealType, Policy>::type value_type;
871+
typedef typename policies::evaluation<RealType, Policy>::type eval_type;
876872
typedef typename policies::normalise<
877873
Policy,
878874
policies::promote_float<false>,
879875
policies::promote_double<false>,
880876
policies::discrete_quantile<>,
881877
policies::assert_undefined<> >::type forwarding_policy;
882-
value_type result = detail::find_t_degrees_of_freedom(
883-
static_cast<value_type>(delta),
884-
static_cast<value_type>(x),
885-
static_cast<value_type>(p),
886-
static_cast<value_type>(1-p),
878+
eval_type result = detail::find_t_degrees_of_freedom(
879+
static_cast<eval_type>(delta),
880+
static_cast<eval_type>(x),
881+
static_cast<eval_type>(p),
882+
static_cast<eval_type>(1-p),
887883
forwarding_policy());
888884
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
889885
result,
@@ -893,18 +889,18 @@ namespace boost
893889
static RealType find_degrees_of_freedom(const complemented3_type<A,B,C>& c)
894890
{
895891
const char* function = "non_central_t<%1%>::find_degrees_of_freedom";
896-
typedef typename policies::evaluation<RealType, Policy>::type value_type;
892+
typedef typename policies::evaluation<RealType, Policy>::type eval_type;
897893
typedef typename policies::normalise<
898894
Policy,
899895
policies::promote_float<false>,
900896
policies::promote_double<false>,
901897
policies::discrete_quantile<>,
902898
policies::assert_undefined<> >::type forwarding_policy;
903-
value_type result = detail::find_t_degrees_of_freedom(
904-
static_cast<value_type>(c.dist),
905-
static_cast<value_type>(c.param1),
906-
static_cast<value_type>(1-c.param2),
907-
static_cast<value_type>(c.param2),
899+
eval_type result = detail::find_t_degrees_of_freedom(
900+
static_cast<eval_type>(c.dist),
901+
static_cast<eval_type>(c.param1),
902+
static_cast<eval_type>(1-c.param2),
903+
static_cast<eval_type>(c.param2),
908904
forwarding_policy());
909905
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
910906
result,
@@ -913,18 +909,18 @@ namespace boost
913909
static RealType find_non_centrality(RealType v, RealType x, RealType p)
914910
{
915911
const char* function = "non_central_t<%1%>::find_t_non_centrality";
916-
typedef typename policies::evaluation<RealType, Policy>::type value_type;
912+
typedef typename policies::evaluation<RealType, Policy>::type eval_type;
917913
typedef typename policies::normalise<
918914
Policy,
919915
policies::promote_float<false>,
920916
policies::promote_double<false>,
921917
policies::discrete_quantile<>,
922918
policies::assert_undefined<> >::type forwarding_policy;
923-
value_type result = detail::find_t_non_centrality(
924-
static_cast<value_type>(v),
925-
static_cast<value_type>(x),
926-
static_cast<value_type>(p),
927-
static_cast<value_type>(1-p),
919+
eval_type result = detail::find_t_non_centrality(
920+
static_cast<eval_type>(v),
921+
static_cast<eval_type>(x),
922+
static_cast<eval_type>(p),
923+
static_cast<eval_type>(1-p),
928924
forwarding_policy());
929925
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
930926
result,
@@ -934,24 +930,23 @@ namespace boost
934930
static RealType find_non_centrality(const complemented3_type<A,B,C>& c)
935931
{
936932
const char* function = "non_central_t<%1%>::find_t_non_centrality";
937-
typedef typename policies::evaluation<RealType, Policy>::type value_type;
933+
typedef typename policies::evaluation<RealType, Policy>::type eval_type;
938934
typedef typename policies::normalise<
939935
Policy,
940936
policies::promote_float<false>,
941937
policies::promote_double<false>,
942938
policies::discrete_quantile<>,
943939
policies::assert_undefined<> >::type forwarding_policy;
944-
value_type result = detail::find_t_non_centrality(
945-
static_cast<value_type>(c.dist),
946-
static_cast<value_type>(c.param1),
947-
static_cast<value_type>(1-c.param2),
948-
static_cast<value_type>(c.param2),
940+
eval_type result = detail::find_t_non_centrality(
941+
static_cast<eval_type>(c.dist),
942+
static_cast<eval_type>(c.param1),
943+
static_cast<eval_type>(1-c.param2),
944+
static_cast<eval_type>(c.param2),
949945
forwarding_policy());
950946
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
951947
result,
952948
function);
953949
}
954-
#endif
955950
private:
956951
// Data member, initialized by constructor.
957952
RealType v; // degrees of freedom

‎test/test_nc_t.hpp‎

Lines changed: 32 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -494,38 +494,60 @@ void quantile_sanity_check(T& data, const char* type_name, const char* test)
494494
}
495495
catch(const boost::math::evaluation_error&) {}
496496
#endif
497-
#if 0
498497
//
499498
// Sanity check degrees-of-freedom finder, don't bother at float
500499
// precision though as there's not enough data in the probability
501500
// values to get back to the correct degrees of freedom or
502501
// non-centrality parameter:
503502
//
504503
try{
504+
Real non_centrality_precision_multiplier = 1;
505+
Real df_precision_multiplier = 1;
506+
if (data[i][0] > 10000)
507+
{
508+
non_centrality_precision_multiplier = 500;
509+
df_precision_multiplier = 500;
510+
}
511+
if (data[i][0] > 1000000000)
512+
{
513+
df_precision_multiplier *= 50; // very little precision left at this point
514+
}
505515
if((data[i][3] < 0.99) && (data[i][3] != 0))
506516
{
507-
BOOST_CHECK_CLOSE_EX(
508-
boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(data[i][1], data[i][2], data[i][3]),
509-
data[i][0], precision, i);
517+
if (data[i][0] < 1.0e12) // no precision above this
518+
{
519+
BOOST_CHECK_CLOSE_EX(
520+
boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(data[i][1], data[i][2], data[i][3]),
521+
data[i][0], precision * df_precision_multiplier, i);
522+
}
510523
BOOST_CHECK_CLOSE_EX(
511524
boost::math::non_central_t_distribution<value_type>::find_non_centrality(data[i][0], data[i][2], data[i][3]),
512-
data[i][1], precision, i);
525+
data[i][1], precision * non_centrality_precision_multiplier, i);
513526
}
514527
if((data[i][4] < 0.99) && (data[i][4] != 0))
515528
{
516-
BOOST_CHECK_CLOSE_EX(
517-
boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(boost::math::complement(data[i][1], data[i][2], data[i][4])),
518-
data[i][0], precision, i);
529+
if (data[i][0] < 1.0e12) // no precision above this
530+
{
531+
BOOST_CHECK_CLOSE_EX(
532+
boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(boost::math::complement(data[i][1], data[i][2], data[i][4])),
533+
data[i][0], precision * df_precision_multiplier, i);
534+
}
519535
BOOST_CHECK_CLOSE_EX(
520536
boost::math::non_central_t_distribution<value_type>::find_non_centrality(boost::math::complement(data[i][0], data[i][2], data[i][4])),
521-
data[i][1], precision, i);
537+
data[i][1], precision * non_centrality_precision_multiplier, i);
522538
}
523539
}
524540
catch(const std::exception& e)
525541
{
526542
BOOST_ERROR(e.what());
527543
}
528-
#endif
544+
// Code coverage:
545+
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);
546+
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);
547+
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);
548+
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);
549+
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);
550+
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);
529551
}
530552
}
531553
#endif

0 commit comments

Comments
 (0)