boost/math/distributions/non_central_t.hpp
// boost\math\distributions\non_central_t.hpp
// Copyright John Maddock 2008.
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt
// or copy at http://www.boost.org/LICENSE_1_0.txt)
#ifndef BOOST_MATH_SPECIAL_NON_CENTRAL_T_HPP
#define BOOST_MATH_SPECIAL_NON_CENTRAL_T_HPP
#include <boost/math/distributions/fwd.hpp>
#include <boost/math/distributions/non_central_beta.hpp> // for nc beta
#include <boost/math/distributions/normal.hpp> // for normal CDF and quantile
#include <boost/math/distributions/students_t.hpp>
#include <boost/math/distributions/detail/generic_quantile.hpp> // quantile
namespace boost
{
namespace math
{
template <class RealType, class Policy>
class non_central_t_distribution;
namespace detail{
template <class T, class Policy>
T non_central_t2_p(T v, T delta, T x, T y, const Policy& pol, T init_val)
{
BOOST_MATH_STD_USING
//
// Variables come first:
//
std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
T errtol = policies::get_epsilon<T, Policy>();
T d2 = delta * delta / 2;
//
// k is the starting point for iteration, and is the
// maximum of the poisson weighting term, we don't
// ever allow k == 0 as this can lead to catastrophic
// cancellation errors later (test case is v = 1621286869049072.3
// delta = 0.16212868690490723, x = 0.86987415482475994).
//
int k = itrunc(d2);
T pois;
if(k == 0) k = 1;
// Starting Poisson weight:
pois = gamma_p_derivative(T(k+1), d2, pol)
* tgamma_delta_ratio(T(k + 1), T(0.5f))
* delta / constants::root_two<T>();
if(pois == 0)
return init_val;
T xterm, beta;
// Recurrence & starting beta terms:
beta = x < y
? detail::ibeta_imp(T(k + 1), T(v / 2), x, pol, false, true, &xterm)
: detail::ibeta_imp(T(v / 2), T(k + 1), y, pol, true, true, &xterm);
xterm *= y / (v / 2 + k);
T poisf(pois), betaf(beta), xtermf(xterm);
T sum = init_val;
if((xterm == 0) && (beta == 0))
return init_val;
//
// Backwards recursion first, this is the stable
// direction for recursion:
//
std::uintmax_t count = 0;
T last_term = 0;
for(int i = k; i >= 0; --i)
{
T term = beta * pois;
sum += term;
// Don't terminate on first term in case we "fixed" k above:
if((fabs(last_term) > fabs(term)) && fabs(term/sum) < errtol)
break;
last_term = term;
pois *= (i + 0.5f) / d2;
beta += xterm;
xterm *= (i) / (x * (v / 2 + i - 1));
++count;
}
last_term = 0;
for(int i = k + 1; ; ++i)
{
poisf *= d2 / (i + 0.5f);
xtermf *= (x * (v / 2 + i - 1)) / (i);
betaf -= xtermf;
T term = poisf * betaf;
sum += term;
if((fabs(last_term) >= fabs(term)) && (fabs(term/sum) < errtol))
break;
last_term = term;
++count;
if(count > max_iter)
{
return policies::raise_evaluation_error(
"cdf(non_central_t_distribution<%1%>, %1%)",
"Series did not converge, closest value was %1%", sum, pol);
}
}
return sum;
}
template <class T, class Policy>
T non_central_t2_q(T v, T delta, T x, T y, const Policy& pol, T init_val)
{
BOOST_MATH_STD_USING
//
// Variables come first:
//
std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
T errtol = boost::math::policies::get_epsilon<T, Policy>();
T d2 = delta * delta / 2;
//
// k is the starting point for iteration, and is the
// maximum of the poisson weighting term, we don't allow
// k == 0 as this can cause catastrophic cancellation errors
// (test case is v = 561908036470413.25, delta = 0.056190803647041321,
// x = 1.6155232703966216):
//
int k = itrunc(d2);
if(k == 0) k = 1;
// Starting Poisson weight:
T pois;
if((k < (int)(max_factorial<T>::value)) && (d2 < tools::log_max_value<T>()) && (log(d2) * k < tools::log_max_value<T>()))
{
//
// For small k we can optimise this calculation by using
// a simpler reduced formula:
//
pois = exp(-d2);
pois *= pow(d2, static_cast<T>(k));
pois /= boost::math::tgamma(T(k + 1 + 0.5), pol);
pois *= delta / constants::root_two<T>();
}
else
{
pois = gamma_p_derivative(T(k+1), d2, pol)
* tgamma_delta_ratio(T(k + 1), T(0.5f))
* delta / constants::root_two<T>();
}
if(pois == 0)
return init_val;
// Recurance term:
T xterm;
T beta;
// Starting beta term:
if(k != 0)
{
beta = x < y
? detail::ibeta_imp(T(k + 1), T(v / 2), x, pol, true, true, &xterm)
: detail::ibeta_imp(T(v / 2), T(k + 1), y, pol, false, true, &xterm);
xterm *= y / (v / 2 + k);
}
else
{
beta = pow(y, v / 2);
xterm = beta;
}
T poisf(pois), betaf(beta), xtermf(xterm);
T sum = init_val;
if((xterm == 0) && (beta == 0))
return init_val;
//
// Fused forward and backwards recursion:
//
std::uintmax_t count = 0;
T last_term = 0;
for(int i = k + 1, j = k; ; ++i, --j)
{
poisf *= d2 / (i + 0.5f);
xtermf *= (x * (v / 2 + i - 1)) / (i);
betaf += xtermf;
T term = poisf * betaf;
if(j >= 0)
{
term += beta * pois;
pois *= (j + 0.5f) / d2;
beta -= xterm;
xterm *= (j) / (x * (v / 2 + j - 1));
}
sum += term;
// Don't terminate on first term in case we "fixed" the value of k above:
if((fabs(last_term) > fabs(term)) && fabs(term/sum) < errtol)
break;
last_term = term;
if(count > max_iter)
{
return policies::raise_evaluation_error(
"cdf(non_central_t_distribution<%1%>, %1%)",
"Series did not converge, closest value was %1%", sum, pol);
}
++count;
}
return sum;
}
template <class T, class Policy>
T non_central_t_cdf(T v, T delta, T t, bool invert, const Policy& pol)
{
BOOST_MATH_STD_USING
if ((boost::math::isinf)(v))
{ // Infinite degrees of freedom, so use normal distribution located at delta.
normal_distribution<T, Policy> n(delta, 1);
return cdf(n, t);
}
//
// Otherwise, for t < 0 we have to use the reflection formula:
if(t < 0)
{
t = -t;
delta = -delta;
invert = !invert;
}
if(fabs(delta / (4 * v)) < policies::get_epsilon<T, Policy>())
{
// Approximate with a Student's T centred on delta,
// the crossover point is based on eq 2.6 from
// "A Comparison of Approximations To Percentiles of the
// Noncentral t-Distribution". H. Sahai and M. M. Ojeda,
// Revista Investigacion Operacional Vol 21, No 2, 2000.
// Original sources referenced in the above are:
// "Some Approximations to the Percentage Points of the Noncentral
// t-Distribution". C. van Eeden. International Statistical Review, 29, 4-31.
// "Continuous Univariate Distributions". N.L. Johnson, S. Kotz and
// N. Balkrishnan. 1995. John Wiley and Sons New York.
T result = cdf(students_t_distribution<T, Policy>(v), t - delta);
return invert ? 1 - result : result;
}
//
// x and y are the corresponding random
// variables for the noncentral beta distribution,
// with y = 1 - x:
//
T x = t * t / (v + t * t);
T y = v / (v + t * t);
T d2 = delta * delta;
T a = 0.5f;
T b = v / 2;
T c = a + b + d2 / 2;
//
// Crossover point for calculating p or q is the same
// as for the noncentral beta:
//
T cross = 1 - (b / c) * (1 + d2 / (2 * c * c));
T result;
if(x < cross)
{
//
// Calculate p:
//
if(x != 0)
{
result = non_central_beta_p(a, b, d2, x, y, pol);
result = non_central_t2_p(v, delta, x, y, pol, result);
result /= 2;
}
else
result = 0;
result += cdf(boost::math::normal_distribution<T, Policy>(), -delta);
}
else
{
//
// Calculate q:
//
invert = !invert;
if(x != 0)
{
result = non_central_beta_q(a, b, d2, x, y, pol);
result = non_central_t2_q(v, delta, x, y, pol, result);
result /= 2;
}
else // x == 0
result = cdf(complement(boost::math::normal_distribution<T, Policy>(), -delta));
}
if(invert)
result = 1 - result;
return result;
}
template <class T, class Policy>
T non_central_t_quantile(const char* function, T v, T delta, T p, T q, const Policy&)
{
BOOST_MATH_STD_USING
// static const char* function = "quantile(non_central_t_distribution<%1%>, %1%)";
// now passed as function
typedef typename policies::evaluation<T, Policy>::type value_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
T r;
if(!detail::check_df_gt0_to_inf(
function,
v, &r, Policy())
||
!detail::check_finite(
function,
delta,
&r,
Policy())
||
!detail::check_probability(
function,
p,
&r,
Policy()))
return r;
value_type guess = 0;
if ( ((boost::math::isinf)(v)) || (v > 1 / boost::math::tools::epsilon<T>()) )
{ // Infinite or very large degrees of freedom, so use normal distribution located at delta.
normal_distribution<T, Policy> n(delta, 1);
if (p < q)
{
return quantile(n, p);
}
else
{
return quantile(complement(n, q));
}
}
else if(v > 3)
{ // Use normal distribution to calculate guess.
value_type mean = (v > 1 / policies::get_epsilon<T, Policy>()) ? delta : delta * sqrt(v / 2) * tgamma_delta_ratio((v - 1) * 0.5f, T(0.5f));
value_type var = (v > 1 / policies::get_epsilon<T, Policy>()) ? value_type(1) : (((delta * delta + 1) * v) / (v - 2) - mean * mean);
if(p < q)
guess = quantile(normal_distribution<value_type, forwarding_policy>(mean, var), p);
else
guess = quantile(complement(normal_distribution<value_type, forwarding_policy>(mean, var), q));
}
//
// We *must* get the sign of the initial guess correct,
// or our root-finder will fail, so double check it now:
//
value_type pzero = non_central_t_cdf(
static_cast<value_type>(v),
static_cast<value_type>(delta),
static_cast<value_type>(0),
!(p < q),
forwarding_policy());
int s;
if(p < q)
s = boost::math::sign(p - pzero);
else
s = boost::math::sign(pzero - q);
if(s != boost::math::sign(guess))
{
guess = static_cast<T>(s);
}
value_type result = detail::generic_quantile(
non_central_t_distribution<value_type, forwarding_policy>(v, delta),
(p < q ? p : q),
guess,
(p >= q),
function);
return policies::checked_narrowing_cast<T, forwarding_policy>(
result,
function);
}
template <class T, class Policy>
T non_central_t2_pdf(T n, T delta, T x, T y, const Policy& pol, T init_val)
{
BOOST_MATH_STD_USING
//
// Variables come first:
//
std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
T errtol = boost::math::policies::get_epsilon<T, Policy>();
T d2 = delta * delta / 2;
//
// k is the starting point for iteration, and is the
// maximum of the poisson weighting term:
//
int k = itrunc(d2);
T pois, xterm;
if(k == 0)
k = 1;
// Starting Poisson weight:
pois = gamma_p_derivative(T(k+1), d2, pol)
* tgamma_delta_ratio(T(k + 1), T(0.5f))
* delta / constants::root_two<T>();
// Starting beta term:
xterm = x < y
? ibeta_derivative(T(k + 1), n / 2, x, pol)
: ibeta_derivative(n / 2, T(k + 1), y, pol);
T poisf(pois), xtermf(xterm);
T sum = init_val;
if((pois == 0) || (xterm == 0))
return init_val;
//
// Backwards recursion first, this is the stable
// direction for recursion:
//
std::uintmax_t count = 0;
for(int i = k; i >= 0; --i)
{
T term = xterm * pois;
sum += term;
if(((fabs(term/sum) < errtol) && (i != k)) || (term == 0))
break;
pois *= (i + 0.5f) / d2;
xterm *= (i) / (x * (n / 2 + i));
++count;
if(count > max_iter)
{
return policies::raise_evaluation_error(
"pdf(non_central_t_distribution<%1%>, %1%)",
"Series did not converge, closest value was %1%", sum, pol);
}
}
for(int i = k + 1; ; ++i)
{
poisf *= d2 / (i + 0.5f);
xtermf *= (x * (n / 2 + i)) / (i);
T term = poisf * xtermf;
sum += term;
if((fabs(term/sum) < errtol) || (term == 0))
break;
++count;
if(count > max_iter)
{
return policies::raise_evaluation_error(
"pdf(non_central_t_distribution<%1%>, %1%)",
"Series did not converge, closest value was %1%", sum, pol);
}
}
return sum;
}
template <class T, class Policy>
T non_central_t_pdf(T n, T delta, T t, const Policy& pol)
{
BOOST_MATH_STD_USING
if ((boost::math::isinf)(n))
{ // Infinite degrees of freedom, so use normal distribution located at delta.
normal_distribution<T, Policy> norm(delta, 1);
return pdf(norm, t);
}
//
// Otherwise, for t < 0 we have to use the reflection formula:
if(t < 0)
{
t = -t;
delta = -delta;
}
if(t == 0)
{
//
// Handle this as a special case, using the formula
// from Weisstein, Eric W.
// "Noncentral Student's t-Distribution."
// From MathWorld--A Wolfram Web Resource.
// http://mathworld.wolfram.com/NoncentralStudentst-Distribution.html
//
// The formula is simplified thanks to the relation
// 1F1(a,b,0) = 1.
//
return tgamma_delta_ratio(n / 2 + 0.5f, T(0.5f))
* sqrt(n / constants::pi<T>())
* exp(-delta * delta / 2) / 2;
}
if(fabs(delta / (4 * n)) < policies::get_epsilon<T, Policy>())
{
// Approximate with a Student's T centred on delta,
// the crossover point is based on eq 2.6 from
// "A Comparison of Approximations To Percentiles of the
// Noncentral t-Distribution". H. Sahai and M. M. Ojeda,
// Revista Investigacion Operacional Vol 21, No 2, 2000.
// Original sources referenced in the above are:
// "Some Approximations to the Percentage Points of the Noncentral
// t-Distribution". C. van Eeden. International Statistical Review, 29, 4-31.
// "Continuous Univariate Distributions". N.L. Johnson, S. Kotz and
// N. Balkrishnan. 1995. John Wiley and Sons New York.
return pdf(students_t_distribution<T, Policy>(n), t - delta);
}
//
// x and y are the corresponding random
// variables for the noncentral beta distribution,
// with y = 1 - x:
//
T x = t * t / (n + t * t);
T y = n / (n + t * t);
T a = 0.5f;
T b = n / 2;
T d2 = delta * delta;
//
// Calculate pdf:
//
T dt = n * t / (n * n + 2 * n * t * t + t * t * t * t);
T result = non_central_beta_pdf(a, b, d2, x, y, pol);
T tol = tools::epsilon<T>() * result * 500;
result = non_central_t2_pdf(n, delta, x, y, pol, result);
if(result <= tol)
result = 0;
result *= dt;
return result;
}
template <class T, class Policy>
T mean(T v, T delta, const Policy& pol)
{
if ((boost::math::isinf)(v))
{
return delta;
}
BOOST_MATH_STD_USING
if (v > 1 / boost::math::tools::epsilon<T>() )
{
//normal_distribution<T, Policy> n(delta, 1);
//return boost::math::mean(n);
return delta;
}
else
{
return delta * sqrt(v / 2) * tgamma_delta_ratio((v - 1) * 0.5f, T(0.5f), pol);
}
// Other moments use mean so using normal distribution is propagated.
}
template <class T, class Policy>
T variance(T v, T delta, const Policy& pol)
{
if ((boost::math::isinf)(v))
{
return 1;
}
if (delta == 0)
{ // == Student's t
return v / (v - 2);
}
T result = ((delta * delta + 1) * v) / (v - 2);
T m = mean(v, delta, pol);
result -= m * m;
return result;
}
template <class T, class Policy>
T skewness(T v, T delta, const Policy& pol)
{
BOOST_MATH_STD_USING
if ((boost::math::isinf)(v))
{
return 0;
}
if(delta == 0)
{ // == Student's t
return 0;
}
T mean = boost::math::detail::mean(v, delta, pol);
T l2 = delta * delta;
T var = ((l2 + 1) * v) / (v - 2) - mean * mean;
T result = -2 * var;
result += v * (l2 + 2 * v - 3) / ((v - 3) * (v - 2));
result *= mean;
result /= pow(var, T(1.5f));
return result;
}
template <class T, class Policy>
T kurtosis_excess(T v, T delta, const Policy& pol)
{
BOOST_MATH_STD_USING
if ((boost::math::isinf)(v))
{
return 3;
}
if (delta == 0)
{ // == Student's t
return 3;
}
T mean = boost::math::detail::mean(v, delta, pol);
T l2 = delta * delta;
T var = ((l2 + 1) * v) / (v - 2) - mean * mean;
T result = -3 * var;
result += v * (l2 * (v + 1) + 3 * (3 * v - 5)) / ((v - 3) * (v - 2));
result *= -mean * mean;
result += v * v * (l2 * l2 + 6 * l2 + 3) / ((v - 4) * (v - 2));
result /= var * var;
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
{
t_degrees_of_freedom_finder(
RealType delta_, RealType x_, RealType p_, bool c)
: delta(delta_), x(x_), p(p_), comp(c) {}
RealType operator()(const RealType& v)
{
non_central_t_distribution<RealType, Policy> d(v, delta);
return comp ?
p - cdf(complement(d, x))
: cdf(d, x) - p;
}
private:
RealType delta;
RealType x;
RealType p;
bool comp;
};
template <class RealType, class Policy>
inline RealType find_t_degrees_of_freedom(
RealType delta, RealType x, RealType p, RealType q, const Policy& pol)
{
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:
//
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());
}
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>();
//
// Pick an initial guess:
//
RealType guess = 200;
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;
if(max_iter >= policies::get_max_root_iterations<Policy>())
{
return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
" or there is no answer to problem. Current best guess is %1%", result, Policy());
}
return result;
}
template <class RealType, class Policy>
struct t_non_centrality_finder
{
t_non_centrality_finder(
RealType v_, RealType x_, RealType p_, bool c)
: v(v_), x(x_), p(p_), comp(c) {}
RealType operator()(const RealType& delta)
{
non_central_t_distribution<RealType, Policy> d(v, delta);
return comp ?
p - cdf(complement(d, x))
: cdf(d, x) - p;
}
private:
RealType v;
RealType x;
RealType p;
bool comp;
};
template <class RealType, class Policy>
inline RealType find_t_non_centrality(
RealType v, RealType x, RealType p, RealType q, const Policy& pol)
{
const char* function = "non_central_t<%1%>::find_t_non_centrality";
if((p == 0) || (q == 0))
{
//
// Can't do a thing if one of p and q is zero:
//
return policies::raise_evaluation_error<RealType>(function,
"Can't find non-centrality parameter when the probability is 0 or 1, only possible answer is %1%",
RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy());
}
t_non_centrality_finder<RealType, Policy> f(v, 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>();
//
// Pick an initial guess that we know is the right side of
// zero:
//
RealType guess;
if(f(0) < 0)
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;
if(max_iter >= policies::get_max_root_iterations<Policy>())
{
return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
" or there is no answer to problem. Current best guess is %1%", result, Policy());
}
return result;
}
#endif
} // namespace detail ======================================================================
template <class RealType = double, class Policy = policies::policy<> >
class non_central_t_distribution
{
public:
typedef RealType value_type;
typedef Policy policy_type;
non_central_t_distribution(RealType v_, RealType lambda) : v(v_), ncp(lambda)
{
const char* function = "boost::math::non_central_t_distribution<%1%>::non_central_t_distribution(%1%,%1%)";
RealType r;
detail::check_df_gt0_to_inf(
function,
v, &r, Policy());
detail::check_finite(
function,
lambda,
&r,
Policy());
} // non_central_t_distribution constructor.
RealType degrees_of_freedom() const
{ // Private data getter function.
return v;
}
RealType non_centrality() const
{ // 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::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),
forwarding_policy());
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
result,
function);
}
template <class A, class B, class C>
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::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),
forwarding_policy());
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
result,
function);
}
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::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),
forwarding_policy());
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
result,
function);
}
template <class A, class B, class C>
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::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),
forwarding_policy());
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
result,
function);
}
#endif
private:
// Data member, initialized by constructor.
RealType v; // degrees of freedom
RealType ncp; // non-centrality parameter
}; // template <class RealType, class Policy> class non_central_t_distribution
typedef non_central_t_distribution<double> non_central_t; // Reserved name of type double.
#ifdef __cpp_deduction_guides
template <class RealType>
non_central_t_distribution(RealType,RealType)->non_central_t_distribution<typename boost::math::tools::promote_args<RealType>::type>;
#endif
// Non-member functions to give properties of the distribution.
template <class RealType, class Policy>
inline const std::pair<RealType, RealType> range(const non_central_t_distribution<RealType, Policy>& /* dist */)
{ // Range of permissible values for random variable k.
using boost::math::tools::max_value;
return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>());
}
template <class RealType, class Policy>
inline const std::pair<RealType, RealType> support(const non_central_t_distribution<RealType, Policy>& /* dist */)
{ // Range of supported values for random variable k.
// This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
using boost::math::tools::max_value;
return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>());
}
template <class RealType, class Policy>
inline RealType mode(const non_central_t_distribution<RealType, Policy>& dist)
{ // mode.
static const char* function = "mode(non_central_t_distribution<%1%> const&)";
RealType v = dist.degrees_of_freedom();
RealType l = dist.non_centrality();
RealType r;
if(!detail::check_df_gt0_to_inf(
function,
v, &r, Policy())
||
!detail::check_finite(
function,
l,
&r,
Policy()))
return (RealType)r;
BOOST_MATH_STD_USING
RealType m = v < 3 ? 0 : detail::mean(v, l, Policy());
RealType var = v < 4 ? 1 : detail::variance(v, l, Policy());
return detail::generic_find_mode(
dist,
m,
function,
sqrt(var));
}
template <class RealType, class Policy>
inline RealType mean(const non_central_t_distribution<RealType, Policy>& dist)
{
BOOST_MATH_STD_USING
const char* function = "mean(const non_central_t_distribution<%1%>&)";
typedef typename policies::evaluation<RealType, Policy>::type value_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
RealType v = dist.degrees_of_freedom();
RealType l = dist.non_centrality();
RealType r;
if(!detail::check_df_gt0_to_inf(
function,
v, &r, Policy())
||
!detail::check_finite(
function,
l,
&r,
Policy()))
return (RealType)r;
if(v <= 1)
return policies::raise_domain_error<RealType>(
function,
"The non-central t distribution has no defined mean for degrees of freedom <= 1: got v=%1%.", v, Policy());
// return l * sqrt(v / 2) * tgamma_delta_ratio((v - 1) * 0.5f, RealType(0.5f));
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
detail::mean(static_cast<value_type>(v), static_cast<value_type>(l), forwarding_policy()), function);
} // mean
template <class RealType, class Policy>
inline RealType variance(const non_central_t_distribution<RealType, Policy>& dist)
{ // variance.
const char* function = "variance(const non_central_t_distribution<%1%>&)";
typedef typename policies::evaluation<RealType, Policy>::type value_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
BOOST_MATH_STD_USING
RealType v = dist.degrees_of_freedom();
RealType l = dist.non_centrality();
RealType r;
if(!detail::check_df_gt0_to_inf(
function,
v, &r, Policy())
||
!detail::check_finite(
function,
l,
&r,
Policy()))
return (RealType)r;
if(v <= 2)
return policies::raise_domain_error<RealType>(
function,
"The non-central t distribution has no defined variance for degrees of freedom <= 2: got v=%1%.", v, Policy());
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
detail::variance(static_cast<value_type>(v), static_cast<value_type>(l), forwarding_policy()), function);
}
// RealType standard_deviation(const non_central_t_distribution<RealType, Policy>& dist)
// standard_deviation provided by derived accessors.
template <class RealType, class Policy>
inline RealType skewness(const non_central_t_distribution<RealType, Policy>& dist)
{ // skewness = sqrt(l).
const char* function = "skewness(const non_central_t_distribution<%1%>&)";
typedef typename policies::evaluation<RealType, Policy>::type value_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
RealType v = dist.degrees_of_freedom();
RealType l = dist.non_centrality();
RealType r;
if(!detail::check_df_gt0_to_inf(
function,
v, &r, Policy())
||
!detail::check_finite(
function,
l,
&r,
Policy()))
return (RealType)r;
if(v <= 3)
return policies::raise_domain_error<RealType>(
function,
"The non-central t distribution has no defined skewness for degrees of freedom <= 3: got v=%1%.", v, Policy());;
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
detail::skewness(static_cast<value_type>(v), static_cast<value_type>(l), forwarding_policy()), function);
}
template <class RealType, class Policy>
inline RealType kurtosis_excess(const non_central_t_distribution<RealType, Policy>& dist)
{
const char* function = "kurtosis_excess(const non_central_t_distribution<%1%>&)";
typedef typename policies::evaluation<RealType, Policy>::type value_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
RealType v = dist.degrees_of_freedom();
RealType l = dist.non_centrality();
RealType r;
if(!detail::check_df_gt0_to_inf(
function,
v, &r, Policy())
||
!detail::check_finite(
function,
l,
&r,
Policy()))
return (RealType)r;
if(v <= 4)
return policies::raise_domain_error<RealType>(
function,
"The non-central t distribution has no defined kurtosis for degrees of freedom <= 4: got v=%1%.", v, Policy());;
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
detail::kurtosis_excess(static_cast<value_type>(v), static_cast<value_type>(l), forwarding_policy()), function);
} // kurtosis_excess
template <class RealType, class Policy>
inline RealType kurtosis(const non_central_t_distribution<RealType, Policy>& dist)
{
return kurtosis_excess(dist) + 3;
}
template <class RealType, class Policy>
inline RealType pdf(const non_central_t_distribution<RealType, Policy>& dist, const RealType& t)
{ // Probability Density/Mass Function.
const char* function = "pdf(non_central_t_distribution<%1%>, %1%)";
typedef typename policies::evaluation<RealType, Policy>::type value_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
RealType v = dist.degrees_of_freedom();
RealType l = dist.non_centrality();
RealType r;
if(!detail::check_df_gt0_to_inf(
function,
v, &r, Policy())
||
!detail::check_finite(
function,
l,
&r,
Policy())
||
!detail::check_x(
function,
t,
&r,
Policy()))
return (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),
static_cast<value_type>(t),
Policy()),
function);
} // pdf
template <class RealType, class Policy>
RealType cdf(const non_central_t_distribution<RealType, Policy>& dist, const RealType& x)
{
const char* function = "boost::math::cdf(non_central_t_distribution<%1%>&, %1%)";
// was const char* function = "boost::math::non_central_t_distribution<%1%>::cdf(%1%)";
typedef typename policies::evaluation<RealType, Policy>::type value_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
RealType v = dist.degrees_of_freedom();
RealType l = dist.non_centrality();
RealType r;
if(!detail::check_df_gt0_to_inf(
function,
v, &r, Policy())
||
!detail::check_finite(
function,
l,
&r,
Policy())
||
!detail::check_x(
function,
x,
&r,
Policy()))
return (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);
//return cdf(normal_distribution<RealType, Policy>(l, 1), x);
}
if(l == 0)
{ // NO non-centrality, so use Student's t instead.
return cdf(students_t_distribution<RealType, Policy>(v), x);
}
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
detail::non_central_t_cdf(
static_cast<value_type>(v),
static_cast<value_type>(l),
static_cast<value_type>(x),
false, Policy()),
function);
} // cdf
template <class RealType, class Policy>
RealType cdf(const complemented2_type<non_central_t_distribution<RealType, Policy>, RealType>& c)
{ // Complemented Cumulative Distribution Function
// was const char* function = "boost::math::non_central_t_distribution<%1%>::cdf(%1%)";
const char* function = "boost::math::cdf(const complement(non_central_t_distribution<%1%>&), %1%)";
typedef typename policies::evaluation<RealType, Policy>::type value_type;
typedef typename policies::normalise<
Policy,
policies::promote_float<false>,
policies::promote_double<false>,
policies::discrete_quantile<>,
policies::assert_undefined<> >::type forwarding_policy;
non_central_t_distribution<RealType, Policy> const& dist = c.dist;
RealType x = c.param;
RealType v = dist.degrees_of_freedom();
RealType l = dist.non_centrality(); // aka delta
RealType r;
if(!detail::check_df_gt0_to_inf(
function,
v, &r, Policy())
||
!detail::check_finite(
function,
l,
&r,
Policy())
||
!detail::check_x(
function,
x,
&r,
Policy()))
return (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);
return cdf(complement(n, x));
}
if(l == 0)
{ // zero non-centrality so use Student's t distribution.
return cdf(complement(students_t_distribution<RealType, Policy>(v), x));
}
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
detail::non_central_t_cdf(
static_cast<value_type>(v),
static_cast<value_type>(l),
static_cast<value_type>(x),
true, Policy()),
function);
} // ccdf
template <class RealType, class Policy>
inline RealType quantile(const non_central_t_distribution<RealType, Policy>& dist, const RealType& p)
{ // Quantile (or Percent Point) function.
static const char* function = "quantile(const non_central_t_distribution<%1%>, %1%)";
RealType v = dist.degrees_of_freedom();
RealType l = dist.non_centrality();
return detail::non_central_t_quantile(function, v, l, p, RealType(1-p), Policy());
} // quantile
template <class RealType, class Policy>
inline RealType quantile(const complemented2_type<non_central_t_distribution<RealType, Policy>, RealType>& c)
{ // Quantile (or Percent Point) function.
static const char* function = "quantile(const complement(non_central_t_distribution<%1%>, %1%))";
non_central_t_distribution<RealType, Policy> const& dist = c.dist;
RealType q = c.param;
RealType v = dist.degrees_of_freedom();
RealType l = dist.non_centrality();
return detail::non_central_t_quantile(function, v, l, RealType(1-q), q, Policy());
} // quantile complement.
} // namespace math
} // namespace boost
// This include must be at the end, *after* the accessors
// for this distribution have been defined, in order to
// keep compilers that support two-phase lookup happy.
#include <boost/math/distributions/detail/derived_accessors.hpp>
#endif // BOOST_MATH_SPECIAL_NON_CENTRAL_T_HPP