Mercurial > hg > sv-dependency-builds
diff any/include/boost/math/distributions/non_central_f.hpp @ 160:cff480c41f97
Add some cross-platform Boost headers
author | Chris Cannam <cannam@all-day-breakfast.com> |
---|---|
date | Sat, 16 Feb 2019 16:31:25 +0000 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/any/include/boost/math/distributions/non_central_f.hpp Sat Feb 16 16:31:25 2019 +0000 @@ -0,0 +1,410 @@ +// boost\math\distributions\non_central_f.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_F_HPP +#define BOOST_MATH_SPECIAL_NON_CENTRAL_F_HPP + +#include <boost/math/distributions/non_central_beta.hpp> +#include <boost/math/distributions/detail/generic_mode.hpp> +#include <boost/math/special_functions/pow.hpp> + +namespace boost +{ + namespace math + { + template <class RealType = double, class Policy = policies::policy<> > + class non_central_f_distribution + { + public: + typedef RealType value_type; + typedef Policy policy_type; + + non_central_f_distribution(RealType v1_, RealType v2_, RealType lambda) : v1(v1_), v2(v2_), ncp(lambda) + { + const char* function = "boost::math::non_central_f_distribution<%1%>::non_central_f_distribution(%1%,%1%)"; + RealType r; + detail::check_df( + function, + v1, &r, Policy()); + detail::check_df( + function, + v2, &r, Policy()); + detail::check_non_centrality( + function, + lambda, + &r, + Policy()); + } // non_central_f_distribution constructor. + + RealType degrees_of_freedom1()const + { + return v1; + } + RealType degrees_of_freedom2()const + { + return v2; + } + RealType non_centrality() const + { // Private data getter function. + return ncp; + } + private: + // Data member, initialized by constructor. + RealType v1; // alpha. + RealType v2; // beta. + RealType ncp; // non-centrality parameter + }; // template <class RealType, class Policy> class non_central_f_distribution + + typedef non_central_f_distribution<double> non_central_f; // Reserved name of type double. + + // Non-member functions to give properties of the distribution. + + template <class RealType, class Policy> + inline const std::pair<RealType, RealType> range(const non_central_f_distribution<RealType, Policy>& /* dist */) + { // Range of permissible values for random variable k. + using boost::math::tools::max_value; + return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); + } + + template <class RealType, class Policy> + inline const std::pair<RealType, RealType> support(const non_central_f_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>(static_cast<RealType>(0), max_value<RealType>()); + } + + template <class RealType, class Policy> + inline RealType mean(const non_central_f_distribution<RealType, Policy>& dist) + { + const char* function = "mean(non_central_f_distribution<%1%> const&)"; + RealType v1 = dist.degrees_of_freedom1(); + RealType v2 = dist.degrees_of_freedom2(); + RealType l = dist.non_centrality(); + RealType r; + if(!detail::check_df( + function, + v1, &r, Policy()) + || + !detail::check_df( + function, + v2, &r, Policy()) + || + !detail::check_non_centrality( + function, + l, + &r, + Policy())) + return r; + if(v2 <= 2) + return policies::raise_domain_error( + function, + "Second degrees of freedom parameter was %1%, but must be > 2 !", + v2, Policy()); + return v2 * (v1 + l) / (v1 * (v2 - 2)); + } // mean + + template <class RealType, class Policy> + inline RealType mode(const non_central_f_distribution<RealType, Policy>& dist) + { // mode. + static const char* function = "mode(non_central_chi_squared_distribution<%1%> const&)"; + + RealType n = dist.degrees_of_freedom1(); + RealType m = dist.degrees_of_freedom2(); + RealType l = dist.non_centrality(); + RealType r; + if(!detail::check_df( + function, + n, &r, Policy()) + || + !detail::check_df( + function, + m, &r, Policy()) + || + !detail::check_non_centrality( + function, + l, + &r, + Policy())) + return r; + RealType guess = m > 2 ? RealType(m * (n + l) / (n * (m - 2))) : RealType(1); + return detail::generic_find_mode( + dist, + guess, + function); + } + + template <class RealType, class Policy> + inline RealType variance(const non_central_f_distribution<RealType, Policy>& dist) + { // variance. + const char* function = "variance(non_central_f_distribution<%1%> const&)"; + RealType n = dist.degrees_of_freedom1(); + RealType m = dist.degrees_of_freedom2(); + RealType l = dist.non_centrality(); + RealType r; + if(!detail::check_df( + function, + n, &r, Policy()) + || + !detail::check_df( + function, + m, &r, Policy()) + || + !detail::check_non_centrality( + function, + l, + &r, + Policy())) + return r; + if(m <= 4) + return policies::raise_domain_error( + function, + "Second degrees of freedom parameter was %1%, but must be > 4 !", + m, Policy()); + RealType result = 2 * m * m * ((n + l) * (n + l) + + (m - 2) * (n + 2 * l)); + result /= (m - 4) * (m - 2) * (m - 2) * n * n; + return result; + } + + // RealType standard_deviation(const non_central_f_distribution<RealType, Policy>& dist) + // standard_deviation provided by derived accessors. + + template <class RealType, class Policy> + inline RealType skewness(const non_central_f_distribution<RealType, Policy>& dist) + { // skewness = sqrt(l). + const char* function = "skewness(non_central_f_distribution<%1%> const&)"; + BOOST_MATH_STD_USING + RealType n = dist.degrees_of_freedom1(); + RealType m = dist.degrees_of_freedom2(); + RealType l = dist.non_centrality(); + RealType r; + if(!detail::check_df( + function, + n, &r, Policy()) + || + !detail::check_df( + function, + m, &r, Policy()) + || + !detail::check_non_centrality( + function, + l, + &r, + Policy())) + return r; + if(m <= 6) + return policies::raise_domain_error( + function, + "Second degrees of freedom parameter was %1%, but must be > 6 !", + m, Policy()); + RealType result = 2 * constants::root_two<RealType>(); + result *= sqrt(m - 4); + result *= (n * (m + n - 2) *(m + 2 * n - 2) + + 3 * (m + n - 2) * (m + 2 * n - 2) * l + + 6 * (m + n - 2) * l * l + 2 * l * l * l); + result /= (m - 6) * pow(n * (m + n - 2) + 2 * (m + n - 2) * l + l * l, RealType(1.5f)); + return result; + } + + template <class RealType, class Policy> + inline RealType kurtosis_excess(const non_central_f_distribution<RealType, Policy>& dist) + { + const char* function = "kurtosis_excess(non_central_f_distribution<%1%> const&)"; + BOOST_MATH_STD_USING + RealType n = dist.degrees_of_freedom1(); + RealType m = dist.degrees_of_freedom2(); + RealType l = dist.non_centrality(); + RealType r; + if(!detail::check_df( + function, + n, &r, Policy()) + || + !detail::check_df( + function, + m, &r, Policy()) + || + !detail::check_non_centrality( + function, + l, + &r, + Policy())) + return r; + if(m <= 8) + return policies::raise_domain_error( + function, + "Second degrees of freedom parameter was %1%, but must be > 8 !", + m, Policy()); + RealType l2 = l * l; + RealType l3 = l2 * l; + RealType l4 = l2 * l2; + RealType result = (3 * (m - 4) * (n * (m + n - 2) + * (4 * (m - 2) * (m - 2) + + (m - 2) * (m + 10) * n + + (10 + m) * n * n) + + 4 * (m + n - 2) * (4 * (m - 2) * (m - 2) + + (m - 2) * (10 + m) * n + + (10 + m) * n * n) * l + 2 * (10 + m) + * (m + n - 2) * (2 * m + 3 * n - 4) * l2 + + 4 * (10 + m) * (-2 + m + n) * l3 + + (10 + m) * l4)) + / + ((-8 + m) * (-6 + m) * boost::math::pow<2>(n * (-2 + m + n) + + 2 * (-2 + m + n) * l + l2)); + return result; + } // kurtosis_excess + + template <class RealType, class Policy> + inline RealType kurtosis(const non_central_f_distribution<RealType, Policy>& dist) + { + return kurtosis_excess(dist) + 3; + } + + template <class RealType, class Policy> + inline RealType pdf(const non_central_f_distribution<RealType, Policy>& dist, const RealType& x) + { // Probability Density/Mass Function. + 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 alpha = dist.degrees_of_freedom1() / 2; + value_type beta = dist.degrees_of_freedom2() / 2; + value_type y = x * alpha / beta; + value_type r = pdf(boost::math::non_central_beta_distribution<value_type, forwarding_policy>(alpha, beta, dist.non_centrality()), y / (1 + y)); + return policies::checked_narrowing_cast<RealType, forwarding_policy>( + r * (dist.degrees_of_freedom1() / dist.degrees_of_freedom2()) / ((1 + y) * (1 + y)), + "pdf(non_central_f_distribution<%1%>, %1%)"); + } // pdf + + template <class RealType, class Policy> + RealType cdf(const non_central_f_distribution<RealType, Policy>& dist, const RealType& x) + { + const char* function = "cdf(const non_central_f_distribution<%1%>&, %1%)"; + RealType r; + if(!detail::check_df( + function, + dist.degrees_of_freedom1(), &r, Policy()) + || + !detail::check_df( + function, + dist.degrees_of_freedom2(), &r, Policy()) + || + !detail::check_non_centrality( + function, + dist.non_centrality(), + &r, + Policy())) + return r; + + if((x < 0) || !(boost::math::isfinite)(x)) + { + return policies::raise_domain_error<RealType>( + function, "Random Variable parameter was %1%, but must be > 0 !", x, Policy()); + } + + RealType alpha = dist.degrees_of_freedom1() / 2; + RealType beta = dist.degrees_of_freedom2() / 2; + RealType y = x * alpha / beta; + RealType c = y / (1 + y); + RealType cp = 1 / (1 + y); + // + // To ensure accuracy, we pass both x and 1-x to the + // non-central beta cdf routine, this ensures accuracy + // even when we compute x to be ~ 1: + // + r = detail::non_central_beta_cdf(c, cp, alpha, beta, + dist.non_centrality(), false, Policy()); + return r; + } // cdf + + template <class RealType, class Policy> + RealType cdf(const complemented2_type<non_central_f_distribution<RealType, Policy>, RealType>& c) + { // Complemented Cumulative Distribution Function + const char* function = "cdf(complement(const non_central_f_distribution<%1%>&, %1%))"; + RealType r; + if(!detail::check_df( + function, + c.dist.degrees_of_freedom1(), &r, Policy()) + || + !detail::check_df( + function, + c.dist.degrees_of_freedom2(), &r, Policy()) + || + !detail::check_non_centrality( + function, + c.dist.non_centrality(), + &r, + Policy())) + return r; + + if((c.param < 0) || !(boost::math::isfinite)(c.param)) + { + return policies::raise_domain_error<RealType>( + function, "Random Variable parameter was %1%, but must be > 0 !", c.param, Policy()); + } + + RealType alpha = c.dist.degrees_of_freedom1() / 2; + RealType beta = c.dist.degrees_of_freedom2() / 2; + RealType y = c.param * alpha / beta; + RealType x = y / (1 + y); + RealType cx = 1 / (1 + y); + // + // To ensure accuracy, we pass both x and 1-x to the + // non-central beta cdf routine, this ensures accuracy + // even when we compute x to be ~ 1: + // + r = detail::non_central_beta_cdf(x, cx, alpha, beta, + c.dist.non_centrality(), true, Policy()); + return r; + } // ccdf + + template <class RealType, class Policy> + inline RealType quantile(const non_central_f_distribution<RealType, Policy>& dist, const RealType& p) + { // Quantile (or Percent Point) function. + RealType alpha = dist.degrees_of_freedom1() / 2; + RealType beta = dist.degrees_of_freedom2() / 2; + RealType x = quantile(boost::math::non_central_beta_distribution<RealType, Policy>(alpha, beta, dist.non_centrality()), p); + if(x == 1) + return policies::raise_overflow_error<RealType>( + "quantile(const non_central_f_distribution<%1%>&, %1%)", + "Result of non central F quantile is too large to represent.", + Policy()); + return (x / (1 - x)) * (dist.degrees_of_freedom2() / dist.degrees_of_freedom1()); + } // quantile + + template <class RealType, class Policy> + inline RealType quantile(const complemented2_type<non_central_f_distribution<RealType, Policy>, RealType>& c) + { // Quantile (or Percent Point) function. + RealType alpha = c.dist.degrees_of_freedom1() / 2; + RealType beta = c.dist.degrees_of_freedom2() / 2; + RealType x = quantile(complement(boost::math::non_central_beta_distribution<RealType, Policy>(alpha, beta, c.dist.non_centrality()), c.param)); + if(x == 1) + return policies::raise_overflow_error<RealType>( + "quantile(complement(const non_central_f_distribution<%1%>&, %1%))", + "Result of non central F quantile is too large to represent.", + Policy()); + return (x / (1 - x)) * (c.dist.degrees_of_freedom2() / c.dist.degrees_of_freedom1()); + } // 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_F_HPP + + +