binomial.hpp

来自「Boost provides free peer-reviewed portab」· HPP 代码 · 共 725 行 · 第 1/2 页

      //
      static RealType find_minimum_number_of_trials(
         RealType k,     // number of events
         RealType p,     // success fraction
         RealType alpha) // risk level
      {
        static const char* function = "boost::math::binomial_distribution<%1%>::find_minimum_number_of_trials";
        // Error checks:
        RealType result;
        if(false == binomial_detail::check_dist_and_k(
           function, k, p, k, &result, Policy())
            &&
           binomial_detail::check_dist_and_prob(
           function, k, p, alpha, &result, Policy()))
        { return result; }

        result = ibetac_invb(k + 1, p, alpha, Policy());  // returns n - k
        return result + k;
      }

      static RealType find_maximum_number_of_trials(
         RealType k,     // number of events
         RealType p,     // success fraction
         RealType alpha) // risk level
      {
        static const char* function = "boost::math::binomial_distribution<%1%>::find_maximum_number_of_trials";
        // Error checks:
        RealType result;
        if(false == binomial_detail::check_dist_and_k(
           function, k, p, k, &result, Policy())
            &&
           binomial_detail::check_dist_and_prob(
           function, k, p, alpha, &result, Policy()))
        { return result; }

        result = ibeta_invb(k + 1, p, alpha, Policy());  // returns n - k
        return result + k;
      }

    private:
        RealType m_n; // Not sure if this shouldn't be an int?
        RealType m_p; // success_fraction
      }; // template <class RealType, class Policy> class binomial_distribution

      typedef binomial_distribution<> binomial;
      // typedef binomial_distribution<double> binomial;
      // IS now included since no longer a name clash with function binomial.
      //typedef binomial_distribution<double> binomial; // Reserved name of type double.

      template <class RealType, class Policy>
      const std::pair<RealType, RealType> range(const binomial_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), dist.trials());
      }

      template <class RealType, class Policy>
      const std::pair<RealType, RealType> support(const binomial_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.
        return std::pair<RealType, RealType>(0,  dist.trials());
      }

      template <class RealType, class Policy>
      inline RealType mean(const binomial_distribution<RealType, Policy>& dist)
      { // Mean of Binomial distribution = np.
        return  dist.trials() * dist.success_fraction();
      } // mean

      template <class RealType, class Policy>
      inline RealType variance(const binomial_distribution<RealType, Policy>& dist)
      { // Variance of Binomial distribution = np(1-p).
        return  dist.trials() * dist.success_fraction() * (1 - dist.success_fraction());
      } // variance

      template <class RealType, class Policy>
      RealType pdf(const binomial_distribution<RealType, Policy>& dist, const RealType& k)
      { // Probability Density/Mass Function.
        BOOST_FPU_EXCEPTION_GUARD

        BOOST_MATH_STD_USING // for ADL of std functions

        RealType n = dist.trials();

        // Error check:
        RealType result;
        if(false == binomial_detail::check_dist_and_k(
           "boost::math::pdf(binomial_distribution<%1%> const&, %1%)",
           n,
           dist.success_fraction(),
           k,
           &result, Policy()))
        {
           return result;
        }

        // Special cases of success_fraction, regardless of k successes and regardless of n trials.
        if (dist.success_fraction() == 0)
        {  // probability of zero successes is 1:
           return static_cast<RealType>(k == 0 ? 1 : 0);
        }
        if (dist.success_fraction() == 1)
        {  // probability of n successes is 1:
           return static_cast<RealType>(k == n ? 1 : 0);
        }
        // k argument may be integral, signed, or unsigned, or floating point.
        // If necessary, it has already been promoted from an integral type.
        if (n == 0)
        {
          return 1; // Probability = 1 = certainty.
        }
        if (k == 0)
        { // binomial coeffic (n 0) = 1,
          // n ^ 0 = 1
          return pow(1 - dist.success_fraction(), n);
        }
        if (k == n)
        { // binomial coeffic (n n) = 1,
          // n ^ 0 = 1
          return pow(dist.success_fraction(), k);  // * pow((1 - dist.success_fraction()), (n - k)) = 1
        }

        // Probability of getting exactly k successes
        // if C(n, k) is the binomial coefficient then:
        //
        // f(k; n,p) = C(n, k) * p^k * (1-p)^(n-k)
        //           = (n!/(k!(n-k)!)) * p^k * (1-p)^(n-k)
        //           = (tgamma(n+1) / (tgamma(k+1)*tgamma(n-k+1))) * p^k * (1-p)^(n-k)
        //           = p^k (1-p)^(n-k) / (beta(k+1, n-k+1) * (n+1))
        //           = ibeta_derivative(k+1, n-k+1, p) / (n+1)
        //
        using boost::math::ibeta_derivative; // a, b, x
        return ibeta_derivative(k+1, n-k+1, dist.success_fraction(), Policy()) / (n+1);

      } // pdf

      template <class RealType, class Policy>
      inline RealType cdf(const binomial_distribution<RealType, Policy>& dist, const RealType& k)
      { // Cumulative Distribution Function Binomial.
        // The random variate k is the number of successes in n trials.
        // k argument may be integral, signed, or unsigned, or floating point.
        // If necessary, it has already been promoted from an integral type.

        // Returns the sum of the terms 0 through k of the Binomial Probability Density/Mass:
        //
        //   i=k
        //   --  ( n )   i      n-i
        //   >   |   |  p  (1-p)
        //   --  ( i )
        //   i=0

        // The terms are not summed directly instead
        // the incomplete beta integral is employed,
        // according to the formula:
        // P = I[1-p]( n-k, k+1).
        //   = 1 - I[p](k + 1, n - k)

        BOOST_MATH_STD_USING // for ADL of std functions

        RealType n = dist.trials();
        RealType p = dist.success_fraction();

        // Error check:
        RealType result;
        if(false == binomial_detail::check_dist_and_k(
           "boost::math::cdf(binomial_distribution<%1%> const&, %1%)",
           n,
           p,
           k,
           &result, Policy()))
        {
           return result;
        }
        if (k == n)
        {
          return 1;
        }

        // Special cases, regardless of k.
        if (p == 0)
        {  // This need explanation:
           // the pdf is zero for all cases except when k == 0.
           // For zero p the probability of zero successes is one.
           // Therefore the cdf is always 1:
           // the probability of k or *fewer* successes is always 1
           // if there are never any successes!
           return 1;
        }
        if (p == 1)
        { // This is correct but needs explanation:
          // when k = 1
          // all the cdf and pdf values are zero *except* when k == n,
          // and that case has been handled above already.
          return 0;
        }
        //
        // P = I[1-p](n - k, k + 1)
        //   = 1 - I[p](k + 1, n - k)
        // Use of ibetac here prevents cancellation errors in calculating
        // 1-p if p is very small, perhaps smaller than machine epsilon.
        //
        // Note that we do not use a finite sum here, since the incomplete
        // beta uses a finite sum internally for integer arguments, so
        // we'll just let it take care of the necessary logic.
        //
        return ibetac(k + 1, n - k, p, Policy());
      } // binomial cdf

      template <class RealType, class Policy>
      inline RealType cdf(const complemented2_type<binomial_distribution<RealType, Policy>, RealType>& c)
      { // Complemented Cumulative Distribution Function Binomial.
        // The random variate k is the number of successes in n trials.
        // k argument may be integral, signed, or unsigned, or floating point.
        // If necessary, it has already been promoted from an integral type.

        // Returns the sum of the terms k+1 through n of the Binomial Probability Density/Mass:
        //
        //   i=n
        //   --  ( n )   i      n-i
        //   >   |   |  p  (1-p)
        //   --  ( i )
        //   i=k+1

        // The terms are not summed directly instead
        // the incomplete beta integral is employed,
        // according to the formula:
        // Q = 1 -I[1-p]( n-k, k+1).
        //   = I[p](k + 1, n - k)

        BOOST_MATH_STD_USING // for ADL of std functions

        RealType const& k = c.param;
        binomial_distribution<RealType, Policy> const& dist = c.dist;
        RealType n = dist.trials();
        RealType p = dist.success_fraction();

        // Error checks:
        RealType result;
        if(false == binomial_detail::check_dist_and_k(
           "boost::math::cdf(binomial_distribution<%1%> const&, %1%)",
           n,
           p,
           k,
           &result, Policy()))
        {
           return result;
        }

        if (k == n)
        { // Probability of greater than n successes is necessarily zero:
          return 0;
        }

        // Special cases, regardless of k.
        if (p == 0)
        {
           // This need explanation: the pdf is zero for all
           // cases except when k == 0.  For zero p the probability
           // of zero successes is one.  Therefore the cdf is always
           // 1: the probability of *more than* k successes is always 0
           // if there are never any successes!
           return 0;
        }
        if (p == 1)
        {
          // This needs explanation, when p = 1
          // we always have n successes, so the probability
          // of more than k successes is 1 as long as k < n.
          // The k == n case has already been handled above.
          return 1;
        }
        //
        // Calculate cdf binomial using the incomplete beta function.
        // Q = 1 -I[1-p](n - k, k + 1)
        //   = I[p](k + 1, n - k)
        // Use of ibeta here prevents cancellation errors in calculating
        // 1-p if p is very small, perhaps smaller than machine epsilon.
        //
        // Note that we do not use a finite sum here, since the incomplete
        // beta uses a finite sum internally for integer arguments, so
        // we'll just let it take care of the necessary logic.
        //
        return ibeta(k + 1, n - k, p, Policy());
      } // binomial cdf

      template <class RealType, class Policy>
      inline RealType quantile(const binomial_distribution<RealType, Policy>& dist, const RealType& p)
      {
         return binomial_detail::quantile_imp(dist, p, 1-p);
      } // quantile

      template <class RealType, class Policy>
      RealType quantile(const complemented2_type<binomial_distribution<RealType, Policy>, RealType>& c)
      {
         return binomial_detail::quantile_imp(c.dist, 1-c.param, c.param);
      } // quantile

      template <class RealType, class Policy>
      inline RealType mode(const binomial_distribution<RealType, Policy>& dist)
      {
         BOOST_MATH_STD_USING // ADL of std functions.
         RealType p = dist.success_fraction();
         RealType n = dist.trials();
         return floor(p * (n + 1));
      }

      template <class RealType, class Policy>
      inline RealType median(const binomial_distribution<RealType, Policy>& dist)
      { // Bounds for the median of the negative binomial distribution
        // VAN DE VEN R. ; WEBER N. C. ;
        // Univ. Sydney, school mathematics statistics, Sydney N.S.W. 2006, AUSTRALIE
        // Metrika  (Metrika)  ISSN 0026-1335   CODEN MTRKA8
        // 1993, vol. 40, no3-4, pp. 185-189 (4 ref.)

        // Bounds for median and 50 percetage point of binomial and negative binomial distribution
        // Metrika, ISSN   0026-1335 (Print) 1435-926X (Online)
        // Volume 41, Number 1 / December, 1994, DOI   10.1007/BF01895303
         BOOST_MATH_STD_USING // ADL of std functions.
         RealType p = dist.success_fraction();
         RealType n = dist.trials();
         // Wikipedia says one of floor(np) -1, floor (np), floor(np) +1
         return floor(p * n); // Chose the middle value.
      }

      template <class RealType, class Policy>
      inline RealType skewness(const binomial_distribution<RealType, Policy>& dist)
      {
         BOOST_MATH_STD_USING // ADL of std functions.
         RealType p = dist.success_fraction();
         RealType n = dist.trials();
         return (1 - 2 * p) / sqrt(n * p * (1 - p));
      }

      template <class RealType, class Policy>
      inline RealType kurtosis(const binomial_distribution<RealType, Policy>& dist)
      {
         RealType p = dist.success_fraction();
         RealType n = dist.trials();
         return 3 - 6 / n + 1 / (n * p * (1 - p));
      }

      template <class RealType, class Policy>
      inline RealType kurtosis_excess(const binomial_distribution<RealType, Policy>& dist)
      {
         RealType p = dist.success_fraction();
         RealType q = 1 - p;
         RealType n = dist.trials();
         return (1 - 6 * p * q) / (n * p * q);
      }

    } // 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_BINOMIAL_HPP