Chris@16: // boost\math\distributions\non_central_f.hpp
Chris@16: 
Chris@16: // Copyright John Maddock 2008.
Chris@16: 
Chris@16: // Use, modification and distribution are subject to the
Chris@16: // Boost Software License, Version 1.0.
Chris@16: // (See accompanying file LICENSE_1_0.txt
Chris@16: // or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@16: 
Chris@16: #ifndef BOOST_MATH_SPECIAL_NON_CENTRAL_F_HPP
Chris@16: #define BOOST_MATH_SPECIAL_NON_CENTRAL_F_HPP
Chris@16: 
Chris@16: #include <boost/math/distributions/non_central_beta.hpp>
Chris@16: #include <boost/math/distributions/detail/generic_mode.hpp>
Chris@16: #include <boost/math/special_functions/pow.hpp>
Chris@16: 
Chris@16: namespace boost
Chris@16: {
Chris@16:    namespace math
Chris@16:    {
Chris@16:       template <class RealType = double, class Policy = policies::policy<> >
Chris@16:       class non_central_f_distribution
Chris@16:       {
Chris@16:       public:
Chris@16:          typedef RealType value_type;
Chris@16:          typedef Policy policy_type;
Chris@16: 
Chris@16:          non_central_f_distribution(RealType v1_, RealType v2_, RealType lambda) : v1(v1_), v2(v2_), ncp(lambda)
Chris@16:          {
Chris@16:             const char* function = "boost::math::non_central_f_distribution<%1%>::non_central_f_distribution(%1%,%1%)";
Chris@16:             RealType r;
Chris@16:             detail::check_df(
Chris@16:                function,
Chris@16:                v1, &r, Policy());
Chris@16:             detail::check_df(
Chris@16:                function,
Chris@16:                v2, &r, Policy());
Chris@16:             detail::check_non_centrality(
Chris@16:                function,
Chris@16:                lambda,
Chris@16:                &r,
Chris@16:                Policy());
Chris@16:          } // non_central_f_distribution constructor.
Chris@16: 
Chris@16:          RealType degrees_of_freedom1()const
Chris@16:          {
Chris@16:             return v1;
Chris@16:          }
Chris@16:          RealType degrees_of_freedom2()const
Chris@16:          {
Chris@16:             return v2;
Chris@16:          }
Chris@16:          RealType non_centrality() const
Chris@16:          { // Private data getter function.
Chris@16:             return ncp;
Chris@16:          }
Chris@16:       private:
Chris@16:          // Data member, initialized by constructor.
Chris@16:          RealType v1;   // alpha.
Chris@16:          RealType v2;   // beta.
Chris@16:          RealType ncp; // non-centrality parameter
Chris@16:       }; // template <class RealType, class Policy> class non_central_f_distribution
Chris@16: 
Chris@16:       typedef non_central_f_distribution<double> non_central_f; // Reserved name of type double.
Chris@16: 
Chris@16:       // Non-member functions to give properties of the distribution.
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       inline const std::pair<RealType, RealType> range(const non_central_f_distribution<RealType, Policy>& /* dist */)
Chris@16:       { // Range of permissible values for random variable k.
Chris@16:          using boost::math::tools::max_value;
Chris@16:          return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());
Chris@16:       }
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       inline const std::pair<RealType, RealType> support(const non_central_f_distribution<RealType, Policy>& /* dist */)
Chris@16:       { // Range of supported values for random variable k.
Chris@16:          // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
Chris@16:          using boost::math::tools::max_value;
Chris@16:          return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());
Chris@16:       }
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       inline RealType mean(const non_central_f_distribution<RealType, Policy>& dist)
Chris@16:       {
Chris@16:          const char* function = "mean(non_central_f_distribution<%1%> const&)";
Chris@16:          RealType v1 = dist.degrees_of_freedom1();
Chris@16:          RealType v2 = dist.degrees_of_freedom2();
Chris@16:          RealType l = dist.non_centrality();
Chris@16:          RealType r;
Chris@16:          if(!detail::check_df(
Chris@16:             function,
Chris@16:             v1, &r, Policy())
Chris@16:                ||
Chris@16:             !detail::check_df(
Chris@16:                function,
Chris@16:                v2, &r, Policy())
Chris@16:                ||
Chris@16:             !detail::check_non_centrality(
Chris@16:                function,
Chris@16:                l,
Chris@16:                &r,
Chris@16:                Policy()))
Chris@16:                return r;
Chris@16:          if(v2 <= 2)
Chris@16:             return policies::raise_domain_error(
Chris@16:                function,
Chris@16:                "Second degrees of freedom parameter was %1%, but must be > 2 !",
Chris@16:                v2, Policy());
Chris@16:          return v2 * (v1 + l) / (v1 * (v2 - 2));
Chris@16:       } // mean
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       inline RealType mode(const non_central_f_distribution<RealType, Policy>& dist)
Chris@16:       { // mode.
Chris@16:          static const char* function = "mode(non_central_chi_squared_distribution<%1%> const&)";
Chris@16: 
Chris@16:          RealType n = dist.degrees_of_freedom1();
Chris@16:          RealType m = dist.degrees_of_freedom2();
Chris@16:          RealType l = dist.non_centrality();
Chris@16:          RealType r;
Chris@16:          if(!detail::check_df(
Chris@16:             function,
Chris@16:             n, &r, Policy())
Chris@16:                ||
Chris@16:             !detail::check_df(
Chris@16:                function,
Chris@16:                m, &r, Policy())
Chris@16:                ||
Chris@16:             !detail::check_non_centrality(
Chris@16:                function,
Chris@16:                l,
Chris@16:                &r,
Chris@16:                Policy()))
Chris@16:                return r;
Chris@101:          RealType guess = m > 2 ? RealType(m * (n + l) / (n * (m - 2))) : RealType(1);
Chris@16:          return detail::generic_find_mode(
Chris@16:             dist,
Chris@101:             guess,
Chris@16:             function);
Chris@16:       }
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       inline RealType variance(const non_central_f_distribution<RealType, Policy>& dist)
Chris@16:       { // variance.
Chris@16:          const char* function = "variance(non_central_f_distribution<%1%> const&)";
Chris@16:          RealType n = dist.degrees_of_freedom1();
Chris@16:          RealType m = dist.degrees_of_freedom2();
Chris@16:          RealType l = dist.non_centrality();
Chris@16:          RealType r;
Chris@16:          if(!detail::check_df(
Chris@16:             function,
Chris@16:             n, &r, Policy())
Chris@16:                ||
Chris@16:             !detail::check_df(
Chris@16:                function,
Chris@16:                m, &r, Policy())
Chris@16:                ||
Chris@16:             !detail::check_non_centrality(
Chris@16:                function,
Chris@16:                l,
Chris@16:                &r,
Chris@16:                Policy()))
Chris@16:                return r;
Chris@16:          if(m <= 4)
Chris@16:             return policies::raise_domain_error(
Chris@16:                function,
Chris@16:                "Second degrees of freedom parameter was %1%, but must be > 4 !",
Chris@16:                m, Policy());
Chris@16:          RealType result = 2 * m * m * ((n + l) * (n + l)
Chris@16:             + (m - 2) * (n + 2 * l));
Chris@16:          result /= (m - 4) * (m - 2) * (m - 2) * n * n;
Chris@16:          return result;
Chris@16:       }
Chris@16: 
Chris@16:       // RealType standard_deviation(const non_central_f_distribution<RealType, Policy>& dist)
Chris@16:       // standard_deviation provided by derived accessors.
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       inline RealType skewness(const non_central_f_distribution<RealType, Policy>& dist)
Chris@16:       { // skewness = sqrt(l).
Chris@16:          const char* function = "skewness(non_central_f_distribution<%1%> const&)";
Chris@16:          BOOST_MATH_STD_USING
Chris@16:          RealType n = dist.degrees_of_freedom1();
Chris@16:          RealType m = dist.degrees_of_freedom2();
Chris@16:          RealType l = dist.non_centrality();
Chris@16:          RealType r;
Chris@16:          if(!detail::check_df(
Chris@16:             function,
Chris@16:             n, &r, Policy())
Chris@16:                ||
Chris@16:             !detail::check_df(
Chris@16:                function,
Chris@16:                m, &r, Policy())
Chris@16:                ||
Chris@16:             !detail::check_non_centrality(
Chris@16:                function,
Chris@16:                l,
Chris@16:                &r,
Chris@16:                Policy()))
Chris@16:                return r;
Chris@16:          if(m <= 6)
Chris@16:             return policies::raise_domain_error(
Chris@16:                function,
Chris@16:                "Second degrees of freedom parameter was %1%, but must be > 6 !",
Chris@16:                m, Policy());
Chris@16:          RealType result = 2 * constants::root_two<RealType>();
Chris@16:          result *= sqrt(m - 4);
Chris@16:          result *= (n * (m + n - 2) *(m + 2 * n - 2)
Chris@16:             + 3 * (m + n - 2) * (m + 2 * n - 2) * l
Chris@16:             + 6 * (m + n - 2) * l * l + 2 * l * l * l);
Chris@16:          result /= (m - 6) * pow(n * (m + n - 2) + 2 * (m + n - 2) * l + l * l, RealType(1.5f));
Chris@16:          return result;
Chris@16:       }
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       inline RealType kurtosis_excess(const non_central_f_distribution<RealType, Policy>& dist)
Chris@16:       {
Chris@16:          const char* function = "kurtosis_excess(non_central_f_distribution<%1%> const&)";
Chris@16:          BOOST_MATH_STD_USING
Chris@16:          RealType n = dist.degrees_of_freedom1();
Chris@16:          RealType m = dist.degrees_of_freedom2();
Chris@16:          RealType l = dist.non_centrality();
Chris@16:          RealType r;
Chris@16:          if(!detail::check_df(
Chris@16:             function,
Chris@16:             n, &r, Policy())
Chris@16:                ||
Chris@16:             !detail::check_df(
Chris@16:                function,
Chris@16:                m, &r, Policy())
Chris@16:                ||
Chris@16:             !detail::check_non_centrality(
Chris@16:                function,
Chris@16:                l,
Chris@16:                &r,
Chris@16:                Policy()))
Chris@16:                return r;
Chris@16:          if(m <= 8)
Chris@16:             return policies::raise_domain_error(
Chris@16:                function,
Chris@16:                "Second degrees of freedom parameter was %1%, but must be > 8 !",
Chris@16:                m, Policy());
Chris@16:          RealType l2 = l * l;
Chris@16:          RealType l3 = l2 * l;
Chris@16:          RealType l4 = l2 * l2;
Chris@16:          RealType result = (3 * (m - 4) * (n * (m + n - 2)
Chris@16:             * (4 * (m - 2) * (m - 2)
Chris@16:             + (m - 2) * (m + 10) * n
Chris@16:             + (10 + m) * n * n)
Chris@16:             + 4 * (m + n - 2) * (4 * (m - 2) * (m - 2)
Chris@16:             + (m - 2) * (10 + m) * n
Chris@16:             + (10 + m) * n * n) * l + 2 * (10 + m)
Chris@16:             * (m + n - 2) * (2 * m + 3 * n - 4) * l2
Chris@16:             + 4 * (10 + m) * (-2 + m + n) * l3
Chris@16:             + (10 + m) * l4))
Chris@16:             /
Chris@16:             ((-8 + m) * (-6 + m) * boost::math::pow<2>(n * (-2 + m + n)
Chris@16:             + 2 * (-2 + m + n) * l + l2));
Chris@16:             return result;
Chris@16:       } // kurtosis_excess
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       inline RealType kurtosis(const non_central_f_distribution<RealType, Policy>& dist)
Chris@16:       {
Chris@16:          return kurtosis_excess(dist) + 3;
Chris@16:       }
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       inline RealType pdf(const non_central_f_distribution<RealType, Policy>& dist, const RealType& x)
Chris@16:       { // Probability Density/Mass Function.
Chris@16:          typedef typename policies::evaluation<RealType, Policy>::type value_type;
Chris@16:          typedef typename policies::normalise<
Chris@16:             Policy,
Chris@16:             policies::promote_float<false>,
Chris@16:             policies::promote_double<false>,
Chris@16:             policies::discrete_quantile<>,
Chris@16:             policies::assert_undefined<> >::type forwarding_policy;
Chris@16: 
Chris@16:          value_type alpha = dist.degrees_of_freedom1() / 2;
Chris@16:          value_type beta = dist.degrees_of_freedom2() / 2;
Chris@16:          value_type y = x * alpha / beta;
Chris@16:          value_type r = pdf(boost::math::non_central_beta_distribution<value_type, forwarding_policy>(alpha, beta, dist.non_centrality()), y / (1 + y));
Chris@16:          return policies::checked_narrowing_cast<RealType, forwarding_policy>(
Chris@16:             r * (dist.degrees_of_freedom1() / dist.degrees_of_freedom2()) / ((1 + y) * (1 + y)),
Chris@16:             "pdf(non_central_f_distribution<%1%>, %1%)");
Chris@16:       } // pdf
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       RealType cdf(const non_central_f_distribution<RealType, Policy>& dist, const RealType& x)
Chris@16:       {
Chris@16:          const char* function = "cdf(const non_central_f_distribution<%1%>&, %1%)";
Chris@16:          RealType r;
Chris@16:          if(!detail::check_df(
Chris@16:             function,
Chris@16:             dist.degrees_of_freedom1(), &r, Policy())
Chris@16:                ||
Chris@16:             !detail::check_df(
Chris@16:                function,
Chris@16:                dist.degrees_of_freedom2(), &r, Policy())
Chris@16:                ||
Chris@16:             !detail::check_non_centrality(
Chris@16:                function,
Chris@16:                dist.non_centrality(),
Chris@16:                &r,
Chris@16:                Policy()))
Chris@16:                return r;
Chris@16: 
Chris@16:          if((x < 0) || !(boost::math::isfinite)(x))
Chris@16:          {
Chris@16:             return policies::raise_domain_error<RealType>(
Chris@16:                function, "Random Variable parameter was %1%, but must be > 0 !", x, Policy());
Chris@16:          }
Chris@16: 
Chris@16:          RealType alpha = dist.degrees_of_freedom1() / 2;
Chris@16:          RealType beta = dist.degrees_of_freedom2() / 2;
Chris@16:          RealType y = x * alpha / beta;
Chris@16:          RealType c = y / (1 + y);
Chris@16:          RealType cp = 1 / (1 + y);
Chris@16:          //
Chris@16:          // To ensure accuracy, we pass both x and 1-x to the
Chris@16:          // non-central beta cdf routine, this ensures accuracy
Chris@16:          // even when we compute x to be ~ 1:
Chris@16:          //
Chris@16:          r = detail::non_central_beta_cdf(c, cp, alpha, beta,
Chris@16:             dist.non_centrality(), false, Policy());
Chris@16:          return r;
Chris@16:       } // cdf
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       RealType cdf(const complemented2_type<non_central_f_distribution<RealType, Policy>, RealType>& c)
Chris@16:       { // Complemented Cumulative Distribution Function
Chris@16:          const char* function = "cdf(complement(const non_central_f_distribution<%1%>&, %1%))";
Chris@16:          RealType r;
Chris@16:          if(!detail::check_df(
Chris@16:             function,
Chris@16:             c.dist.degrees_of_freedom1(), &r, Policy())
Chris@16:                ||
Chris@16:             !detail::check_df(
Chris@16:                function,
Chris@16:                c.dist.degrees_of_freedom2(), &r, Policy())
Chris@16:                ||
Chris@16:             !detail::check_non_centrality(
Chris@16:                function,
Chris@16:                c.dist.non_centrality(),
Chris@16:                &r,
Chris@16:                Policy()))
Chris@16:                return r;
Chris@16: 
Chris@16:          if((c.param < 0) || !(boost::math::isfinite)(c.param))
Chris@16:          {
Chris@16:             return policies::raise_domain_error<RealType>(
Chris@16:                function, "Random Variable parameter was %1%, but must be > 0 !", c.param, Policy());
Chris@16:          }
Chris@16: 
Chris@16:          RealType alpha = c.dist.degrees_of_freedom1() / 2;
Chris@16:          RealType beta = c.dist.degrees_of_freedom2() / 2;
Chris@16:          RealType y = c.param * alpha / beta;
Chris@16:          RealType x = y / (1 + y);
Chris@16:          RealType cx = 1 / (1 + y);
Chris@16:          //
Chris@16:          // To ensure accuracy, we pass both x and 1-x to the
Chris@16:          // non-central beta cdf routine, this ensures accuracy
Chris@16:          // even when we compute x to be ~ 1:
Chris@16:          //
Chris@16:          r = detail::non_central_beta_cdf(x, cx, alpha, beta,
Chris@16:             c.dist.non_centrality(), true, Policy());
Chris@16:          return r;
Chris@16:       } // ccdf
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       inline RealType quantile(const non_central_f_distribution<RealType, Policy>& dist, const RealType& p)
Chris@16:       { // Quantile (or Percent Point) function.
Chris@16:          RealType alpha = dist.degrees_of_freedom1() / 2;
Chris@16:          RealType beta = dist.degrees_of_freedom2() / 2;
Chris@16:          RealType x = quantile(boost::math::non_central_beta_distribution<RealType, Policy>(alpha, beta, dist.non_centrality()), p);
Chris@16:          if(x == 1)
Chris@16:             return policies::raise_overflow_error<RealType>(
Chris@16:                "quantile(const non_central_f_distribution<%1%>&, %1%)",
Chris@16:                "Result of non central F quantile is too large to represent.",
Chris@16:                Policy());
Chris@16:          return (x / (1 - x)) * (dist.degrees_of_freedom2() / dist.degrees_of_freedom1());
Chris@16:       } // quantile
Chris@16: 
Chris@16:       template <class RealType, class Policy>
Chris@16:       inline RealType quantile(const complemented2_type<non_central_f_distribution<RealType, Policy>, RealType>& c)
Chris@16:       { // Quantile (or Percent Point) function.
Chris@16:          RealType alpha = c.dist.degrees_of_freedom1() / 2;
Chris@16:          RealType beta = c.dist.degrees_of_freedom2() / 2;
Chris@16:          RealType x = quantile(complement(boost::math::non_central_beta_distribution<RealType, Policy>(alpha, beta, c.dist.non_centrality()), c.param));
Chris@16:          if(x == 1)
Chris@16:             return policies::raise_overflow_error<RealType>(
Chris@16:                "quantile(complement(const non_central_f_distribution<%1%>&, %1%))",
Chris@16:                "Result of non central F quantile is too large to represent.",
Chris@16:                Policy());
Chris@16:          return (x / (1 - x)) * (c.dist.degrees_of_freedom2() / c.dist.degrees_of_freedom1());
Chris@16:       } // quantile complement.
Chris@16: 
Chris@16:    } // namespace math
Chris@16: } // namespace boost
Chris@16: 
Chris@16: // This include must be at the end, *after* the accessors
Chris@16: // for this distribution have been defined, in order to
Chris@16: // keep compilers that support two-phase lookup happy.
Chris@16: #include <boost/math/distributions/detail/derived_accessors.hpp>
Chris@16: 
Chris@16: #endif // BOOST_MATH_SPECIAL_NON_CENTRAL_F_HPP
Chris@16: 
Chris@16: 
Chris@16: