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