math/distributions/non_central_f.hpp

0001 // boost\math\distributions\non_central_f.hpp
0002
0003 // Copyright John Maddock 2008.
0004
0005 // Use, modification and distribution are subject to the
0006 // Boost Software License, Version 1.0.
0007 // (See accompanying file LICENSE_1_0.txt
0008 // or copy at http://www.boost.org/LICENSE_1_0.txt)
0009
0010 #ifndef BOOST_MATH_SPECIAL_NON_CENTRAL_F_HPP
0011 #define BOOST_MATH_SPECIAL_NON_CENTRAL_F_HPP
0012
0013 #include <boost/math/distributions/non_central_beta.hpp>
0014 #include <boost/math/distributions/detail/generic_mode.hpp>
0015 #include <boost/math/special_functions/pow.hpp>
0016
0017 namespace boost
0018 {
0019    namespace math
0020    {
0021       template <class RealType = double, class Policy = policies::policy<> >
0022       class non_central_f_distribution
0023       {
0024       public:
0025          typedef RealType value_type;
0026          typedef Policy policy_type;
0027
0028          non_central_f_distribution(RealType v1_, RealType v2_, RealType lambda) : v1(v1_), v2(v2_), ncp(lambda)
0029          {
0030             const char* function = "boost::math::non_central_f_distribution<%1%>::non_central_f_distribution(%1%,%1%)";
0031             RealType r;
0032             detail::check_df(
0033                function,
0034                v1, &r, Policy());
0035             detail::check_df(
0036                function,
0037                v2, &r, Policy());
0038             detail::check_non_centrality(
0039                function,
0040                lambda,
0041                &r,
0042                Policy());
0043          } // non_central_f_distribution constructor.
0044
0045          RealType degrees_of_freedom1()const
0046          {
0047             return v1;
0048          }
0049          RealType degrees_of_freedom2()const
0050          {
0051             return v2;
0052          }
0053          RealType non_centrality() const
0054          { // Private data getter function.
0055             return ncp;
0056          }
0057       private:
0058          // Data member, initialized by constructor.
0059          RealType v1;   // alpha.
0060          RealType v2;   // beta.
0061          RealType ncp; // non-centrality parameter
0062       }; // template <class RealType, class Policy> class non_central_f_distribution
0063
0064       typedef non_central_f_distribution<double> non_central_f; // Reserved name of type double.
0065
0066       #ifdef __cpp_deduction_guides
0067       template <class RealType>
0068       non_central_f_distribution(RealType,RealType,RealType)->non_central_f_distribution<typename boost::math::tools::promote_args<RealType>::type>;
0069       #endif
0070
0071       // Non-member functions to give properties of the distribution.
0072
0073       template <class RealType, class Policy>
0074       inline const std::pair<RealType, RealType> range(const non_central_f_distribution<RealType, Policy>& /* dist */)
0075       { // Range of permissible values for random variable k.
0076          using boost::math::tools::max_value;
0077          return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());
0078       }
0079
0080       template <class RealType, class Policy>
0081       inline const std::pair<RealType, RealType> support(const non_central_f_distribution<RealType, Policy>& /* dist */)
0082       { // Range of supported values for random variable k.
0083          // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
0084          using boost::math::tools::max_value;
0085          return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());
0086       }
0087
0088       template <class RealType, class Policy>
0089       inline RealType mean(const non_central_f_distribution<RealType, Policy>& dist)
0090       {
0091          const char* function = "mean(non_central_f_distribution<%1%> const&)";
0092          RealType v1 = dist.degrees_of_freedom1();
0093          RealType v2 = dist.degrees_of_freedom2();
0094          RealType l = dist.non_centrality();
0095          RealType r;
0096          if(!detail::check_df(
0097             function,
0098             v1, &r, Policy())
0099                ||
0100             !detail::check_df(
0101                function,
0102                v2, &r, Policy())
0103                ||
0104             !detail::check_non_centrality(
0105                function,
0106                l,
0107                &r,
0108                Policy()))
0109                   return r;
0110          if(v2 <= 2)
0111             return policies::raise_domain_error(
0112                function,
0113                "Second degrees of freedom parameter was %1%, but must be > 2 !",
0114                v2, Policy());
0115          return v2 * (v1 + l) / (v1 * (v2 - 2));
0116       } // mean
0117
0118       template <class RealType, class Policy>
0119       inline RealType mode(const non_central_f_distribution<RealType, Policy>& dist)
0120       { // mode.
0121          static const char* function = "mode(non_central_chi_squared_distribution<%1%> const&)";
0122
0123          RealType n = dist.degrees_of_freedom1();
0124          RealType m = dist.degrees_of_freedom2();
0125          RealType l = dist.non_centrality();
0126          RealType r;
0127          if(!detail::check_df(
0128             function,
0129             n, &r, Policy())
0130                ||
0131             !detail::check_df(
0132                function,
0133                m, &r, Policy())
0134                ||
0135             !detail::check_non_centrality(
0136                function,
0137                l,
0138                &r,
0139                Policy()))
0140                   return r;
0141          RealType guess = m > 2 ? RealType(m * (n + l) / (n * (m - 2))) : RealType(1);
0142          return detail::generic_find_mode(
0143             dist,
0144             guess,
0145             function);
0146       }
0147
0148       template <class RealType, class Policy>
0149       inline RealType variance(const non_central_f_distribution<RealType, Policy>& dist)
0150       { // variance.
0151          const char* function = "variance(non_central_f_distribution<%1%> const&)";
0152          RealType n = dist.degrees_of_freedom1();
0153          RealType m = dist.degrees_of_freedom2();
0154          RealType l = dist.non_centrality();
0155          RealType r;
0156          if(!detail::check_df(
0157             function,
0158             n, &r, Policy())
0159                ||
0160             !detail::check_df(
0161                function,
0162                m, &r, Policy())
0163                ||
0164             !detail::check_non_centrality(
0165                function,
0166                l,
0167                &r,
0168                Policy()))
0169                   return r;
0170          if(m <= 4)
0171             return policies::raise_domain_error(
0172                function,
0173                "Second degrees of freedom parameter was %1%, but must be > 4 !",
0174                m, Policy());
0175          RealType result = 2 * m * m * ((n + l) * (n + l)
0176             + (m - 2) * (n + 2 * l));
0177          result /= (m - 4) * (m - 2) * (m - 2) * n * n;
0178          return result;
0179       }
0180
0181       // RealType standard_deviation(const non_central_f_distribution<RealType, Policy>& dist)
0182       // standard_deviation provided by derived accessors.
0183
0184       template <class RealType, class Policy>
0185       inline RealType skewness(const non_central_f_distribution<RealType, Policy>& dist)
0186       { // skewness = sqrt(l).
0187          const char* function = "skewness(non_central_f_distribution<%1%> const&)";
0188          BOOST_MATH_STD_USING
0189          RealType n = dist.degrees_of_freedom1();
0190          RealType m = dist.degrees_of_freedom2();
0191          RealType l = dist.non_centrality();
0192          RealType r;
0193          if(!detail::check_df(
0194             function,
0195             n, &r, Policy())
0196                ||
0197             !detail::check_df(
0198                function,
0199                m, &r, Policy())
0200                ||
0201             !detail::check_non_centrality(
0202                function,
0203                l,
0204                &r,
0205                Policy()))
0206                   return r;
0207          if(m <= 6)
0208             return policies::raise_domain_error(
0209                function,
0210                "Second degrees of freedom parameter was %1%, but must be > 6 !",
0211                m, Policy());
0212          RealType result = 2 * constants::root_two<RealType>();
0213          result *= sqrt(m - 4);
0214          result *= (n * (m + n - 2) *(m + 2 * n - 2)
0215             + 3 * (m + n - 2) * (m + 2 * n - 2) * l
0216             + 6 * (m + n - 2) * l * l + 2 * l * l * l);
0217          result /= (m - 6) * pow(n * (m + n - 2) + 2 * (m + n - 2) * l + l * l, RealType(1.5f));
0218          return result;
0219       }
0220
0221       template <class RealType, class Policy>
0222       inline RealType kurtosis_excess(const non_central_f_distribution<RealType, Policy>& dist)
0223       {
0224          const char* function = "kurtosis_excess(non_central_f_distribution<%1%> const&)";
0225          BOOST_MATH_STD_USING
0226          RealType n = dist.degrees_of_freedom1();
0227          RealType m = dist.degrees_of_freedom2();
0228          RealType l = dist.non_centrality();
0229          RealType r;
0230          if(!detail::check_df(
0231             function,
0232             n, &r, Policy())
0233                ||
0234             !detail::check_df(
0235                function,
0236                m, &r, Policy())
0237                ||
0238             !detail::check_non_centrality(
0239                function,
0240                l,
0241                &r,
0242                Policy()))
0243                   return r;
0244          if(m <= 8)
0245             return policies::raise_domain_error(
0246                function,
0247                "Second degrees of freedom parameter was %1%, but must be > 8 !",
0248                m, Policy());
0249          RealType l2 = l * l;
0250          RealType l3 = l2 * l;
0251          RealType l4 = l2 * l2;
0252          RealType result = (3 * (m - 4) * (n * (m + n - 2)
0253             * (4 * (m - 2) * (m - 2)
0254             + (m - 2) * (m + 10) * n
0255             + (10 + m) * n * n)
0256             + 4 * (m + n - 2) * (4 * (m - 2) * (m - 2)
0257             + (m - 2) * (10 + m) * n
0258             + (10 + m) * n * n) * l + 2 * (10 + m)
0259             * (m + n - 2) * (2 * m + 3 * n - 4) * l2
0260             + 4 * (10 + m) * (-2 + m + n) * l3
0261             + (10 + m) * l4))
0262             /
0263             ((-8 + m) * (-6 + m) * boost::math::pow<2>(n * (-2 + m + n)
0264             + 2 * (-2 + m + n) * l + l2));
0265             return result;
0266       } // kurtosis_excess
0267
0268       template <class RealType, class Policy>
0269       inline RealType kurtosis(const non_central_f_distribution<RealType, Policy>& dist)
0270       {
0271          return kurtosis_excess(dist) + 3;
0272       }
0273
0274       template <class RealType, class Policy>
0275       inline RealType pdf(const non_central_f_distribution<RealType, Policy>& dist, const RealType& x)
0276       { // Probability Density/Mass Function.
0277          typedef typename policies::evaluation<RealType, Policy>::type value_type;
0278          typedef typename policies::normalise<
0279             Policy,
0280             policies::promote_float<false>,
0281             policies::promote_double<false>,
0282             policies::discrete_quantile<>,
0283             policies::assert_undefined<> >::type forwarding_policy;
0284
0285          value_type alpha = dist.degrees_of_freedom1() / 2;
0286          value_type beta = dist.degrees_of_freedom2() / 2;
0287          value_type y = x * alpha / beta;
0288          value_type r = pdf(boost::math::non_central_beta_distribution<value_type, forwarding_policy>(alpha, beta, dist.non_centrality()), y / (1 + y));
0289          return policies::checked_narrowing_cast<RealType, forwarding_policy>(
0290             r * (dist.degrees_of_freedom1() / dist.degrees_of_freedom2()) / ((1 + y) * (1 + y)),
0291             "pdf(non_central_f_distribution<%1%>, %1%)");
0292       } // pdf
0293
0294       template <class RealType, class Policy>
0295       RealType cdf(const non_central_f_distribution<RealType, Policy>& dist, const RealType& x)
0296       {
0297          const char* function = "cdf(const non_central_f_distribution<%1%>&, %1%)";
0298          RealType r;
0299          if(!detail::check_df(
0300             function,
0301             dist.degrees_of_freedom1(), &r, Policy())
0302                ||
0303             !detail::check_df(
0304                function,
0305                dist.degrees_of_freedom2(), &r, Policy())
0306                ||
0307             !detail::check_non_centrality(
0308                function,
0309                dist.non_centrality(),
0310                &r,
0311                Policy()))
0312                   return r;
0313
0314          if((x < 0) || !(boost::math::isfinite)(x))
0315          {
0316             return policies::raise_domain_error<RealType>(
0317                function, "Random Variable parameter was %1%, but must be > 0 !", x, Policy());
0318          }
0319
0320          RealType alpha = dist.degrees_of_freedom1() / 2;
0321          RealType beta = dist.degrees_of_freedom2() / 2;
0322          RealType y = x * alpha / beta;
0323          RealType c = y / (1 + y);
0324          RealType cp = 1 / (1 + y);
0325          //
0326          // To ensure accuracy, we pass both x and 1-x to the
0327          // non-central beta cdf routine, this ensures accuracy
0328          // even when we compute x to be ~ 1:
0329          //
0330          r = detail::non_central_beta_cdf(c, cp, alpha, beta,
0331             dist.non_centrality(), false, Policy());
0332          return r;
0333       } // cdf
0334
0335       template <class RealType, class Policy>
0336       RealType cdf(const complemented2_type<non_central_f_distribution<RealType, Policy>, RealType>& c)
0337       { // Complemented Cumulative Distribution Function
0338          const char* function = "cdf(complement(const non_central_f_distribution<%1%>&, %1%))";
0339          RealType r;
0340          if(!detail::check_df(
0341             function,
0342             c.dist.degrees_of_freedom1(), &r, Policy())
0343                ||
0344             !detail::check_df(
0345                function,
0346                c.dist.degrees_of_freedom2(), &r, Policy())
0347                ||
0348             !detail::check_non_centrality(
0349                function,
0350                c.dist.non_centrality(),
0351                &r,
0352                Policy()))
0353                   return r;
0354
0355          if((c.param < 0) || !(boost::math::isfinite)(c.param))
0356          {
0357             return policies::raise_domain_error<RealType>(
0358                function, "Random Variable parameter was %1%, but must be > 0 !", c.param, Policy());
0359          }
0360
0361          RealType alpha = c.dist.degrees_of_freedom1() / 2;
0362          RealType beta = c.dist.degrees_of_freedom2() / 2;
0363          RealType y = c.param * alpha / beta;
0364          RealType x = y / (1 + y);
0365          RealType cx = 1 / (1 + y);
0366          //
0367          // To ensure accuracy, we pass both x and 1-x to the
0368          // non-central beta cdf routine, this ensures accuracy
0369          // even when we compute x to be ~ 1:
0370          //
0371          r = detail::non_central_beta_cdf(x, cx, alpha, beta,
0372             c.dist.non_centrality(), true, Policy());
0373          return r;
0374       } // ccdf
0375
0376       template <class RealType, class Policy>
0377       inline RealType quantile(const non_central_f_distribution<RealType, Policy>& dist, const RealType& p)
0378       { // Quantile (or Percent Point) function.
0379          RealType alpha = dist.degrees_of_freedom1() / 2;
0380          RealType beta = dist.degrees_of_freedom2() / 2;
0381          RealType x = quantile(boost::math::non_central_beta_distribution<RealType, Policy>(alpha, beta, dist.non_centrality()), p);
0382          if(x == 1)
0383             return policies::raise_overflow_error<RealType>(
0384                "quantile(const non_central_f_distribution<%1%>&, %1%)",
0385                "Result of non central F quantile is too large to represent.",
0386                Policy());
0387          return (x / (1 - x)) * (dist.degrees_of_freedom2() / dist.degrees_of_freedom1());
0388       } // quantile
0389
0390       template <class RealType, class Policy>
0391       inline RealType quantile(const complemented2_type<non_central_f_distribution<RealType, Policy>, RealType>& c)
0392       { // Quantile (or Percent Point) function.
0393          RealType alpha = c.dist.degrees_of_freedom1() / 2;
0394          RealType beta = c.dist.degrees_of_freedom2() / 2;
0395          RealType x = quantile(complement(boost::math::non_central_beta_distribution<RealType, Policy>(alpha, beta, c.dist.non_centrality()), c.param));
0396          if(x == 1)
0397             return policies::raise_overflow_error<RealType>(
0398                "quantile(complement(const non_central_f_distribution<%1%>&, %1%))",
0399                "Result of non central F quantile is too large to represent.",
0400                Policy());
0401          return (x / (1 - x)) * (c.dist.degrees_of_freedom2() / c.dist.degrees_of_freedom1());
0402       } // quantile complement.
0403
0404    } // namespace math
0405 } // namespace boost
0406
0407 // This include must be at the end, *after* the accessors
0408 // for this distribution have been defined, in order to
0409 // keep compilers that support two-phase lookup happy.
0410 #include <boost/math/distributions/detail/derived_accessors.hpp>
0411
0412 #endif // BOOST_MATH_SPECIAL_NON_CENTRAL_F_HPP
0413
0414
0415