boost/math/distributions/non_central_f.hpp
// 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>(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>(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; return detail::generic_find_mode( dist, m * (n + l) / (n * (m - 2)), 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