SV platform dependency builds

// boost\math\distributions\binomial.hpp

// Copyright John Maddock 2006.

// Copyright Paul A. Bristow 2007.

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

// http://en.wikipedia.org/wiki/binomial_distribution

// Binomial distribution is the discrete probability distribution of

// the number (k) of successes, in a sequence of

// n independent (yes or no, success or failure) Bernoulli trials.

// It expresses the probability of a number of events occurring in a fixed time

// if these events occur with a known average rate (probability of success),

// and are independent of the time since the last event.

// The number of cars that pass through a certain point on a road during a given period of time.

// The number of spelling mistakes a secretary makes while typing a single page.

// The number of phone calls at a call center per minute.

// The number of times a web server is accessed per minute.

// The number of light bulbs that burn out in a certain amount of time.

// The number of roadkill found per unit length of road

// http://en.wikipedia.org/wiki/binomial_distribution

// Given a sample of N measured values k[i],

// we wish to estimate the value of the parameter x (mean)

// of the binomial population from which the sample was drawn.

// To calculate the maximum likelihood value = 1/N sum i = 1 to N of k[i]

// Also may want a function for EXACTLY k.

// And probability that there are EXACTLY k occurrences is

// exp(-x) * pow(x, k) / factorial(k)

// where x is expected occurrences (mean) during the given interval.

// For example, if events occur, on average, every 4 min,

// and we are interested in number of events occurring in 10 min,

// then x = 10/4 = 2.5

// http://www.itl.nist.gov/div898/handbook/eda/section3/eda366i.htm

// The binomial distribution is used when there are

// exactly two mutually exclusive outcomes of a trial.

// These outcomes are appropriately labeled "success" and "failure".

// The binomial distribution is used to obtain

// the probability of observing x successes in N trials,

// with the probability of success on a single trial denoted by p.

// The binomial distribution assumes that p is fixed for all trials.

// P(x, p, n) = n!/(x! * (n-x)!) * p^x * (1-p)^(n-x)

// http://mathworld.wolfram.com/BinomialCoefficient.html

// The binomial coefficient (n; k) is the number of ways of picking

// k unordered outcomes from n possibilities,

// also known as a combination or combinatorial number.

// The symbols _nC_k and (n; k) are used to denote a binomial coefficient,

// and are sometimes read as "n choose k."

// (n; k) therefore gives the number of k-subsets  possible out of a set of n distinct items.

// For example:

//  The 2-subsets of {1,2,3,4} are the six pairs {1,2}, {1,3}, {1,4}, {2,3}, {2,4}, and {3,4}, so (4; 2)==6.

// http://functions.wolfram.com/GammaBetaErf/Binomial/ for evaluation.

// But note that the binomial distribution

// (like others including the poisson, negative binomial & Bernoulli)

// is strictly defined as a discrete function: only integral values of k are envisaged.

// However because of the method of calculation using a continuous gamma function,

// it is convenient to treat it as if a continous function,

// and permit non-integral values of k.

// To enforce the strict mathematical model, users should use floor or ceil functions

// on k outside this function to ensure that k is integral.

#ifndef BOOST_MATH_SPECIAL_BINOMIAL_HPP

#define BOOST_MATH_SPECIAL_BINOMIAL_HPP

#include <boost/math/distributions/fwd.hpp>

#include <boost/math/special_functions/beta.hpp> // for incomplete beta.

#include <boost/math/distributions/complement.hpp> // complements

#include <boost/math/distributions/detail/common_error_handling.hpp> // error checks

#include <boost/math/distributions/detail/inv_discrete_quantile.hpp> // error checks

#include <boost/math/special_functions/fpclassify.hpp> // isnan.

#include <boost/math/tools/roots.hpp> // for root finding.

#include <utility>

namespace boost

{

  namespace math

  {

     template <class RealType, class Policy>

     class binomial_distribution;

     namespace binomial_detail{

        // common error checking routines for binomial distribution functions:

        template <class RealType, class Policy>

        inline bool check_N(const char* function, const RealType& N, RealType* result, const Policy& pol)

        {

           if((N < 0) || !(boost::math::isfinite)(N))

           {

               *result = policies::raise_domain_error<RealType>(

                  function,

                  "Number of Trials argument is %1%, but must be >= 0 !", N, pol);

               return false;

           }

           return true;

        }

        template <class RealType, class Policy>

        inline bool check_success_fraction(const char* function, const RealType& p, RealType* result, const Policy& pol)

        {

           if((p < 0) || (p > 1) || !(boost::math::isfinite)(p))

           {

               *result = policies::raise_domain_error<RealType>(

                  function,

                  "Success fraction argument is %1%, but must be >= 0 and <= 1 !", p, pol);

               return false;

           }

           return true;

        }

        template <class RealType, class Policy>

        inline bool check_dist(const char* function, const RealType& N, const RealType& p, RealType* result, const Policy& pol)

        {

           return check_success_fraction(

              function, p, result, pol)

              && check_N(

               function, N, result, pol);

        }

        template <class RealType, class Policy>

        inline bool check_dist_and_k(const char* function, const RealType& N, const RealType& p, RealType k, RealType* result, const Policy& pol)

        {

           if(check_dist(function, N, p, result, pol) == false)

              return false;

           if((k < 0) || !(boost::math::isfinite)(k))

           {

               *result = policies::raise_domain_error<RealType>(

                  function,

                  "Number of Successes argument is %1%, but must be >= 0 !", k, pol);

               return false;

           }

           if(k > N)

           {

               *result = policies::raise_domain_error<RealType>(

                  function,

                  "Number of Successes argument is %1%, but must be <= Number of Trials !", k, pol);

               return false;

           }

           return true;

        }

        template <class RealType, class Policy>

        inline bool check_dist_and_prob(const char* function, const RealType& N, RealType p, RealType prob, RealType* result, const Policy& pol)

        {

           if((check_dist(function, N, p, result, pol) && detail::check_probability(function, prob, result, pol)) == false)

              return false;

           return true;

        }

         template <class T, class Policy>

         T inverse_binomial_cornish_fisher(T n, T sf, T p, T q, const Policy& pol)

         {

            BOOST_MATH_STD_USING

            // mean:

            T m = n * sf;

            // standard deviation:

            T sigma = sqrt(n * sf * (1 - sf));

            // skewness

            T sk = (1 - 2 * sf) / sigma;

            // kurtosis:

            // T k = (1 - 6 * sf * (1 - sf) ) / (n * sf * (1 - sf));

            // Get the inverse of a std normal distribution:

            T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p, pol) * constants::root_two<T>();

            // Set the sign:

            if(p < 0.5)

               x = -x;

            T x2 = x * x;

            // w is correction term due to skewness

            T w = x + sk * (x2 - 1) / 6;

            /*

            // Add on correction due to kurtosis.

            // Disabled for now, seems to make things worse?

            //

            if(n >= 10)

               w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36;

               */

            w = m + sigma * w;

            if(w < tools::min_value<T>())

               return sqrt(tools::min_value<T>());

            if(w > n)

               return n;

            return w;

         }

      template <class RealType, class Policy>

      RealType quantile_imp(const binomial_distribution<RealType, Policy>& dist, const RealType& p, const RealType& q, bool comp)

      { // Quantile or Percent Point Binomial function.

        // Return the number of expected successes k,

        // for a given probability p.

        //

        // Error checks:

        BOOST_MATH_STD_USING  // ADL of std names

        RealType result = 0;

        RealType trials = dist.trials();

        RealType success_fraction = dist.success_fraction();

        if(false == binomial_detail::check_dist_and_prob(

           "boost::math::quantile(binomial_distribution<%1%> const&, %1%)",

           trials,

           success_fraction,

           p,

           &result, Policy()))

        {

           return result;

        }

        // Special cases:

        //

        if(p == 0)

        {  // There may actually be no answer to this question,

           // since the probability of zero successes may be non-zero,

           // but zero is the best we can do:

           return 0;

        }

        if(p == 1)

        {  // Probability of n or fewer successes is always one,

           // so n is the most sensible answer here:

           return trials;

        }

        if (p <= pow(1 - success_fraction, trials))

        { // p <= pdf(dist, 0) == cdf(dist, 0)

          return 0; // So the only reasonable result is zero.

        } // And root finder would fail otherwise.

        if(success_fraction == 1)

        {  // our formulae break down in this case:

           return p > 0.5f ? trials : 0;

        }

        // Solve for quantile numerically:

        //

        RealType guess = binomial_detail::inverse_binomial_cornish_fisher(trials, success_fraction, p, q, Policy());

        RealType factor = 8;

        if(trials > 100)

           factor = 1.01f; // guess is pretty accurate

        else if((trials > 10) && (trials - 1 > guess) && (guess > 3))

           factor = 1.15f; // less accurate but OK.

        else if(trials < 10)

        {

           // pretty inaccurate guess in this area:

           if(guess > trials / 64)

           {

              guess = trials / 4;

              factor = 2;

           }

           else

              guess = trials / 1024;

        }

        else

           factor = 2; // trials largish, but in far tails.

        typedef typename Policy::discrete_quantile_type discrete_quantile_type;

        boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();

        return detail::inverse_discrete_quantile(

            dist,

            comp ? q : p,

            comp,

            guess,

            factor,

            RealType(1),

            discrete_quantile_type(),

            max_iter);

      } // quantile

     }

    template <class RealType = double, class Policy = policies::policy<> >

    class binomial_distribution

    {

    public:

      typedef RealType value_type;

      typedef Policy policy_type;

      binomial_distribution(RealType n = 1, RealType p = 0.5) : m_n(n), m_p(p)

      { // Default n = 1 is the Bernoulli distribution

        // with equal probability of 'heads' or 'tails.

         RealType r;

         binomial_detail::check_dist(

            "boost::math::binomial_distribution<%1%>::binomial_distribution",

            m_n,

            m_p,

            &r, Policy());

      } // binomial_distribution constructor.

      RealType success_fraction() const

      { // Probability.

        return m_p;

      }

      RealType trials() const

      { // Total number of trials.

        return m_n;

      }

      enum interval_type{

         clopper_pearson_exact_interval,

         jeffreys_prior_interval

      };

      //

      // Estimation of the success fraction parameter.

      // The best estimate is actually simply successes/trials,

      // these functions are used

      // to obtain confidence intervals for the success fraction.

      //

      static RealType find_lower_bound_on_p(

         RealType trials,

         RealType successes,

         RealType probability,

         interval_type t = clopper_pearson_exact_interval)

      {

        static const char* function = "boost::math::binomial_distribution<%1%>::find_lower_bound_on_p";

        // Error checks:

        RealType result = 0;

        if(false == binomial_detail::check_dist_and_k(

           function, trials, RealType(0), successes, &result, Policy())

            &&

           binomial_detail::check_dist_and_prob(

           function, trials, RealType(0), probability, &result, Policy()))

        { return result; }

        if(successes == 0)

           return 0;

        // NOTE!!! The Clopper Pearson formula uses "successes" not

        // "successes+1" as usual to get the lower bound,

        // see http://www.itl.nist.gov/div898/handbook/prc/section2/prc241.htm

        return (t == clopper_pearson_exact_interval) ? ibeta_inv(successes, trials - successes + 1, probability, static_cast<RealType*>(0), Policy())

           : ibeta_inv(successes + 0.5f, trials - successes + 0.5f, probability, static_cast<RealType*>(0), Policy());

      }

      static RealType find_upper_bound_on_p(

         RealType trials,

         RealType successes,

         RealType probability,

         interval_type t = clopper_pearson_exact_interval)

      {

        static const char* function = "boost::math::binomial_distribution<%1%>::find_upper_bound_on_p";

        // Error checks:

        RealType result = 0;

        if(false == binomial_detail::check_dist_and_k(

           function, trials, RealType(0), successes, &result, Policy())

            &&

           binomial_detail::check_dist_and_prob(

           function, trials, RealType(0), probability, &result, Policy()))

        { return result; }

        if(trials == successes)

           return 1;

        return (t == clopper_pearson_exact_interval) ? ibetac_inv(successes + 1, trials - successes, probability, static_cast<RealType*>(0), Policy())

           : ibetac_inv(successes + 0.5f, trials - successes + 0.5f, probability, static_cast<RealType*>(0), Policy());

      }

      // Estimate number of trials parameter:

      //

      // "How many trials do I need to be P% sure of seeing k events?"

      //    or

      // "How many trials can I have to be P% sure of seeing fewer than k events?"

      //

      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 = 0;

        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 = 0;

        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>(static_cast<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 = 0; // initialization silences some compiler warnings

        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 = 0;

        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 = 0;

        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, RealType(1-p), false);

      } // quantile

      template <class RealType, class Policy>

      RealType quantile(const complemented2_type<binomial_distribution<RealType, Policy>, RealType>& c)

      {

         return binomial_detail::quantile_imp(c.dist, RealType(1-c.param), c.param, true);

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