Chris@16: // Copyright 2008 Gautam Sewani
Chris@16: // Copyright 2008 John Maddock
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_DISTRIBUTIONS_HYPERGEOMETRIC_HPP
Chris@16: #define BOOST_MATH_DISTRIBUTIONS_HYPERGEOMETRIC_HPP
Chris@16: 
Chris@16: #include <boost/math/distributions/detail/common_error_handling.hpp>
Chris@16: #include <boost/math/distributions/complement.hpp>
Chris@16: #include <boost/math/distributions/detail/hypergeometric_pdf.hpp>
Chris@16: #include <boost/math/distributions/detail/hypergeometric_cdf.hpp>
Chris@16: #include <boost/math/distributions/detail/hypergeometric_quantile.hpp>
Chris@16: #include <boost/math/special_functions/fpclassify.hpp>
Chris@16: 
Chris@16: 
Chris@16: namespace boost { namespace math {
Chris@16: 
Chris@16:    template <class RealType = double, class Policy = policies::policy<> >
Chris@16:    class hypergeometric_distribution
Chris@16:    {
Chris@16:    public:
Chris@16:       typedef RealType value_type;
Chris@16:       typedef Policy policy_type;
Chris@16: 
Chris@16:       hypergeometric_distribution(unsigned r, unsigned n, unsigned N) // Constructor.
Chris@16:          : m_n(n), m_N(N), m_r(r)
Chris@16:       {
Chris@16:          static const char* function = "boost::math::hypergeometric_distribution<%1%>::hypergeometric_distribution";
Chris@16:          RealType ret;
Chris@16:          check_params(function, &ret);
Chris@16:       }
Chris@16:       // Accessor functions.
Chris@16:       unsigned total()const
Chris@16:       {
Chris@16:          return m_N;
Chris@16:       }
Chris@16: 
Chris@16:       unsigned defective()const
Chris@16:       {
Chris@16:          return m_n;
Chris@16:       }
Chris@16: 
Chris@16:       unsigned sample_count()const
Chris@16:       {
Chris@16:          return m_r;
Chris@16:       }
Chris@16: 
Chris@16:       bool check_params(const char* function, RealType* result)const
Chris@16:       {
Chris@16:          if(m_r > m_N)
Chris@16:          {
Chris@16:             *result = boost::math::policies::raise_domain_error<RealType>(
Chris@16:                function, "Parameter r out of range: must be <= N but got %1%", static_cast<RealType>(m_r), Policy());
Chris@16:             return false;
Chris@16:          }
Chris@16:          if(m_n > m_N)
Chris@16:          {
Chris@16:             *result = boost::math::policies::raise_domain_error<RealType>(
Chris@16:                function, "Parameter n out of range: must be <= N but got %1%", static_cast<RealType>(m_n), Policy());
Chris@16:             return false;
Chris@16:          }
Chris@16:          return true;
Chris@16:       }
Chris@16:       bool check_x(unsigned x, const char* function, RealType* result)const
Chris@16:       {
Chris@16:          if(x < static_cast<unsigned>((std::max)(0, (int)(m_n + m_r) - (int)(m_N))))
Chris@16:          {
Chris@16:             *result = boost::math::policies::raise_domain_error<RealType>(
Chris@16:                function, "Random variable out of range: must be > 0 and > m + r - N but got %1%", static_cast<RealType>(x), Policy());
Chris@16:             return false;
Chris@16:          }
Chris@16:          if(x > (std::min)(m_r, m_n))
Chris@16:          {
Chris@16:             *result = boost::math::policies::raise_domain_error<RealType>(
Chris@16:                function, "Random variable out of range: must be less than both n and r but got %1%", static_cast<RealType>(x), Policy());
Chris@16:             return false;
Chris@16:          }
Chris@16:          return true;
Chris@16:       }
Chris@16: 
Chris@16:    private:
Chris@16:       // Data members:
Chris@16:       unsigned m_n;  // number of "defective" items
Chris@16:       unsigned m_N; // number of "total" items
Chris@16:       unsigned m_r; // number of items picked
Chris@16: 
Chris@16:    }; // class hypergeometric_distribution
Chris@16: 
Chris@16:    typedef hypergeometric_distribution<double> hypergeometric;
Chris@16: 
Chris@16:    template <class RealType, class Policy>
Chris@16:    inline const std::pair<unsigned, unsigned> range(const hypergeometric_distribution<RealType, Policy>& dist)
Chris@16:    { // Range of permissible values for random variable x.
Chris@16: #ifdef BOOST_MSVC
Chris@16: #  pragma warning(push)
Chris@16: #  pragma warning(disable:4267)
Chris@16: #endif
Chris@16:       unsigned r = dist.sample_count();
Chris@16:       unsigned n = dist.defective();
Chris@16:       unsigned N = dist.total();
Chris@16:       unsigned l = static_cast<unsigned>((std::max)(0, (int)(n + r) - (int)(N)));
Chris@16:       unsigned u = (std::min)(r, n);
Chris@16:       return std::pair<unsigned, unsigned>(l, u);
Chris@16: #ifdef BOOST_MSVC
Chris@16: #  pragma warning(pop)
Chris@16: #endif
Chris@16:    }
Chris@16: 
Chris@16:    template <class RealType, class Policy>
Chris@16:    inline const std::pair<unsigned, unsigned> support(const hypergeometric_distribution<RealType, Policy>& d)
Chris@16:    { 
Chris@16:       return range(d);
Chris@16:    }
Chris@16: 
Chris@16:    template <class RealType, class Policy>
Chris@16:    inline RealType pdf(const hypergeometric_distribution<RealType, Policy>& dist, const unsigned& x)
Chris@16:    {
Chris@16:       static const char* function = "boost::math::pdf(const hypergeometric_distribution<%1%>&, const %1%&)";
Chris@16:       RealType result = 0;
Chris@16:       if(!dist.check_params(function, &result))
Chris@16:          return result;
Chris@16:       if(!dist.check_x(x, function, &result))
Chris@16:          return result;
Chris@16: 
Chris@16:       return boost::math::detail::hypergeometric_pdf<RealType>(
Chris@16:          x, dist.sample_count(), dist.defective(), dist.total(), Policy());
Chris@16:    }
Chris@16: 
Chris@16:    template <class RealType, class Policy, class U>
Chris@16:    inline RealType pdf(const hypergeometric_distribution<RealType, Policy>& dist, const U& x)
Chris@16:    {
Chris@16:       BOOST_MATH_STD_USING
Chris@16:       static const char* function = "boost::math::pdf(const hypergeometric_distribution<%1%>&, const %1%&)";
Chris@16:       RealType r = static_cast<RealType>(x);
Chris@16:       unsigned u = itrunc(r, typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type());
Chris@16:       if(u != r)
Chris@16:       {
Chris@16:          return boost::math::policies::raise_domain_error<RealType>(
Chris@16:             function, "Random variable out of range: must be an integer but got %1%", r, Policy());
Chris@16:       }
Chris@16:       return pdf(dist, u);
Chris@16:    }
Chris@16: 
Chris@16:    template <class RealType, class Policy>
Chris@16:    inline RealType cdf(const hypergeometric_distribution<RealType, Policy>& dist, const unsigned& x)
Chris@16:    {
Chris@16:       static const char* function = "boost::math::cdf(const hypergeometric_distribution<%1%>&, const %1%&)";
Chris@16:       RealType result = 0;
Chris@16:       if(!dist.check_params(function, &result))
Chris@16:          return result;
Chris@16:       if(!dist.check_x(x, function, &result))
Chris@16:          return result;
Chris@16: 
Chris@16:       return boost::math::detail::hypergeometric_cdf<RealType>(
Chris@16:          x, dist.sample_count(), dist.defective(), dist.total(), false, Policy());
Chris@16:    }
Chris@16: 
Chris@16:    template <class RealType, class Policy, class U>
Chris@16:    inline RealType cdf(const hypergeometric_distribution<RealType, Policy>& dist, const U& x)
Chris@16:    {
Chris@16:       BOOST_MATH_STD_USING
Chris@16:       static const char* function = "boost::math::cdf(const hypergeometric_distribution<%1%>&, const %1%&)";
Chris@16:       RealType r = static_cast<RealType>(x);
Chris@16:       unsigned u = itrunc(r, typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type());
Chris@16:       if(u != r)
Chris@16:       {
Chris@16:          return boost::math::policies::raise_domain_error<RealType>(
Chris@16:             function, "Random variable out of range: must be an integer but got %1%", r, Policy());
Chris@16:       }
Chris@16:       return cdf(dist, u);
Chris@16:    }
Chris@16: 
Chris@16:    template <class RealType, class Policy>
Chris@16:    inline RealType cdf(const complemented2_type<hypergeometric_distribution<RealType, Policy>, unsigned>& c)
Chris@16:    {
Chris@16:       static const char* function = "boost::math::cdf(const hypergeometric_distribution<%1%>&, const %1%&)";
Chris@16:       RealType result = 0;
Chris@16:       if(!c.dist.check_params(function, &result))
Chris@16:          return result;
Chris@16:       if(!c.dist.check_x(c.param, function, &result))
Chris@16:          return result;
Chris@16: 
Chris@16:       return boost::math::detail::hypergeometric_cdf<RealType>(
Chris@16:          c.param, c.dist.sample_count(), c.dist.defective(), c.dist.total(), true, Policy());
Chris@16:    }
Chris@16: 
Chris@16:    template <class RealType, class Policy, class U>
Chris@16:    inline RealType cdf(const complemented2_type<hypergeometric_distribution<RealType, Policy>, U>& c)
Chris@16:    {
Chris@16:       BOOST_MATH_STD_USING
Chris@16:       static const char* function = "boost::math::cdf(const hypergeometric_distribution<%1%>&, const %1%&)";
Chris@16:       RealType r = static_cast<RealType>(c.param);
Chris@16:       unsigned u = itrunc(r, typename policies::normalise<Policy, policies::rounding_error<policies::ignore_error> >::type());
Chris@16:       if(u != r)
Chris@16:       {
Chris@16:          return boost::math::policies::raise_domain_error<RealType>(
Chris@16:             function, "Random variable out of range: must be an integer but got %1%", r, Policy());
Chris@16:       }
Chris@16:       return cdf(complement(c.dist, u));
Chris@16:    }
Chris@16: 
Chris@16:    template <class RealType, class Policy>
Chris@16:    inline RealType quantile(const hypergeometric_distribution<RealType, Policy>& dist, const RealType& p)
Chris@16:    {
Chris@16:       BOOST_MATH_STD_USING // for ADL of std functions
Chris@16: 
Chris@16:          // Checking function argument
Chris@16:          RealType result = 0;
Chris@16:       const char* function = "boost::math::quantile(const hypergeometric_distribution<%1%>&, %1%)";
Chris@16:       if (false == dist.check_params(function, &result)) return result;
Chris@16:       if(false == detail::check_probability(function, p, &result, Policy())) return result;
Chris@16: 
Chris@16:       return static_cast<RealType>(detail::hypergeometric_quantile(p, RealType(1 - p), dist.sample_count(), dist.defective(), dist.total(), Policy()));
Chris@16:    } // quantile
Chris@16: 
Chris@16:    template <class RealType, class Policy>
Chris@16:    inline RealType quantile(const complemented2_type<hypergeometric_distribution<RealType, Policy>, RealType>& c)
Chris@16:    {
Chris@16:       BOOST_MATH_STD_USING // for ADL of std functions
Chris@16: 
Chris@16:       // Checking function argument
Chris@16:       RealType result = 0;
Chris@16:       const char* function = "quantile(const complemented2_type<hypergeometric_distribution<%1%>, %1%>&)";
Chris@16:       if (false == c.dist.check_params(function, &result)) return result;
Chris@16:       if(false == detail::check_probability(function, c.param, &result, Policy())) return result;
Chris@16: 
Chris@16:       return static_cast<RealType>(detail::hypergeometric_quantile(RealType(1 - c.param), c.param, c.dist.sample_count(), c.dist.defective(), c.dist.total(), Policy()));
Chris@16:    } // quantile
Chris@16: 
Chris@16:    template <class RealType, class Policy>
Chris@16:    inline RealType mean(const hypergeometric_distribution<RealType, Policy>& dist)
Chris@16:    {
Chris@16:       return static_cast<RealType>(dist.sample_count() * dist.defective()) / dist.total();
Chris@16:    } // RealType mean(const hypergeometric_distribution<RealType, Policy>& dist)
Chris@16: 
Chris@16:    template <class RealType, class Policy>
Chris@16:    inline RealType variance(const hypergeometric_distribution<RealType, Policy>& dist)
Chris@16:    {
Chris@16:       RealType r = static_cast<RealType>(dist.sample_count());
Chris@16:       RealType n = static_cast<RealType>(dist.defective());
Chris@16:       RealType N = static_cast<RealType>(dist.total());
Chris@16:       return r * (n / N) * (1 - n / N) * (N - r) / (N - 1);
Chris@16:    } // RealType variance(const hypergeometric_distribution<RealType, Policy>& dist)
Chris@16: 
Chris@16:    template <class RealType, class Policy>
Chris@16:    inline RealType mode(const hypergeometric_distribution<RealType, Policy>& dist)
Chris@16:    {
Chris@16:       BOOST_MATH_STD_USING
Chris@16:       RealType r = static_cast<RealType>(dist.sample_count());
Chris@16:       RealType n = static_cast<RealType>(dist.defective());
Chris@16:       RealType N = static_cast<RealType>(dist.total());
Chris@16:       return floor((r + 1) * (n + 1) / (N + 2));
Chris@16:    }
Chris@16: 
Chris@16:    template <class RealType, class Policy>
Chris@16:    inline RealType skewness(const hypergeometric_distribution<RealType, Policy>& dist)
Chris@16:    {
Chris@16:       BOOST_MATH_STD_USING
Chris@16:       RealType r = static_cast<RealType>(dist.sample_count());
Chris@16:       RealType n = static_cast<RealType>(dist.defective());
Chris@16:       RealType N = static_cast<RealType>(dist.total());
Chris@16:       return (N - 2 * n) * sqrt(N - 1) * (N - 2 * r) / (sqrt(n * r * (N - n) * (N - r)) * (N - 2));
Chris@16:    } // RealType skewness(const hypergeometric_distribution<RealType, Policy>& dist)
Chris@16: 
Chris@16:    template <class RealType, class Policy>
Chris@16:    inline RealType kurtosis_excess(const hypergeometric_distribution<RealType, Policy>& dist)
Chris@16:    {
Chris@16:       RealType r = static_cast<RealType>(dist.sample_count());
Chris@16:       RealType n = static_cast<RealType>(dist.defective());
Chris@16:       RealType N = static_cast<RealType>(dist.total());
Chris@16:       RealType t1 = N * N * (N - 1) / (r * (N - 2) * (N - 3) * (N - r));
Chris@16:       RealType t2 = (N * (N + 1) - 6 * N * (N - r)) / (n * (N - n))
Chris@16:          + 3 * r * (N - r) * (N + 6) / (N * N) - 6;
Chris@16:       return t1 * t2;
Chris@16:    } // RealType kurtosis_excess(const hypergeometric_distribution<RealType, Policy>& dist)
Chris@16: 
Chris@16:    template <class RealType, class Policy>
Chris@16:    inline RealType kurtosis(const hypergeometric_distribution<RealType, Policy>& dist)
Chris@16:    {
Chris@16:       return kurtosis_excess(dist) + 3;
Chris@16:    } // RealType kurtosis_excess(const hypergeometric_distribution<RealType, Policy>& dist)
Chris@16: }} // namespaces
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 // include guard