File indexing completed on 2025-07-06 08:18:24
0001
0002
0003
0004
0005
0006
0007
0008
0009 #ifndef BOOST_STATS_STUDENTS_T_HPP
0010 #define BOOST_STATS_STUDENTS_T_HPP
0011
0012
0013
0014
0015 #include <boost/math/distributions/fwd.hpp>
0016 #include <boost/math/special_functions/beta.hpp> // for ibeta(a, b, x).
0017 #include <boost/math/special_functions/digamma.hpp>
0018 #include <boost/math/distributions/complement.hpp>
0019 #include <boost/math/distributions/detail/common_error_handling.hpp>
0020 #include <boost/math/distributions/normal.hpp>
0021
0022 #include <utility>
0023
0024 #ifdef _MSC_VER
0025 # pragma warning(push)
0026 # pragma warning(disable: 4702)
0027 #endif
0028
0029 namespace boost { namespace math {
0030
0031 template <class RealType = double, class Policy = policies::policy<> >
0032 class students_t_distribution
0033 {
0034 public:
0035 typedef RealType value_type;
0036 typedef Policy policy_type;
0037
0038 students_t_distribution(RealType df) : df_(df)
0039 {
0040 RealType result;
0041 detail::check_df_gt0_to_inf(
0042 "boost::math::students_t_distribution<%1%>::students_t_distribution", df_, &result, Policy());
0043 }
0044
0045 RealType degrees_of_freedom()const
0046 {
0047 return df_;
0048 }
0049
0050
0051 static RealType find_degrees_of_freedom(
0052 RealType difference_from_mean,
0053 RealType alpha,
0054 RealType beta,
0055 RealType sd,
0056 RealType hint = 100);
0057
0058 private:
0059
0060 RealType df_;
0061 };
0062
0063 typedef students_t_distribution<double> students_t;
0064
0065 #ifdef __cpp_deduction_guides
0066 template <class RealType>
0067 students_t_distribution(RealType)->students_t_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0068 #endif
0069
0070 template <class RealType, class Policy>
0071 inline const std::pair<RealType, RealType> range(const students_t_distribution<RealType, Policy>& )
0072 {
0073
0074 using boost::math::tools::max_value;
0075
0076 return std::pair<RealType, RealType>(((::std::numeric_limits<RealType>::is_specialized & ::std::numeric_limits<RealType>::has_infinity) ? -std::numeric_limits<RealType>::infinity() : -max_value<RealType>()), ((::std::numeric_limits<RealType>::is_specialized & ::std::numeric_limits<RealType>::has_infinity) ? +std::numeric_limits<RealType>::infinity() : +max_value<RealType>()));
0077 }
0078
0079 template <class RealType, class Policy>
0080 inline const std::pair<RealType, RealType> support(const students_t_distribution<RealType, Policy>& )
0081 {
0082
0083
0084 using boost::math::tools::max_value;
0085
0086 return std::pair<RealType, RealType>(((::std::numeric_limits<RealType>::is_specialized & ::std::numeric_limits<RealType>::has_infinity) ? -std::numeric_limits<RealType>::infinity() : -max_value<RealType>()), ((::std::numeric_limits<RealType>::is_specialized & ::std::numeric_limits<RealType>::has_infinity) ? +std::numeric_limits<RealType>::infinity() : +max_value<RealType>()));
0087 }
0088
0089 template <class RealType, class Policy>
0090 inline RealType pdf(const students_t_distribution<RealType, Policy>& dist, const RealType& x)
0091 {
0092 BOOST_FPU_EXCEPTION_GUARD
0093 BOOST_MATH_STD_USING
0094
0095 RealType error_result;
0096 if(false == detail::check_x_not_NaN(
0097 "boost::math::pdf(const students_t_distribution<%1%>&, %1%)", x, &error_result, Policy()))
0098 return error_result;
0099 RealType df = dist.degrees_of_freedom();
0100 if(false == detail::check_df_gt0_to_inf(
0101 "boost::math::pdf(const students_t_distribution<%1%>&, %1%)", df, &error_result, Policy()))
0102 return error_result;
0103
0104 RealType result;
0105 if ((boost::math::isinf)(x))
0106 {
0107 result = static_cast<RealType>(0);
0108 return result;
0109 }
0110 RealType limit = policies::get_epsilon<RealType, Policy>();
0111
0112
0113 limit = static_cast<RealType>(1) / limit;
0114
0115 if (df > limit)
0116 {
0117
0118 normal_distribution<RealType, Policy> n(0, 1);
0119 result = pdf(n, x);
0120 }
0121 else
0122 {
0123 RealType basem1 = x * x / df;
0124 if(basem1 < 0.125)
0125 {
0126 result = exp(-boost::math::log1p(basem1, Policy()) * (1+df) / 2);
0127 }
0128 else
0129 {
0130 result = pow(1 / (1 + basem1), (df + 1) / 2);
0131 }
0132 result /= sqrt(df) * boost::math::beta(df / 2, RealType(0.5f), Policy());
0133 }
0134 return result;
0135 }
0136
0137 template <class RealType, class Policy>
0138 inline RealType cdf(const students_t_distribution<RealType, Policy>& dist, const RealType& x)
0139 {
0140 RealType error_result;
0141
0142 RealType df = dist.degrees_of_freedom();
0143 if (false == detail::check_df_gt0_to_inf(
0144 "boost::math::cdf(const students_t_distribution<%1%>&, %1%)", df, &error_result, Policy()))
0145 {
0146 return error_result;
0147 }
0148
0149 if(false == detail::check_x_not_NaN(
0150 "boost::math::cdf(const students_t_distribution<%1%>&, %1%)", x, &error_result, Policy()))
0151 {
0152 return error_result;
0153 }
0154 if (x == 0)
0155 {
0156 return static_cast<RealType>(0.5);
0157 }
0158 if ((boost::math::isinf)(x))
0159 {
0160 return ((x < 0) ? static_cast<RealType>(0) : static_cast<RealType>(1));
0161 }
0162
0163 RealType limit = policies::get_epsilon<RealType, Policy>();
0164
0165
0166 limit = static_cast<RealType>(1) / limit;
0167
0168 if (df > limit)
0169 {
0170
0171 normal_distribution<RealType, Policy> n(0, 1);
0172 RealType result = cdf(n, x);
0173 return result;
0174 }
0175 else
0176 {
0177
0178
0179
0180
0181
0182
0183
0184
0185
0186
0187
0188
0189
0190
0191
0192
0193
0194
0195 RealType x2 = x * x;
0196 RealType probability;
0197 if(df > 2 * x2)
0198 {
0199 RealType z = x2 / (df + x2);
0200 probability = ibetac(static_cast<RealType>(0.5), df / 2, z, Policy()) / 2;
0201 }
0202 else
0203 {
0204 RealType z = df / (df + x2);
0205 probability = ibeta(df / 2, static_cast<RealType>(0.5), z, Policy()) / 2;
0206 }
0207 return (x > 0 ? 1 - probability : probability);
0208 }
0209 }
0210
0211 template <class RealType, class Policy>
0212 inline RealType quantile(const students_t_distribution<RealType, Policy>& dist, const RealType& p)
0213 {
0214 BOOST_MATH_STD_USING
0215
0216
0217 RealType probability = p;
0218
0219
0220 RealType df = dist.degrees_of_freedom();
0221 static const char* function = "boost::math::quantile(const students_t_distribution<%1%>&, %1%)";
0222 RealType error_result;
0223 if(false == (detail::check_df_gt0_to_inf(
0224 function, df, &error_result, Policy())
0225 && detail::check_probability(function, probability, &error_result, Policy())))
0226 return error_result;
0227
0228 if (probability == 0)
0229 return -policies::raise_overflow_error<RealType>(function, 0, Policy());
0230 if (probability == 1)
0231 return policies::raise_overflow_error<RealType>(function, 0, Policy());
0232 if (probability == static_cast<RealType>(0.5))
0233 return 0;
0234
0235 #if 0
0236
0237
0238
0239
0240 probability = (probability > 0.5) ? 1 - probability : probability;
0241 RealType t, x, y;
0242 x = ibeta_inv(degrees_of_freedom / 2, RealType(0.5), 2 * probability, &y);
0243 if(degrees_of_freedom * y > tools::max_value<RealType>() * x)
0244 t = tools::overflow_error<RealType>(function);
0245 else
0246 t = sqrt(degrees_of_freedom * y / x);
0247
0248
0249
0250 if(p < 0.5)
0251 t = -t;
0252
0253 return t;
0254 #endif
0255
0256
0257
0258
0259
0260
0261
0262 return boost::math::detail::fast_students_t_quantile(df, probability, Policy());
0263 }
0264
0265 template <class RealType, class Policy>
0266 inline RealType cdf(const complemented2_type<students_t_distribution<RealType, Policy>, RealType>& c)
0267 {
0268 return cdf(c.dist, -c.param);
0269 }
0270
0271 template <class RealType, class Policy>
0272 inline RealType quantile(const complemented2_type<students_t_distribution<RealType, Policy>, RealType>& c)
0273 {
0274 return -quantile(c.dist, c.param);
0275 }
0276
0277
0278
0279
0280 namespace detail{
0281
0282
0283
0284 template <class RealType, class Policy>
0285 struct sample_size_func
0286 {
0287 sample_size_func(RealType a, RealType b, RealType s, RealType d)
0288 : alpha(a), beta(b), ratio(s*s/(d*d)) {}
0289
0290 RealType operator()(const RealType& df)
0291 {
0292 if(df <= tools::min_value<RealType>())
0293 {
0294 return 1;
0295 }
0296 students_t_distribution<RealType, Policy> t(df);
0297 RealType qa = quantile(complement(t, alpha));
0298 RealType qb = quantile(complement(t, beta));
0299 qa += qb;
0300 qa *= qa;
0301 qa *= ratio;
0302 qa -= (df + 1);
0303 return qa;
0304 }
0305 RealType alpha, beta, ratio;
0306 };
0307
0308 }
0309
0310 template <class RealType, class Policy>
0311 RealType students_t_distribution<RealType, Policy>::find_degrees_of_freedom(
0312 RealType difference_from_mean,
0313 RealType alpha,
0314 RealType beta,
0315 RealType sd,
0316 RealType hint)
0317 {
0318 static const char* function = "boost::math::students_t_distribution<%1%>::find_degrees_of_freedom";
0319
0320
0321
0322 RealType error_result;
0323 if(false == detail::check_probability(
0324 function, alpha, &error_result, Policy())
0325 && detail::check_probability(function, beta, &error_result, Policy()))
0326 return error_result;
0327
0328 if(hint <= 0)
0329 hint = 1;
0330
0331 detail::sample_size_func<RealType, Policy> f(alpha, beta, sd, difference_from_mean);
0332 tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
0333 std::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
0334 std::pair<RealType, RealType> r = tools::bracket_and_solve_root(f, hint, RealType(2), false, tol, max_iter, Policy());
0335 RealType result = r.first + (r.second - r.first) / 2;
0336 if(max_iter >= policies::get_max_root_iterations<Policy>())
0337 {
0338 return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time: either there is no answer to how many degrees of freedom are required"
0339 " or the answer is infinite. Current best guess is %1%", result, Policy());
0340 }
0341 return result;
0342 }
0343
0344 template <class RealType, class Policy>
0345 inline RealType mode(const students_t_distribution<RealType, Policy>& )
0346 {
0347
0348 return 0;
0349 }
0350
0351 template <class RealType, class Policy>
0352 inline RealType median(const students_t_distribution<RealType, Policy>& )
0353 {
0354
0355 return 0;
0356 }
0357
0358
0359
0360 template <class RealType, class Policy>
0361 inline RealType mean(const students_t_distribution<RealType, Policy>& dist)
0362 {
0363 RealType df = dist.degrees_of_freedom();
0364 if(((boost::math::isnan)(df)) || (df <= 1) )
0365 {
0366 return policies::raise_domain_error<RealType>(
0367 "boost::math::mean(students_t_distribution<%1%> const&, %1%)",
0368 "Mean is undefined for degrees of freedom < 1 but got %1%.", df, Policy());
0369 return std::numeric_limits<RealType>::quiet_NaN();
0370 }
0371 return 0;
0372 }
0373
0374 template <class RealType, class Policy>
0375 inline RealType variance(const students_t_distribution<RealType, Policy>& dist)
0376 {
0377
0378 RealType df = dist.degrees_of_freedom();
0379 if ((boost::math::isnan)(df) || (df <= 2))
0380 {
0381 return policies::raise_domain_error<RealType>(
0382 "boost::math::variance(students_t_distribution<%1%> const&, %1%)",
0383 "variance is undefined for degrees of freedom <= 2, but got %1%.",
0384 df, Policy());
0385 return std::numeric_limits<RealType>::quiet_NaN();
0386 }
0387 if ((boost::math::isinf)(df))
0388 {
0389 return 1;
0390 }
0391 RealType limit = policies::get_epsilon<RealType, Policy>();
0392
0393
0394 limit = static_cast<RealType>(1) / limit;
0395
0396 if (df > limit)
0397 {
0398 return 1;
0399 }
0400 else
0401 {
0402 return df / (df - 2);
0403 }
0404 }
0405
0406 template <class RealType, class Policy>
0407 inline RealType skewness(const students_t_distribution<RealType, Policy>& dist)
0408 {
0409 RealType df = dist.degrees_of_freedom();
0410 if( ((boost::math::isnan)(df)) || (dist.degrees_of_freedom() <= 3))
0411 {
0412 return policies::raise_domain_error<RealType>(
0413 "boost::math::skewness(students_t_distribution<%1%> const&, %1%)",
0414 "Skewness is undefined for degrees of freedom <= 3, but got %1%.",
0415 dist.degrees_of_freedom(), Policy());
0416 return std::numeric_limits<RealType>::quiet_NaN();
0417 }
0418 return 0;
0419 }
0420
0421 template <class RealType, class Policy>
0422 inline RealType kurtosis(const students_t_distribution<RealType, Policy>& dist)
0423 {
0424 RealType df = dist.degrees_of_freedom();
0425 if(((boost::math::isnan)(df)) || (df <= 4))
0426 {
0427 return policies::raise_domain_error<RealType>(
0428 "boost::math::kurtosis(students_t_distribution<%1%> const&, %1%)",
0429 "Kurtosis is undefined for degrees of freedom <= 4, but got %1%.",
0430 df, Policy());
0431 return std::numeric_limits<RealType>::quiet_NaN();
0432 }
0433 if ((boost::math::isinf)(df))
0434 {
0435 return 3;
0436 }
0437 RealType limit = policies::get_epsilon<RealType, Policy>();
0438
0439
0440 limit = static_cast<RealType>(1) / limit;
0441
0442 if (df > limit)
0443 {
0444 return 3;
0445 }
0446 else
0447 {
0448
0449 return 6 / (df - 4) + 3;
0450 }
0451 }
0452
0453 template <class RealType, class Policy>
0454 inline RealType kurtosis_excess(const students_t_distribution<RealType, Policy>& dist)
0455 {
0456
0457
0458 RealType df = dist.degrees_of_freedom();
0459 if(((boost::math::isnan)(df)) || (df <= 4))
0460 {
0461 return policies::raise_domain_error<RealType>(
0462 "boost::math::kurtosis_excess(students_t_distribution<%1%> const&, %1%)",
0463 "Kurtosis_excess is undefined for degrees of freedom <= 4, but got %1%.",
0464 df, Policy());
0465 return std::numeric_limits<RealType>::quiet_NaN();
0466 }
0467 if ((boost::math::isinf)(df))
0468 {
0469 return 0;
0470 }
0471 RealType limit = policies::get_epsilon<RealType, Policy>();
0472
0473
0474 limit = static_cast<RealType>(1) / limit;
0475
0476 if (df > limit)
0477 {
0478 return 0;
0479 }
0480 else
0481 {
0482 return 6 / (df - 4);
0483 }
0484 }
0485
0486 template <class RealType, class Policy>
0487 inline RealType entropy(const students_t_distribution<RealType, Policy>& dist)
0488 {
0489 using std::log;
0490 using std::sqrt;
0491 RealType v = dist.degrees_of_freedom();
0492 RealType vp1 = (v+1)/2;
0493 RealType vd2 = v/2;
0494
0495 return vp1*(digamma(vp1) - digamma(vd2)) + log(sqrt(v)*beta(vd2, RealType(1)/RealType(2)));
0496 }
0497
0498 }
0499 }
0500
0501 #ifdef _MSC_VER
0502 # pragma warning(pop)
0503 #endif
0504
0505
0506
0507
0508 #include <boost/math/distributions/detail/derived_accessors.hpp>
0509
0510 #endif