ceph/src/boost/libs/math/include/boost/math/distributions/non_central_f.hpp

   1 // boost\math\distributions\non_central_f.hpp
   2
   3 // Copyright John Maddock 2008.
   4
   5 // Use, modification and distribution are subject to the
   6 // Boost Software License, Version 1.0.
   7 // (See accompanying file LICENSE_1_0.txt
   8 // or copy at http://www.boost.org/LICENSE_1_0.txt)
   9
  10 #ifndef BOOST_MATH_SPECIAL_NON_CENTRAL_F_HPP
  11 #define BOOST_MATH_SPECIAL_NON_CENTRAL_F_HPP
  12
  13 #include <boost/math/distributions/non_central_beta.hpp>
  14 #include <boost/math/distributions/detail/generic_mode.hpp>
  15 #include <boost/math/special_functions/pow.hpp>
  16
  17 namespace boost
  18 {
  19    namespace math
  20    {
  21       template <class RealType = double, class Policy = policies::policy<> >
  22       class non_central_f_distribution
  23       {
  24       public:
  25          typedef RealType value_type;
  26          typedef Policy policy_type;
  27
  28          non_central_f_distribution(RealType v1_, RealType v2_, RealType lambda) : v1(v1_), v2(v2_), ncp(lambda)
  29          {
  30             const char* function = "boost::math::non_central_f_distribution<%1%>::non_central_f_distribution(%1%,%1%)";
  31             RealType r;
  32             detail::check_df(
  33                function,
  34                v1, &r, Policy());
  35             detail::check_df(
  36                function,
  37                v2, &r, Policy());
  38             detail::check_non_centrality(
  39                function,
  40                lambda,
  41                &r,
  42                Policy());
  43          } // non_central_f_distribution constructor.
  44
  45          RealType degrees_of_freedom1()const
  46          {
  47             return v1;
  48          }
  49          RealType degrees_of_freedom2()const
  50          {
  51             return v2;
  52          }
  53          RealType non_centrality() const
  54          { // Private data getter function.
  55             return ncp;
  56          }
  57       private:
  58          // Data member, initialized by constructor.
  59          RealType v1;   // alpha.
  60          RealType v2;   // beta.
  61          RealType ncp; // non-centrality parameter
  62       }; // template <class RealType, class Policy> class non_central_f_distribution
  63
  64       typedef non_central_f_distribution<double> non_central_f; // Reserved name of type double.
  65
  66       // Non-member functions to give properties of the distribution.
  67
  68       template <class RealType, class Policy>
  69       inline const std::pair<RealType, RealType> range(const non_central_f_distribution<RealType, Policy>& /* dist */)
  70       { // Range of permissible values for random variable k.
  71          using boost::math::tools::max_value;
  72          return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());
  73       }
  74
  75       template <class RealType, class Policy>
  76       inline const std::pair<RealType, RealType> support(const non_central_f_distribution<RealType, Policy>& /* dist */)
  77       { // Range of supported values for random variable k.
  78          // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
  79          using boost::math::tools::max_value;
  80          return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());
  81       }
  82
  83       template <class RealType, class Policy>
  84       inline RealType mean(const non_central_f_distribution<RealType, Policy>& dist)
  85       {
  86          const char* function = "mean(non_central_f_distribution<%1%> const&)";
  87          RealType v1 = dist.degrees_of_freedom1();
  88          RealType v2 = dist.degrees_of_freedom2();
  89          RealType l = dist.non_centrality();
  90          RealType r;
  91          if(!detail::check_df(
  92             function,
  93             v1, &r, Policy())
  94                ||
  95             !detail::check_df(
  96                function,
  97                v2, &r, Policy())
  98                ||
  99             !detail::check_non_centrality(
 100                function,
 101                l,
 102                &r,
 103                Policy()))
 104                return r;
 105          if(v2 <= 2)
 106             return policies::raise_domain_error(
 107                function,
 108                "Second degrees of freedom parameter was %1%, but must be > 2 !",
 109                v2, Policy());
 110          return v2 * (v1 + l) / (v1 * (v2 - 2));
 111       } // mean
 112
 113       template <class RealType, class Policy>
 114       inline RealType mode(const non_central_f_distribution<RealType, Policy>& dist)
 115       { // mode.
 116          static const char* function = "mode(non_central_chi_squared_distribution<%1%> const&)";
 117
 118          RealType n = dist.degrees_of_freedom1();
 119          RealType m = dist.degrees_of_freedom2();
 120          RealType l = dist.non_centrality();
 121          RealType r;
 122          if(!detail::check_df(
 123             function,
 124             n, &r, Policy())
 125                ||
 126             !detail::check_df(
 127                function,
 128                m, &r, Policy())
 129                ||
 130             !detail::check_non_centrality(
 131                function,
 132                l,
 133                &r,
 134                Policy()))
 135                return r;
 136          RealType guess = m > 2 ? RealType(m * (n + l) / (n * (m - 2))) : RealType(1);
 137          return detail::generic_find_mode(
 138             dist,
 139             guess,
 140             function);
 141       }
 142
 143       template <class RealType, class Policy>
 144       inline RealType variance(const non_central_f_distribution<RealType, Policy>& dist)
 145       { // variance.
 146          const char* function = "variance(non_central_f_distribution<%1%> const&)";
 147          RealType n = dist.degrees_of_freedom1();
 148          RealType m = dist.degrees_of_freedom2();
 149          RealType l = dist.non_centrality();
 150          RealType r;
 151          if(!detail::check_df(
 152             function,
 153             n, &r, Policy())
 154                ||
 155             !detail::check_df(
 156                function,
 157                m, &r, Policy())
 158                ||
 159             !detail::check_non_centrality(
 160                function,
 161                l,
 162                &r,
 163                Policy()))
 164                return r;
 165          if(m <= 4)
 166             return policies::raise_domain_error(
 167                function,
 168                "Second degrees of freedom parameter was %1%, but must be > 4 !",
 169                m, Policy());
 170          RealType result = 2 * m * m * ((n + l) * (n + l)
 171             + (m - 2) * (n + 2 * l));
 172          result /= (m - 4) * (m - 2) * (m - 2) * n * n;
 173          return result;
 174       }
 175
 176       // RealType standard_deviation(const non_central_f_distribution<RealType, Policy>& dist)
 177       // standard_deviation provided by derived accessors.
 178
 179       template <class RealType, class Policy>
 180       inline RealType skewness(const non_central_f_distribution<RealType, Policy>& dist)
 181       { // skewness = sqrt(l).
 182          const char* function = "skewness(non_central_f_distribution<%1%> const&)";
 183          BOOST_MATH_STD_USING
 184          RealType n = dist.degrees_of_freedom1();
 185          RealType m = dist.degrees_of_freedom2();
 186          RealType l = dist.non_centrality();
 187          RealType r;
 188          if(!detail::check_df(
 189             function,
 190             n, &r, Policy())
 191                ||
 192             !detail::check_df(
 193                function,
 194                m, &r, Policy())
 195                ||
 196             !detail::check_non_centrality(
 197                function,
 198                l,
 199                &r,
 200                Policy()))
 201                return r;
 202          if(m <= 6)
 203             return policies::raise_domain_error(
 204                function,
 205                "Second degrees of freedom parameter was %1%, but must be > 6 !",
 206                m, Policy());
 207          RealType result = 2 * constants::root_two<RealType>();
 208          result *= sqrt(m - 4);
 209          result *= (n * (m + n - 2) *(m + 2 * n - 2)
 210             + 3 * (m + n - 2) * (m + 2 * n - 2) * l
 211             + 6 * (m + n - 2) * l * l + 2 * l * l * l);
 212          result /= (m - 6) * pow(n * (m + n - 2) + 2 * (m + n - 2) * l + l * l, RealType(1.5f));
 213          return result;
 214       }
 215
 216       template <class RealType, class Policy>
 217       inline RealType kurtosis_excess(const non_central_f_distribution<RealType, Policy>& dist)
 218       {
 219          const char* function = "kurtosis_excess(non_central_f_distribution<%1%> const&)";
 220          BOOST_MATH_STD_USING
 221          RealType n = dist.degrees_of_freedom1();
 222          RealType m = dist.degrees_of_freedom2();
 223          RealType l = dist.non_centrality();
 224          RealType r;
 225          if(!detail::check_df(
 226             function,
 227             n, &r, Policy())
 228                ||
 229             !detail::check_df(
 230                function,
 231                m, &r, Policy())
 232                ||
 233             !detail::check_non_centrality(
 234                function,
 235                l,
 236                &r,
 237                Policy()))
 238                return r;
 239          if(m <= 8)
 240             return policies::raise_domain_error(
 241                function,
 242                "Second degrees of freedom parameter was %1%, but must be > 8 !",
 243                m, Policy());
 244          RealType l2 = l * l;
 245          RealType l3 = l2 * l;
 246          RealType l4 = l2 * l2;
 247          RealType result = (3 * (m - 4) * (n * (m + n - 2)
 248             * (4 * (m - 2) * (m - 2)
 249             + (m - 2) * (m + 10) * n
 250             + (10 + m) * n * n)
 251             + 4 * (m + n - 2) * (4 * (m - 2) * (m - 2)
 252             + (m - 2) * (10 + m) * n
 253             + (10 + m) * n * n) * l + 2 * (10 + m)
 254             * (m + n - 2) * (2 * m + 3 * n - 4) * l2
 255             + 4 * (10 + m) * (-2 + m + n) * l3
 256             + (10 + m) * l4))
 257             /
 258             ((-8 + m) * (-6 + m) * boost::math::pow<2>(n * (-2 + m + n)
 259             + 2 * (-2 + m + n) * l + l2));
 260             return result;
 261       } // kurtosis_excess
 262
 263       template <class RealType, class Policy>
 264       inline RealType kurtosis(const non_central_f_distribution<RealType, Policy>& dist)
 265       {
 266          return kurtosis_excess(dist) + 3;
 267       }
 268
 269       template <class RealType, class Policy>
 270       inline RealType pdf(const non_central_f_distribution<RealType, Policy>& dist, const RealType& x)
 271       { // Probability Density/Mass Function.
 272          typedef typename policies::evaluation<RealType, Policy>::type value_type;
 273          typedef typename policies::normalise<
 274             Policy,
 275             policies::promote_float<false>,
 276             policies::promote_double<false>,
 277             policies::discrete_quantile<>,
 278             policies::assert_undefined<> >::type forwarding_policy;
 279
 280          value_type alpha = dist.degrees_of_freedom1() / 2;
 281          value_type beta = dist.degrees_of_freedom2() / 2;
 282          value_type y = x * alpha / beta;
 283          value_type r = pdf(boost::math::non_central_beta_distribution<value_type, forwarding_policy>(alpha, beta, dist.non_centrality()), y / (1 + y));
 284          return policies::checked_narrowing_cast<RealType, forwarding_policy>(
 285             r * (dist.degrees_of_freedom1() / dist.degrees_of_freedom2()) / ((1 + y) * (1 + y)),
 286             "pdf(non_central_f_distribution<%1%>, %1%)");
 287       } // pdf
 288
 289       template <class RealType, class Policy>
 290       RealType cdf(const non_central_f_distribution<RealType, Policy>& dist, const RealType& x)
 291       {
 292          const char* function = "cdf(const non_central_f_distribution<%1%>&, %1%)";
 293          RealType r;
 294          if(!detail::check_df(
 295             function,
 296             dist.degrees_of_freedom1(), &r, Policy())
 297                ||
 298             !detail::check_df(
 299                function,
 300                dist.degrees_of_freedom2(), &r, Policy())
 301                ||
 302             !detail::check_non_centrality(
 303                function,
 304                dist.non_centrality(),
 305                &r,
 306                Policy()))
 307                return r;
 308
 309          if((x < 0) || !(boost::math::isfinite)(x))
 310          {
 311             return policies::raise_domain_error<RealType>(
 312                function, "Random Variable parameter was %1%, but must be > 0 !", x, Policy());
 313          }
 314
 315          RealType alpha = dist.degrees_of_freedom1() / 2;
 316          RealType beta = dist.degrees_of_freedom2() / 2;
 317          RealType y = x * alpha / beta;
 318          RealType c = y / (1 + y);
 319          RealType cp = 1 / (1 + y);
 320          //
 321          // To ensure accuracy, we pass both x and 1-x to the
 322          // non-central beta cdf routine, this ensures accuracy
 323          // even when we compute x to be ~ 1:
 324          //
 325          r = detail::non_central_beta_cdf(c, cp, alpha, beta,
 326             dist.non_centrality(), false, Policy());
 327          return r;
 328       } // cdf
 329
 330       template <class RealType, class Policy>
 331       RealType cdf(const complemented2_type<non_central_f_distribution<RealType, Policy>, RealType>& c)
 332       { // Complemented Cumulative Distribution Function
 333          const char* function = "cdf(complement(const non_central_f_distribution<%1%>&, %1%))";
 334          RealType r;
 335          if(!detail::check_df(
 336             function,
 337             c.dist.degrees_of_freedom1(), &r, Policy())
 338                ||
 339             !detail::check_df(
 340                function,
 341                c.dist.degrees_of_freedom2(), &r, Policy())
 342                ||
 343             !detail::check_non_centrality(
 344                function,
 345                c.dist.non_centrality(),
 346                &r,
 347                Policy()))
 348                return r;
 349
 350          if((c.param < 0) || !(boost::math::isfinite)(c.param))
 351          {
 352             return policies::raise_domain_error<RealType>(
 353                function, "Random Variable parameter was %1%, but must be > 0 !", c.param, Policy());
 354          }
 355
 356          RealType alpha = c.dist.degrees_of_freedom1() / 2;
 357          RealType beta = c.dist.degrees_of_freedom2() / 2;
 358          RealType y = c.param * alpha / beta;
 359          RealType x = y / (1 + y);
 360          RealType cx = 1 / (1 + y);
 361          //
 362          // To ensure accuracy, we pass both x and 1-x to the
 363          // non-central beta cdf routine, this ensures accuracy
 364          // even when we compute x to be ~ 1:
 365          //
 366          r = detail::non_central_beta_cdf(x, cx, alpha, beta,
 367             c.dist.non_centrality(), true, Policy());
 368          return r;
 369       } // ccdf
 370
 371       template <class RealType, class Policy>
 372       inline RealType quantile(const non_central_f_distribution<RealType, Policy>& dist, const RealType& p)
 373       { // Quantile (or Percent Point) function.
 374          RealType alpha = dist.degrees_of_freedom1() / 2;
 375          RealType beta = dist.degrees_of_freedom2() / 2;
 376          RealType x = quantile(boost::math::non_central_beta_distribution<RealType, Policy>(alpha, beta, dist.non_centrality()), p);
 377          if(x == 1)
 378             return policies::raise_overflow_error<RealType>(
 379                "quantile(const non_central_f_distribution<%1%>&, %1%)",
 380                "Result of non central F quantile is too large to represent.",
 381                Policy());
 382          return (x / (1 - x)) * (dist.degrees_of_freedom2() / dist.degrees_of_freedom1());
 383       } // quantile
 384
 385       template <class RealType, class Policy>
 386       inline RealType quantile(const complemented2_type<non_central_f_distribution<RealType, Policy>, RealType>& c)
 387       { // Quantile (or Percent Point) function.
 388          RealType alpha = c.dist.degrees_of_freedom1() / 2;
 389          RealType beta = c.dist.degrees_of_freedom2() / 2;
 390          RealType x = quantile(complement(boost::math::non_central_beta_distribution<RealType, Policy>(alpha, beta, c.dist.non_centrality()), c.param));
 391          if(x == 1)
 392             return policies::raise_overflow_error<RealType>(
 393                "quantile(complement(const non_central_f_distribution<%1%>&, %1%))",
 394                "Result of non central F quantile is too large to represent.",
 395                Policy());
 396          return (x / (1 - x)) * (c.dist.degrees_of_freedom2() / c.dist.degrees_of_freedom1());
 397       } // quantile complement.
 398
 399    } // namespace math
 400 } // namespace boost
 401
 402 // This include must be at the end, *after* the accessors
 403 // for this distribution have been defined, in order to
 404 // keep compilers that support two-phase lookup happy.
 405 #include <boost/math/distributions/detail/derived_accessors.hpp>
 406
 407 #endif // BOOST_MATH_SPECIAL_NON_CENTRAL_F_HPP
 408
 409
 410