| // 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; |
| 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 |
| |
| |
| |