Chris@102: //  Copyright 2014 Marco Guazzone (marco.guazzone@gmail.com)
Chris@102: //
Chris@102: //  Use, modification and distribution are subject to the
Chris@102: //  Boost Software License, Version 1.0. (See accompanying file
Chris@102: //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
Chris@102: //
Chris@102: // This module implements the Hyper-Exponential distribution.
Chris@102: //
Chris@102: // References:
Chris@102: // - "Queueing Theory in Manufacturing Systems Analysis and Design" by H.T. Papadopolous, C. Heavey and J. Browne (Chapman & Hall/CRC, 1993)
Chris@102: // - http://reference.wolfram.com/language/ref/HyperexponentialDistribution.html
Chris@102: // - http://en.wikipedia.org/wiki/Hyperexponential_distribution
Chris@102: //
Chris@102: 
Chris@102: #ifndef BOOST_MATH_DISTRIBUTIONS_HYPEREXPONENTIAL_HPP
Chris@102: #define BOOST_MATH_DISTRIBUTIONS_HYPEREXPONENTIAL_HPP
Chris@102: 
Chris@102: 
Chris@102: #include <boost/config.hpp>
Chris@102: #include <boost/math/distributions/complement.hpp>
Chris@102: #include <boost/math/distributions/detail/common_error_handling.hpp>
Chris@102: #include <boost/math/distributions/exponential.hpp>
Chris@102: #include <boost/math/policies/policy.hpp>
Chris@102: #include <boost/math/special_functions/fpclassify.hpp>
Chris@102: #include <boost/math/tools/precision.hpp>
Chris@102: #include <boost/math/tools/roots.hpp>
Chris@102: #include <boost/range/begin.hpp>
Chris@102: #include <boost/range/end.hpp>
Chris@102: #include <boost/range/size.hpp>
Chris@102: #include <boost/type_traits/has_pre_increment.hpp>
Chris@102: #include <cstddef>
Chris@102: #include <iterator>
Chris@102: #include <limits>
Chris@102: #include <numeric>
Chris@102: #include <utility>
Chris@102: #include <vector>
Chris@102: 
Chris@102: #if !defined(BOOST_NO_CXX11_HDR_INITIALIZER_LIST)
Chris@102: # include <initializer_list>
Chris@102: #endif
Chris@102: 
Chris@102: #ifdef _MSC_VER
Chris@102: # pragma warning (push)
Chris@102: # pragma warning(disable:4127) // conditional expression is constant
Chris@102: # pragma warning(disable:4389) // '==' : signed/unsigned mismatch in test_tools
Chris@102: #endif // _MSC_VER
Chris@102: 
Chris@102: namespace boost { namespace math {
Chris@102: 
Chris@102: namespace detail {
Chris@102: 
Chris@102: template <typename Dist>
Chris@102: typename Dist::value_type generic_quantile(const Dist& dist, const typename Dist::value_type& p, const typename Dist::value_type& guess, bool comp, const char* function);
Chris@102: 
Chris@102: } // Namespace detail
Chris@102: 
Chris@102: 
Chris@102: template <typename RealT, typename PolicyT>
Chris@102: class hyperexponential_distribution;
Chris@102: 
Chris@102: 
Chris@102: namespace /*<unnamed>*/ { namespace hyperexp_detail {
Chris@102: 
Chris@102: template <typename T>
Chris@102: void normalize(std::vector<T>& v)
Chris@102: {
Chris@102:    if(!v.size())
Chris@102:       return;  // Our error handlers will get this later
Chris@102:     const T sum = std::accumulate(v.begin(), v.end(), static_cast<T>(0));
Chris@102:     T final_sum = 0;
Chris@102:     const typename std::vector<T>::iterator end = --v.end();
Chris@102:     for (typename std::vector<T>::iterator it = v.begin();
Chris@102:          it != end;
Chris@102:          ++it)
Chris@102:     {
Chris@102:         *it /= sum;
Chris@102:         final_sum += *it;
Chris@102:     }
Chris@102:     *end = 1 - final_sum;  // avoids round off errors, ensures the probs really do sum to 1.
Chris@102: }
Chris@102: 
Chris@102: template <typename RealT, typename PolicyT>
Chris@102: bool check_probabilities(char const* function, std::vector<RealT> const& probabilities, RealT* presult, PolicyT const& pol)
Chris@102: {
Chris@102:     BOOST_MATH_STD_USING
Chris@102:     const std::size_t n = probabilities.size();
Chris@102:     RealT sum = 0;
Chris@102:     for (std::size_t i = 0; i < n; ++i)
Chris@102:     {
Chris@102:         if (probabilities[i] < 0
Chris@102:             || probabilities[i] > 1
Chris@102:             || !(boost::math::isfinite)(probabilities[i]))
Chris@102:         {
Chris@102:             *presult = policies::raise_domain_error<RealT>(function,
Chris@102:                                                            "The elements of parameter \"probabilities\" must be >= 0 and <= 1, but at least one of them was: %1%.",
Chris@102:                                                            probabilities[i],
Chris@102:                                                            pol);
Chris@102:             return false;
Chris@102:         }
Chris@102:         sum += probabilities[i];
Chris@102:     }
Chris@102: 
Chris@102:     //
Chris@102:     // We try to keep phase probabilities correctly normalized in the distribution constructors,
Chris@102:     // however in practice we have to allow for a very slight divergence from a sum of exactly 1:
Chris@102:     //
Chris@102:     if (fabs(sum - 1) > tools::epsilon<RealT>() * 2)
Chris@102:     {
Chris@102:         *presult = policies::raise_domain_error<RealT>(function,
Chris@102:                                                        "The elements of parameter \"probabilities\" must sum to 1, but their sum is: %1%.",
Chris@102:                                                        sum,
Chris@102:                                                        pol);
Chris@102:         return false;
Chris@102:     }
Chris@102: 
Chris@102:     return true;
Chris@102: }
Chris@102: 
Chris@102: template <typename RealT, typename PolicyT>
Chris@102: bool check_rates(char const* function, std::vector<RealT> const& rates, RealT* presult, PolicyT const& pol)
Chris@102: {
Chris@102:     const std::size_t n = rates.size();
Chris@102:     for (std::size_t i = 0; i < n; ++i)
Chris@102:     {
Chris@102:         if (rates[i] <= 0
Chris@102:             || !(boost::math::isfinite)(rates[i]))
Chris@102:         {
Chris@102:             *presult = policies::raise_domain_error<RealT>(function,
Chris@102:                                                            "The elements of parameter \"rates\" must be > 0, but at least one of them is: %1%.",
Chris@102:                                                            rates[i],
Chris@102:                                                            pol);
Chris@102:             return false;
Chris@102:         }
Chris@102:     }
Chris@102:     return true;
Chris@102: }
Chris@102: 
Chris@102: template <typename RealT, typename PolicyT>
Chris@102: bool check_dist(char const* function, std::vector<RealT> const& probabilities, std::vector<RealT> const& rates, RealT* presult, PolicyT const& pol)
Chris@102: {
Chris@102:     BOOST_MATH_STD_USING
Chris@102:     if (probabilities.size() != rates.size())
Chris@102:     {
Chris@102:         *presult = policies::raise_domain_error<RealT>(function,
Chris@102:                                                        "The parameters \"probabilities\" and \"rates\" must have the same length, but their size differ by: %1%.",
Chris@102:                                                        fabs(static_cast<RealT>(probabilities.size())-static_cast<RealT>(rates.size())),
Chris@102:                                                        pol);
Chris@102:         return false;
Chris@102:     }
Chris@102: 
Chris@102:     return check_probabilities(function, probabilities, presult, pol)
Chris@102:            && check_rates(function, rates, presult, pol);
Chris@102: }
Chris@102: 
Chris@102: template <typename RealT, typename PolicyT>
Chris@102: bool check_x(char const* function, RealT x, RealT* presult, PolicyT const& pol)
Chris@102: {
Chris@102:     if (x < 0 || (boost::math::isnan)(x))
Chris@102:     {
Chris@102:         *presult = policies::raise_domain_error<RealT>(function, "The random variable must be >= 0, but is: %1%.", x, pol);
Chris@102:         return false;
Chris@102:     }
Chris@102:     return true;
Chris@102: }
Chris@102: 
Chris@102: template <typename RealT, typename PolicyT>
Chris@102: bool check_probability(char const* function, RealT p, RealT* presult, PolicyT const& pol)
Chris@102: {
Chris@102:     if (p < 0 || p > 1 || (boost::math::isnan)(p))
Chris@102:     {
Chris@102:         *presult = policies::raise_domain_error<RealT>(function, "The probability be >= 0 and <= 1, but is: %1%.", p, pol);
Chris@102:         return false;
Chris@102:     }
Chris@102:     return true;
Chris@102: }
Chris@102: 
Chris@102: template <typename RealT, typename PolicyT>
Chris@102: RealT quantile_impl(hyperexponential_distribution<RealT, PolicyT> const& dist, RealT const& p, bool comp)
Chris@102: {
Chris@102:     // Don't have a closed form so try to numerically solve the inverse CDF...
Chris@102: 
Chris@102:     typedef typename policies::evaluation<RealT, PolicyT>::type value_type;
Chris@102:     typedef typename policies::normalise<PolicyT,
Chris@102:                                          policies::promote_float<false>,
Chris@102:                                          policies::promote_double<false>,
Chris@102:                                          policies::discrete_quantile<>,
Chris@102:                                          policies::assert_undefined<> >::type forwarding_policy;
Chris@102: 
Chris@102:     static const char* function = comp ? "boost::math::quantile(const boost::math::complemented2_type<boost::math::hyperexponential_distribution<%1%>, %1%>&)"
Chris@102:                                        : "boost::math::quantile(const boost::math::hyperexponential_distribution<%1%>&, %1%)";
Chris@102: 
Chris@102:     RealT result = 0;
Chris@102: 
Chris@102:     if (!check_probability(function, p, &result, PolicyT()))
Chris@102:     {
Chris@102:         return result;
Chris@102:     }
Chris@102: 
Chris@102:     const std::size_t n = dist.num_phases();
Chris@102:     const std::vector<RealT> probs = dist.probabilities();
Chris@102:     const std::vector<RealT> rates = dist.rates();
Chris@102: 
Chris@102:     // A possible (but inaccurate) approximation is given below, where the
Chris@102:     // quantile is given by the weighted sum of exponential quantiles:
Chris@102:     RealT guess = 0;
Chris@102:     if (comp)
Chris@102:     {
Chris@102:         for (std::size_t i = 0; i < n; ++i)
Chris@102:         {
Chris@102:             const exponential_distribution<RealT,PolicyT> exp(rates[i]);
Chris@102: 
Chris@102:             guess += probs[i]*quantile(complement(exp, p));
Chris@102:         }
Chris@102:     }
Chris@102:     else
Chris@102:     {
Chris@102:         for (std::size_t i = 0; i < n; ++i)
Chris@102:         {
Chris@102:             const exponential_distribution<RealT,PolicyT> exp(rates[i]);
Chris@102: 
Chris@102:             guess += probs[i]*quantile(exp, p);
Chris@102:         }
Chris@102:     }
Chris@102: 
Chris@102:     // Fast return in case the Hyper-Exponential is essentially an Exponential
Chris@102:     if (n == 1)
Chris@102:     {
Chris@102:         return guess;
Chris@102:     }
Chris@102: 
Chris@102:     value_type q;
Chris@102:     q = detail::generic_quantile(hyperexponential_distribution<RealT,forwarding_policy>(probs, rates),
Chris@102:                                  p,
Chris@102:                                  guess,
Chris@102:                                  comp,
Chris@102:                                  function);
Chris@102: 
Chris@102:     result = policies::checked_narrowing_cast<RealT,forwarding_policy>(q, function);
Chris@102: 
Chris@102:     return result;
Chris@102: }
Chris@102: 
Chris@102: }} // Namespace <unnamed>::hyperexp_detail
Chris@102: 
Chris@102: 
Chris@102: template <typename RealT = double, typename PolicyT = policies::policy<> >
Chris@102: class hyperexponential_distribution
Chris@102: {
Chris@102:     public: typedef RealT value_type;
Chris@102:     public: typedef PolicyT policy_type;
Chris@102: 
Chris@102: 
Chris@102:     public: hyperexponential_distribution()
Chris@102:     : probs_(1, 1),
Chris@102:       rates_(1, 1)
Chris@102:     {
Chris@102:         RealT err;
Chris@102:         hyperexp_detail::check_dist("boost::math::hyperexponential_distribution<%1%>::hyperexponential_distribution",
Chris@102:                                     probs_,
Chris@102:                                     rates_,
Chris@102:                                     &err,
Chris@102:                                     PolicyT());
Chris@102:     }
Chris@102: 
Chris@102:     // Four arg constructor: no ambiguity here, the arguments must be two pairs of iterators:
Chris@102:     public: template <typename ProbIterT, typename RateIterT>
Chris@102:             hyperexponential_distribution(ProbIterT prob_first, ProbIterT prob_last,
Chris@102:                                           RateIterT rate_first, RateIterT rate_last)
Chris@102:     : probs_(prob_first, prob_last),
Chris@102:       rates_(rate_first, rate_last)
Chris@102:     {
Chris@102:         hyperexp_detail::normalize(probs_);
Chris@102:         RealT err;
Chris@102:         hyperexp_detail::check_dist("boost::math::hyperexponential_distribution<%1%>::hyperexponential_distribution",
Chris@102:                                     probs_,
Chris@102:                                     rates_,
Chris@102:                                     &err,
Chris@102:                                     PolicyT());
Chris@102:     }
Chris@102: 
Chris@102:     // Two arg constructor from 2 ranges, we SFINAE this out of existance if
Chris@102:     // either argument type is incrementable as in that case the type is
Chris@102:     // probably an iterator:
Chris@102:     public: template <typename ProbRangeT, typename RateRangeT>
Chris@102:             hyperexponential_distribution(ProbRangeT const& prob_range,
Chris@102:                                           RateRangeT const& rate_range,
Chris@102:                                           typename boost::disable_if_c<boost::has_pre_increment<ProbRangeT>::value || boost::has_pre_increment<RateRangeT>::value>::type* = 0)
Chris@102:     : probs_(boost::begin(prob_range), boost::end(prob_range)),
Chris@102:       rates_(boost::begin(rate_range), boost::end(rate_range))
Chris@102:     {
Chris@102:         hyperexp_detail::normalize(probs_);
Chris@102: 
Chris@102:         RealT err;
Chris@102:         hyperexp_detail::check_dist("boost::math::hyperexponential_distribution<%1%>::hyperexponential_distribution",
Chris@102:                                     probs_,
Chris@102:                                     rates_,
Chris@102:                                     &err,
Chris@102:                                     PolicyT());
Chris@102:     }
Chris@102: 
Chris@102:     // Two arg constructor for a pair of iterators: we SFINAE this out of
Chris@102:     // existance if neither argument types are incrementable.
Chris@102:     // Note that we allow different argument types here to allow for
Chris@102:     // construction from an array plus a pointer into that array.
Chris@102:     public: template <typename RateIterT, typename RateIterT2>
Chris@102:             hyperexponential_distribution(RateIterT const& rate_first, 
Chris@102:                                           RateIterT2 const& rate_last, 
Chris@102:                                           typename boost::enable_if_c<boost::has_pre_increment<RateIterT>::value || boost::has_pre_increment<RateIterT2>::value>::type* = 0)
Chris@102:     : probs_(std::distance(rate_first, rate_last), 1), // will be normalized below
Chris@102:       rates_(rate_first, rate_last)
Chris@102:     {
Chris@102:         hyperexp_detail::normalize(probs_);
Chris@102: 
Chris@102:         RealT err;
Chris@102:         hyperexp_detail::check_dist("boost::math::hyperexponential_distribution<%1%>::hyperexponential_distribution",
Chris@102:                                     probs_,
Chris@102:                                     rates_,
Chris@102:                                     &err,
Chris@102:                                     PolicyT());
Chris@102:     }
Chris@102: 
Chris@102: #if !defined(BOOST_NO_CXX11_HDR_INITIALIZER_LIST)
Chris@102:       // Initializer list constructor: allows for construction from array literals:
Chris@102: public: hyperexponential_distribution(std::initializer_list<RealT> l1, std::initializer_list<RealT> l2)
Chris@102:       : probs_(l1.begin(), l1.end()),
Chris@102:         rates_(l2.begin(), l2.end())
Chris@102:       {
Chris@102:          hyperexp_detail::normalize(probs_);
Chris@102: 
Chris@102:          RealT err;
Chris@102:          hyperexp_detail::check_dist("boost::math::hyperexponential_distribution<%1%>::hyperexponential_distribution",
Chris@102:             probs_,
Chris@102:             rates_,
Chris@102:             &err,
Chris@102:             PolicyT());
Chris@102:       }
Chris@102: 
Chris@102: public: hyperexponential_distribution(std::initializer_list<RealT> l1)
Chris@102:       : probs_(l1.size(), 1),
Chris@102:         rates_(l1.begin(), l1.end())
Chris@102:       {
Chris@102:          hyperexp_detail::normalize(probs_);
Chris@102: 
Chris@102:          RealT err;
Chris@102:          hyperexp_detail::check_dist("boost::math::hyperexponential_distribution<%1%>::hyperexponential_distribution",
Chris@102:             probs_,
Chris@102:             rates_,
Chris@102:             &err,
Chris@102:             PolicyT());
Chris@102:       }
Chris@102: #endif // !defined(BOOST_NO_CXX11_HDR_INITIALIZER_LIST)
Chris@102: 
Chris@102:     // Single argument constructor: argument must be a range.
Chris@102:     public: template <typename RateRangeT>
Chris@102:     hyperexponential_distribution(RateRangeT const& rate_range)
Chris@102:     : probs_(boost::size(rate_range), 1), // will be normalized below
Chris@102:       rates_(boost::begin(rate_range), boost::end(rate_range))
Chris@102:     {
Chris@102:         hyperexp_detail::normalize(probs_);
Chris@102: 
Chris@102:         RealT err;
Chris@102:         hyperexp_detail::check_dist("boost::math::hyperexponential_distribution<%1%>::hyperexponential_distribution",
Chris@102:                                     probs_,
Chris@102:                                     rates_,
Chris@102:                                     &err,
Chris@102:                                     PolicyT());
Chris@102:     }
Chris@102: 
Chris@102:     public: std::vector<RealT> probabilities() const
Chris@102:     {
Chris@102:         return probs_;
Chris@102:     }
Chris@102: 
Chris@102:     public: std::vector<RealT> rates() const
Chris@102:     {
Chris@102:         return rates_;
Chris@102:     }
Chris@102: 
Chris@102:     public: std::size_t num_phases() const
Chris@102:     {
Chris@102:         return rates_.size();
Chris@102:     }
Chris@102: 
Chris@102: 
Chris@102:     private: std::vector<RealT> probs_;
Chris@102:     private: std::vector<RealT> rates_;
Chris@102: }; // class hyperexponential_distribution
Chris@102: 
Chris@102: 
Chris@102: // Convenient type synonym for double.
Chris@102: typedef hyperexponential_distribution<double> hyperexponential;
Chris@102: 
Chris@102: 
Chris@102: // Range of permissible values for random variable x
Chris@102: template <typename RealT, typename PolicyT>
Chris@102: std::pair<RealT,RealT> range(hyperexponential_distribution<RealT,PolicyT> const&)
Chris@102: {
Chris@102:     if (std::numeric_limits<RealT>::has_infinity)
Chris@102:     {
Chris@102:         return std::make_pair(static_cast<RealT>(0), std::numeric_limits<RealT>::infinity()); // 0 to +inf.
Chris@102:     }
Chris@102: 
Chris@102:     return std::make_pair(static_cast<RealT>(0), tools::max_value<RealT>()); // 0 to +<max value>
Chris@102: }
Chris@102: 
Chris@102: // Range of supported values for random variable x.
Chris@102: // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
Chris@102: template <typename RealT, typename PolicyT>
Chris@102: std::pair<RealT,RealT> support(hyperexponential_distribution<RealT,PolicyT> const&)
Chris@102: {
Chris@102:     return std::make_pair(tools::min_value<RealT>(), tools::max_value<RealT>()); // <min value> to +<max value>.
Chris@102: }
Chris@102: 
Chris@102: template <typename RealT, typename PolicyT>
Chris@102: RealT pdf(hyperexponential_distribution<RealT, PolicyT> const& dist, RealT const& x)
Chris@102: {
Chris@102:     BOOST_MATH_STD_USING
Chris@102:     RealT result = 0;
Chris@102: 
Chris@102:     if (!hyperexp_detail::check_x("boost::math::pdf(const boost::math::hyperexponential_distribution<%1%>&, %1%)", x, &result, PolicyT()))
Chris@102:     {
Chris@102:         return result;
Chris@102:     }
Chris@102: 
Chris@102:     const std::size_t n = dist.num_phases();
Chris@102:     const std::vector<RealT> probs = dist.probabilities();
Chris@102:     const std::vector<RealT> rates = dist.rates();
Chris@102: 
Chris@102:     for (std::size_t i = 0; i < n; ++i)
Chris@102:     {
Chris@102:         const exponential_distribution<RealT,PolicyT> exp(rates[i]);
Chris@102: 
Chris@102:         result += probs[i]*pdf(exp, x);
Chris@102:         //result += probs[i]*rates[i]*exp(-rates[i]*x);
Chris@102:     }
Chris@102: 
Chris@102:     return result;
Chris@102: }
Chris@102: 
Chris@102: template <typename RealT, typename PolicyT>
Chris@102: RealT cdf(hyperexponential_distribution<RealT, PolicyT> const& dist, RealT const& x)
Chris@102: {
Chris@102:     RealT result = 0;
Chris@102: 
Chris@102:     if (!hyperexp_detail::check_x("boost::math::cdf(const boost::math::hyperexponential_distribution<%1%>&, %1%)", x, &result, PolicyT()))
Chris@102:     {
Chris@102:         return result;
Chris@102:     }
Chris@102: 
Chris@102:     const std::size_t n = dist.num_phases();
Chris@102:     const std::vector<RealT> probs = dist.probabilities();
Chris@102:     const std::vector<RealT> rates = dist.rates();
Chris@102: 
Chris@102:     for (std::size_t i = 0; i < n; ++i)
Chris@102:     {
Chris@102:         const exponential_distribution<RealT,PolicyT> exp(rates[i]);
Chris@102: 
Chris@102:         result += probs[i]*cdf(exp, x);
Chris@102:     }
Chris@102: 
Chris@102:     return result;
Chris@102: }
Chris@102: 
Chris@102: template <typename RealT, typename PolicyT>
Chris@102: RealT quantile(hyperexponential_distribution<RealT, PolicyT> const& dist, RealT const& p)
Chris@102: {
Chris@102:     return hyperexp_detail::quantile_impl(dist, p , false);
Chris@102: }
Chris@102: 
Chris@102: template <typename RealT, typename PolicyT>
Chris@102: RealT cdf(complemented2_type<hyperexponential_distribution<RealT,PolicyT>, RealT> const& c)
Chris@102: {
Chris@102:     RealT const& x = c.param;
Chris@102:     hyperexponential_distribution<RealT,PolicyT> const& dist = c.dist;
Chris@102: 
Chris@102:     RealT result = 0;
Chris@102: 
Chris@102:     if (!hyperexp_detail::check_x("boost::math::cdf(boost::math::complemented2_type<const boost::math::hyperexponential_distribution<%1%>&, %1%>)", x, &result, PolicyT()))
Chris@102:     {
Chris@102:         return result;
Chris@102:     }
Chris@102: 
Chris@102:     const std::size_t n = dist.num_phases();
Chris@102:     const std::vector<RealT> probs = dist.probabilities();
Chris@102:     const std::vector<RealT> rates = dist.rates();
Chris@102: 
Chris@102:     for (std::size_t i = 0; i < n; ++i)
Chris@102:     {
Chris@102:         const exponential_distribution<RealT,PolicyT> exp(rates[i]);
Chris@102: 
Chris@102:         result += probs[i]*cdf(complement(exp, x));
Chris@102:     }
Chris@102: 
Chris@102:     return result;
Chris@102: }
Chris@102: 
Chris@102: 
Chris@102: template <typename RealT, typename PolicyT>
Chris@102: RealT quantile(complemented2_type<hyperexponential_distribution<RealT, PolicyT>, RealT> const& c)
Chris@102: {
Chris@102:     RealT const& p = c.param;
Chris@102:     hyperexponential_distribution<RealT,PolicyT> const& dist = c.dist;
Chris@102: 
Chris@102:     return hyperexp_detail::quantile_impl(dist, p , true);
Chris@102: }
Chris@102: 
Chris@102: template <typename RealT, typename PolicyT>
Chris@102: RealT mean(hyperexponential_distribution<RealT, PolicyT> const& dist)
Chris@102: {
Chris@102:     RealT result = 0;
Chris@102: 
Chris@102:     const std::size_t n = dist.num_phases();
Chris@102:     const std::vector<RealT> probs = dist.probabilities();
Chris@102:     const std::vector<RealT> rates = dist.rates();
Chris@102: 
Chris@102:     for (std::size_t i = 0; i < n; ++i)
Chris@102:     {
Chris@102:         const exponential_distribution<RealT,PolicyT> exp(rates[i]);
Chris@102: 
Chris@102:         result += probs[i]*mean(exp);
Chris@102:     }
Chris@102: 
Chris@102:     return result;
Chris@102: }
Chris@102: 
Chris@102: template <typename RealT, typename PolicyT>
Chris@102: RealT variance(hyperexponential_distribution<RealT, PolicyT> const& dist)
Chris@102: {
Chris@102:     RealT result = 0;
Chris@102: 
Chris@102:     const std::size_t n = dist.num_phases();
Chris@102:     const std::vector<RealT> probs = dist.probabilities();
Chris@102:     const std::vector<RealT> rates = dist.rates();
Chris@102: 
Chris@102:     for (std::size_t i = 0; i < n; ++i)
Chris@102:     {
Chris@102:         result += probs[i]/(rates[i]*rates[i]);
Chris@102:     }
Chris@102: 
Chris@102:     const RealT mean = boost::math::mean(dist);
Chris@102: 
Chris@102:     result = 2*result-mean*mean;
Chris@102: 
Chris@102:     return result;
Chris@102: }
Chris@102: 
Chris@102: template <typename RealT, typename PolicyT>
Chris@102: RealT skewness(hyperexponential_distribution<RealT,PolicyT> const& dist)
Chris@102: {
Chris@102:     BOOST_MATH_STD_USING
Chris@102:     const std::size_t n = dist.num_phases();
Chris@102:     const std::vector<RealT> probs = dist.probabilities();
Chris@102:     const std::vector<RealT> rates = dist.rates();
Chris@102: 
Chris@102:     RealT s1 = 0; // \sum_{i=1}^n \frac{p_i}{\lambda_i}
Chris@102:     RealT s2 = 0; // \sum_{i=1}^n \frac{p_i}{\lambda_i^2}
Chris@102:     RealT s3 = 0; // \sum_{i=1}^n \frac{p_i}{\lambda_i^3}
Chris@102:     for (std::size_t i = 0; i < n; ++i)
Chris@102:     {
Chris@102:         const RealT p = probs[i];
Chris@102:         const RealT r = rates[i];
Chris@102:         const RealT r2 = r*r;
Chris@102:         const RealT r3 = r2*r;
Chris@102: 
Chris@102:         s1 += p/r;
Chris@102:         s2 += p/r2;
Chris@102:         s3 += p/r3;
Chris@102:     }
Chris@102: 
Chris@102:     const RealT s1s1 = s1*s1;
Chris@102: 
Chris@102:     const RealT num = (6*s3 - (3*(2*s2 - s1s1) + s1s1)*s1);
Chris@102:     const RealT den = (2*s2 - s1s1);
Chris@102: 
Chris@102:     return num / pow(den, static_cast<RealT>(1.5));
Chris@102: }
Chris@102: 
Chris@102: template <typename RealT, typename PolicyT>
Chris@102: RealT kurtosis(hyperexponential_distribution<RealT,PolicyT> const& dist)
Chris@102: {
Chris@102:     const std::size_t n = dist.num_phases();
Chris@102:     const std::vector<RealT> probs = dist.probabilities();
Chris@102:     const std::vector<RealT> rates = dist.rates();
Chris@102: 
Chris@102:     RealT s1 = 0; // \sum_{i=1}^n \frac{p_i}{\lambda_i}
Chris@102:     RealT s2 = 0; // \sum_{i=1}^n \frac{p_i}{\lambda_i^2}
Chris@102:     RealT s3 = 0; // \sum_{i=1}^n \frac{p_i}{\lambda_i^3}
Chris@102:     RealT s4 = 0; // \sum_{i=1}^n \frac{p_i}{\lambda_i^4}
Chris@102:     for (std::size_t i = 0; i < n; ++i)
Chris@102:     {
Chris@102:         const RealT p = probs[i];
Chris@102:         const RealT r = rates[i];
Chris@102:         const RealT r2 = r*r;
Chris@102:         const RealT r3 = r2*r;
Chris@102:         const RealT r4 = r3*r;
Chris@102: 
Chris@102:         s1 += p/r;
Chris@102:         s2 += p/r2;
Chris@102:         s3 += p/r3;
Chris@102:         s4 += p/r4;
Chris@102:     }
Chris@102: 
Chris@102:     const RealT s1s1 = s1*s1;
Chris@102: 
Chris@102:     const RealT num = (24*s4 - 24*s3*s1 + 3*(2*(2*s2 - s1s1) + s1s1)*s1s1);
Chris@102:     const RealT den = (2*s2 - s1s1);
Chris@102: 
Chris@102:     return num/(den*den);
Chris@102: }
Chris@102: 
Chris@102: template <typename RealT, typename PolicyT>
Chris@102: RealT kurtosis_excess(hyperexponential_distribution<RealT,PolicyT> const& dist)
Chris@102: {
Chris@102:     return kurtosis(dist) - 3;
Chris@102: }
Chris@102: 
Chris@102: template <typename RealT, typename PolicyT>
Chris@102: RealT mode(hyperexponential_distribution<RealT,PolicyT> const& /*dist*/)
Chris@102: {
Chris@102:     return 0;
Chris@102: }
Chris@102: 
Chris@102: }} // namespace boost::math
Chris@102: 
Chris@102: #ifdef BOOST_MSVC
Chris@102: #pragma warning (pop)
Chris@102: #endif
Chris@102: // This include must be at the end, *after* the accessors
Chris@102: // for this distribution have been defined, in order to
Chris@102: // keep compilers that support two-phase lookup happy.
Chris@102: #include <boost/math/distributions/detail/derived_accessors.hpp>
Chris@102: #include <boost/math/distributions/detail/generic_quantile.hpp>
Chris@102: 
Chris@102: #endif // BOOST_MATH_DISTRIBUTIONS_HYPEREXPONENTIAL