642 lines · plain
1// Copyright 2014 Marco Guazzone (marco.guazzone@gmail.com)2//3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6//7// This module implements the Hyper-Exponential distribution.8//9// References:10// - "Queueing Theory in Manufacturing Systems Analysis and Design" by H.T. Papadopolous, C. Heavey and J. Browne (Chapman & Hall/CRC, 1993)11// - http://reference.wolfram.com/language/ref/HyperexponentialDistribution.html12// - http://en.wikipedia.org/wiki/Hyperexponential_distribution13//14 15#ifndef BOOST_MATH_DISTRIBUTIONS_HYPEREXPONENTIAL_HPP16#define BOOST_MATH_DISTRIBUTIONS_HYPEREXPONENTIAL_HPP17 18#include <boost/math/tools/cxx03_warn.hpp>19#include <boost/math/distributions/complement.hpp>20#include <boost/math/distributions/detail/common_error_handling.hpp>21#include <boost/math/distributions/exponential.hpp>22#include <boost/math/policies/policy.hpp>23#include <boost/math/special_functions/fpclassify.hpp>24#include <boost/math/tools/precision.hpp>25#include <boost/math/tools/roots.hpp>26#include <boost/math/tools/is_detected.hpp>27#include <cstddef>28#include <iterator>29#include <limits>30#include <numeric>31#include <utility>32#include <vector>33#include <type_traits>34#include <initializer_list>35 36 37#ifdef _MSC_VER38# pragma warning (push)39# pragma warning(disable:4127) // conditional expression is constant40# pragma warning(disable:4389) // '==' : signed/unsigned mismatch in test_tools41#endif // _MSC_VER42 43namespace boost { namespace math {44 45namespace detail {46 47template <typename Dist>48typename 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);49 50} // Namespace detail51 52 53template <typename RealT, typename PolicyT>54class hyperexponential_distribution;55 56 57namespace /*<unnamed>*/ { namespace hyperexp_detail {58 59template <typename T>60void normalize(std::vector<T>& v)61{62 if(!v.size())63 return; // Our error handlers will get this later64 const T sum = std::accumulate(v.begin(), v.end(), static_cast<T>(0));65 T final_sum = 0;66 const typename std::vector<T>::iterator end = --v.end();67 for (typename std::vector<T>::iterator it = v.begin();68 it != end;69 ++it)70 {71 *it /= sum;72 final_sum += *it;73 }74 *end = 1 - final_sum; // avoids round off errors, ensures the probs really do sum to 1.75}76 77template <typename RealT, typename PolicyT>78bool check_probabilities(char const* function, std::vector<RealT> const& probabilities, RealT* presult, PolicyT const& pol)79{80 BOOST_MATH_STD_USING81 const std::size_t n = probabilities.size();82 RealT sum = 0;83 for (std::size_t i = 0; i < n; ++i)84 {85 if (probabilities[i] < 086 || probabilities[i] > 187 || !(boost::math::isfinite)(probabilities[i]))88 {89 *presult = policies::raise_domain_error<RealT>(function,90 "The elements of parameter \"probabilities\" must be >= 0 and <= 1, but at least one of them was: %1%.",91 probabilities[i],92 pol);93 return false;94 }95 sum += probabilities[i];96 }97 98 //99 // We try to keep phase probabilities correctly normalized in the distribution constructors,100 // however in practice we have to allow for a very slight divergence from a sum of exactly 1:101 //102 if (fabs(sum - 1) > tools::epsilon<RealT>() * 2)103 {104 *presult = policies::raise_domain_error<RealT>(function,105 "The elements of parameter \"probabilities\" must sum to 1, but their sum is: %1%.",106 sum,107 pol);108 return false;109 }110 111 return true;112}113 114template <typename RealT, typename PolicyT>115bool check_rates(char const* function, std::vector<RealT> const& rates, RealT* presult, PolicyT const& pol)116{117 const std::size_t n = rates.size();118 for (std::size_t i = 0; i < n; ++i)119 {120 if (rates[i] <= 0121 || !(boost::math::isfinite)(rates[i]))122 {123 *presult = policies::raise_domain_error<RealT>(function,124 "The elements of parameter \"rates\" must be > 0, but at least one of them is: %1%.",125 rates[i],126 pol);127 return false;128 }129 }130 return true;131}132 133template <typename RealT, typename PolicyT>134bool check_dist(char const* function, std::vector<RealT> const& probabilities, std::vector<RealT> const& rates, RealT* presult, PolicyT const& pol)135{136 BOOST_MATH_STD_USING137 if (probabilities.size() != rates.size())138 {139 *presult = policies::raise_domain_error<RealT>(function,140 R"(The parameters "probabilities" and "rates" must have the same length, but their size differ by: %1%.)",141 fabs(static_cast<RealT>(probabilities.size())-static_cast<RealT>(rates.size())),142 pol);143 return false;144 }145 146 return check_probabilities(function, probabilities, presult, pol)147 && check_rates(function, rates, presult, pol);148}149 150template <typename RealT, typename PolicyT>151bool check_x(char const* function, RealT x, RealT* presult, PolicyT const& pol)152{153 if (x < 0 || (boost::math::isnan)(x))154 {155 *presult = policies::raise_domain_error<RealT>(function, "The random variable must be >= 0, but is: %1%.", x, pol);156 return false;157 }158 return true;159}160 161template <typename RealT, typename PolicyT>162bool check_probability(char const* function, RealT p, RealT* presult, PolicyT const& pol)163{164 if (p < 0 || p > 1 || (boost::math::isnan)(p))165 {166 *presult = policies::raise_domain_error<RealT>(function, "The probability be >= 0 and <= 1, but is: %1%.", p, pol);167 return false;168 }169 return true;170}171 172template <typename RealT, typename PolicyT>173RealT quantile_impl(hyperexponential_distribution<RealT, PolicyT> const& dist, RealT const& p, bool comp)174{175 // Don't have a closed form so try to numerically solve the inverse CDF...176 177 typedef typename policies::evaluation<RealT, PolicyT>::type value_type;178 typedef typename policies::normalise<PolicyT,179 policies::promote_float<false>,180 policies::promote_double<false>,181 policies::discrete_quantile<>,182 policies::assert_undefined<> >::type forwarding_policy;183 184 static const char* function = comp ? "boost::math::quantile(const boost::math::complemented2_type<boost::math::hyperexponential_distribution<%1%>, %1%>&)"185 : "boost::math::quantile(const boost::math::hyperexponential_distribution<%1%>&, %1%)";186 187 RealT result = 0;188 189 if (!check_probability(function, p, &result, PolicyT()))190 {191 return result;192 }193 194 const std::size_t n = dist.num_phases();195 const std::vector<RealT> probs = dist.probabilities();196 const std::vector<RealT> rates = dist.rates();197 198 // A possible (but inaccurate) approximation is given below, where the199 // quantile is given by the weighted sum of exponential quantiles:200 RealT guess = 0;201 if (comp)202 {203 for (std::size_t i = 0; i < n; ++i)204 {205 const exponential_distribution<RealT,PolicyT> exp(rates[i]);206 207 guess += probs[i]*quantile(complement(exp, p));208 }209 }210 else211 {212 for (std::size_t i = 0; i < n; ++i)213 {214 const exponential_distribution<RealT,PolicyT> exp(rates[i]);215 216 guess += probs[i]*quantile(exp, p);217 }218 }219 220 // Fast return in case the Hyper-Exponential is essentially an Exponential221 if (n == 1)222 {223 return guess;224 }225 226 value_type q;227 q = detail::generic_quantile(hyperexponential_distribution<RealT,forwarding_policy>(probs, rates),228 p,229 guess,230 comp,231 function);232 233 result = policies::checked_narrowing_cast<RealT,forwarding_policy>(q, function);234 235 return result;236}237 238}} // Namespace <unnamed>::hyperexp_detail239 240 241template <typename RealT = double, typename PolicyT = policies::policy<> >242class hyperexponential_distribution243{244 public: typedef RealT value_type;245 public: typedef PolicyT policy_type;246 247 248 public: hyperexponential_distribution()249 : probs_(1, 1),250 rates_(1, 1)251 {252 RealT err;253 hyperexp_detail::check_dist("boost::math::hyperexponential_distribution<%1%>::hyperexponential_distribution",254 probs_,255 rates_,256 &err,257 PolicyT());258 }259 260 // Four arg constructor: no ambiguity here, the arguments must be two pairs of iterators:261 public: template <typename ProbIterT, typename RateIterT>262 hyperexponential_distribution(ProbIterT prob_first, ProbIterT prob_last,263 RateIterT rate_first, RateIterT rate_last)264 : probs_(prob_first, prob_last),265 rates_(rate_first, rate_last)266 {267 hyperexp_detail::normalize(probs_);268 RealT err;269 hyperexp_detail::check_dist("boost::math::hyperexponential_distribution<%1%>::hyperexponential_distribution",270 probs_,271 rates_,272 &err,273 PolicyT());274 }275 private: template <typename T, typename = void>276 struct is_iterator277 {278 static constexpr bool value = false;279 };280 281 template <typename T>282 struct is_iterator<T, boost::math::tools::void_t<typename std::iterator_traits<T>::difference_type>>283 {284 // std::iterator_traits<T>::difference_type returns void for invalid types285 static constexpr bool value = !std::is_same<typename std::iterator_traits<T>::difference_type, void>::value;286 };287 288 // Two arg constructor from 2 ranges, we SFINAE this out of existence if289 // either argument type is incrementable as in that case the type is290 // probably an iterator:291 public: template <typename ProbRangeT, typename RateRangeT, 292 typename std::enable_if<!is_iterator<ProbRangeT>::value && 293 !is_iterator<RateRangeT>::value, bool>::type = true>294 hyperexponential_distribution(ProbRangeT const& prob_range,295 RateRangeT const& rate_range)296 : probs_(std::begin(prob_range), std::end(prob_range)),297 rates_(std::begin(rate_range), std::end(rate_range))298 {299 hyperexp_detail::normalize(probs_);300 301 RealT err;302 hyperexp_detail::check_dist("boost::math::hyperexponential_distribution<%1%>::hyperexponential_distribution",303 probs_,304 rates_,305 &err,306 PolicyT());307 }308 309 // Two arg constructor for a pair of iterators: we SFINAE this out of310 // existence if neither argument types are incrementable.311 // Note that we allow different argument types here to allow for312 // construction from an array plus a pointer into that array.313 public: template <typename RateIterT, typename RateIterT2, 314 typename std::enable_if<is_iterator<RateIterT>::value || 315 is_iterator<RateIterT2>::value, bool>::type = true>316 hyperexponential_distribution(RateIterT const& rate_first, 317 RateIterT2 const& rate_last)318 : probs_(std::distance(rate_first, rate_last), 1), // will be normalized below319 rates_(rate_first, rate_last)320 {321 hyperexp_detail::normalize(probs_);322 323 RealT err;324 hyperexp_detail::check_dist("boost::math::hyperexponential_distribution<%1%>::hyperexponential_distribution",325 probs_,326 rates_,327 &err,328 PolicyT());329 }330 331 // Initializer list constructor: allows for construction from array literals:332public: hyperexponential_distribution(std::initializer_list<RealT> l1, std::initializer_list<RealT> l2)333 : probs_(l1.begin(), l1.end()),334 rates_(l2.begin(), l2.end())335 {336 hyperexp_detail::normalize(probs_);337 338 RealT err;339 hyperexp_detail::check_dist("boost::math::hyperexponential_distribution<%1%>::hyperexponential_distribution",340 probs_,341 rates_,342 &err,343 PolicyT());344 }345 346public: hyperexponential_distribution(std::initializer_list<RealT> l1)347 : probs_(l1.size(), 1),348 rates_(l1.begin(), l1.end())349 {350 hyperexp_detail::normalize(probs_);351 352 RealT err;353 hyperexp_detail::check_dist("boost::math::hyperexponential_distribution<%1%>::hyperexponential_distribution",354 probs_,355 rates_,356 &err,357 PolicyT());358 }359 360 // Single argument constructor: argument must be a range.361 public: template <typename RateRangeT>362 hyperexponential_distribution(RateRangeT const& rate_range)363 : probs_(std::distance(std::begin(rate_range), std::end(rate_range)), 1), // will be normalized below364 rates_(std::begin(rate_range), std::end(rate_range))365 {366 hyperexp_detail::normalize(probs_);367 368 RealT err;369 hyperexp_detail::check_dist("boost::math::hyperexponential_distribution<%1%>::hyperexponential_distribution",370 probs_,371 rates_,372 &err,373 PolicyT());374 }375 376 public: std::vector<RealT> probabilities() const377 {378 return probs_;379 }380 381 public: std::vector<RealT> rates() const382 {383 return rates_;384 }385 386 public: std::size_t num_phases() const387 {388 return rates_.size();389 }390 391 392 private: std::vector<RealT> probs_;393 private: std::vector<RealT> rates_;394}; // class hyperexponential_distribution395 396 397// Convenient type synonym for double.398typedef hyperexponential_distribution<double> hyperexponential;399 400 401// Range of permissible values for random variable x402template <typename RealT, typename PolicyT>403std::pair<RealT,RealT> range(hyperexponential_distribution<RealT,PolicyT> const&)404{405 if (std::numeric_limits<RealT>::has_infinity)406 {407 return std::make_pair(static_cast<RealT>(0), std::numeric_limits<RealT>::infinity()); // 0 to +inf.408 }409 410 return std::make_pair(static_cast<RealT>(0), tools::max_value<RealT>()); // 0 to +<max value>411}412 413// Range of supported values for random variable x.414// This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.415template <typename RealT, typename PolicyT>416std::pair<RealT,RealT> support(hyperexponential_distribution<RealT,PolicyT> const&)417{418 return std::make_pair(tools::min_value<RealT>(), tools::max_value<RealT>()); // <min value> to +<max value>.419}420 421template <typename RealT, typename PolicyT>422RealT pdf(hyperexponential_distribution<RealT, PolicyT> const& dist, RealT const& x)423{424 BOOST_MATH_STD_USING425 RealT result = 0;426 427 if (!hyperexp_detail::check_x("boost::math::pdf(const boost::math::hyperexponential_distribution<%1%>&, %1%)", x, &result, PolicyT()))428 {429 return result;430 }431 432 const std::size_t n = dist.num_phases();433 const std::vector<RealT> probs = dist.probabilities();434 const std::vector<RealT> rates = dist.rates();435 436 for (std::size_t i = 0; i < n; ++i)437 {438 const exponential_distribution<RealT,PolicyT> exp(rates[i]);439 440 result += probs[i]*pdf(exp, x);441 //result += probs[i]*rates[i]*exp(-rates[i]*x);442 }443 444 return result;445}446 447template <typename RealT, typename PolicyT>448RealT cdf(hyperexponential_distribution<RealT, PolicyT> const& dist, RealT const& x)449{450 RealT result = 0;451 452 if (!hyperexp_detail::check_x("boost::math::cdf(const boost::math::hyperexponential_distribution<%1%>&, %1%)", x, &result, PolicyT()))453 {454 return result;455 }456 457 const std::size_t n = dist.num_phases();458 const std::vector<RealT> probs = dist.probabilities();459 const std::vector<RealT> rates = dist.rates();460 461 for (std::size_t i = 0; i < n; ++i)462 {463 const exponential_distribution<RealT,PolicyT> exp(rates[i]);464 465 result += probs[i]*cdf(exp, x);466 }467 468 return result;469}470 471template <typename RealT, typename PolicyT>472RealT quantile(hyperexponential_distribution<RealT, PolicyT> const& dist, RealT const& p)473{474 return hyperexp_detail::quantile_impl(dist, p , false);475}476 477template <typename RealT, typename PolicyT>478RealT cdf(complemented2_type<hyperexponential_distribution<RealT,PolicyT>, RealT> const& c)479{480 RealT const& x = c.param;481 hyperexponential_distribution<RealT,PolicyT> const& dist = c.dist;482 483 RealT result = 0;484 485 if (!hyperexp_detail::check_x("boost::math::cdf(boost::math::complemented2_type<const boost::math::hyperexponential_distribution<%1%>&, %1%>)", x, &result, PolicyT()))486 {487 return result;488 }489 490 const std::size_t n = dist.num_phases();491 const std::vector<RealT> probs = dist.probabilities();492 const std::vector<RealT> rates = dist.rates();493 494 for (std::size_t i = 0; i < n; ++i)495 {496 const exponential_distribution<RealT,PolicyT> exp(rates[i]);497 498 result += probs[i]*cdf(complement(exp, x));499 }500 501 return result;502}503 504 505template <typename RealT, typename PolicyT>506RealT quantile(complemented2_type<hyperexponential_distribution<RealT, PolicyT>, RealT> const& c)507{508 RealT const& p = c.param;509 hyperexponential_distribution<RealT,PolicyT> const& dist = c.dist;510 511 return hyperexp_detail::quantile_impl(dist, p , true);512}513 514template <typename RealT, typename PolicyT>515RealT mean(hyperexponential_distribution<RealT, PolicyT> const& dist)516{517 RealT result = 0;518 519 const std::size_t n = dist.num_phases();520 const std::vector<RealT> probs = dist.probabilities();521 const std::vector<RealT> rates = dist.rates();522 523 for (std::size_t i = 0; i < n; ++i)524 {525 const exponential_distribution<RealT,PolicyT> exp(rates[i]);526 527 result += probs[i]*mean(exp);528 }529 530 return result;531}532 533template <typename RealT, typename PolicyT>534RealT variance(hyperexponential_distribution<RealT, PolicyT> const& dist)535{536 RealT result = 0;537 538 const std::size_t n = dist.num_phases();539 const std::vector<RealT> probs = dist.probabilities();540 const std::vector<RealT> rates = dist.rates();541 542 for (std::size_t i = 0; i < n; ++i)543 {544 result += probs[i]/(rates[i]*rates[i]);545 }546 547 const RealT mean = boost::math::mean(dist);548 549 result = 2*result-mean*mean;550 551 return result;552}553 554template <typename RealT, typename PolicyT>555RealT skewness(hyperexponential_distribution<RealT,PolicyT> const& dist)556{557 BOOST_MATH_STD_USING558 const std::size_t n = dist.num_phases();559 const std::vector<RealT> probs = dist.probabilities();560 const std::vector<RealT> rates = dist.rates();561 562 RealT s1 = 0; // \sum_{i=1}^n \frac{p_i}{\lambda_i}563 RealT s2 = 0; // \sum_{i=1}^n \frac{p_i}{\lambda_i^2}564 RealT s3 = 0; // \sum_{i=1}^n \frac{p_i}{\lambda_i^3}565 for (std::size_t i = 0; i < n; ++i)566 {567 const RealT p = probs[i];568 const RealT r = rates[i];569 const RealT r2 = r*r;570 const RealT r3 = r2*r;571 572 s1 += p/r;573 s2 += p/r2;574 s3 += p/r3;575 }576 577 const RealT s1s1 = s1*s1;578 579 const RealT num = (6*s3 - (3*(2*s2 - s1s1) + s1s1)*s1);580 const RealT den = (2*s2 - s1s1);581 582 return num / pow(den, static_cast<RealT>(1.5));583}584 585template <typename RealT, typename PolicyT>586RealT kurtosis(hyperexponential_distribution<RealT,PolicyT> const& dist)587{588 const std::size_t n = dist.num_phases();589 const std::vector<RealT> probs = dist.probabilities();590 const std::vector<RealT> rates = dist.rates();591 592 RealT s1 = 0; // \sum_{i=1}^n \frac{p_i}{\lambda_i}593 RealT s2 = 0; // \sum_{i=1}^n \frac{p_i}{\lambda_i^2}594 RealT s3 = 0; // \sum_{i=1}^n \frac{p_i}{\lambda_i^3}595 RealT s4 = 0; // \sum_{i=1}^n \frac{p_i}{\lambda_i^4}596 for (std::size_t i = 0; i < n; ++i)597 {598 const RealT p = probs[i];599 const RealT r = rates[i];600 const RealT r2 = r*r;601 const RealT r3 = r2*r;602 const RealT r4 = r3*r;603 604 s1 += p/r;605 s2 += p/r2;606 s3 += p/r3;607 s4 += p/r4;608 }609 610 const RealT s1s1 = s1*s1;611 612 const RealT num = (24*s4 - 24*s3*s1 + 3*(2*(2*s2 - s1s1) + s1s1)*s1s1);613 const RealT den = (2*s2 - s1s1);614 615 return num/(den*den);616}617 618template <typename RealT, typename PolicyT>619RealT kurtosis_excess(hyperexponential_distribution<RealT,PolicyT> const& dist)620{621 return kurtosis(dist) - 3;622}623 624template <typename RealT, typename PolicyT>625RealT mode(hyperexponential_distribution<RealT,PolicyT> const& /*dist*/)626{627 return 0;628}629 630}} // namespace boost::math631 632#ifdef _MSC_VER633#pragma warning (pop)634#endif635// This include must be at the end, *after* the accessors636// for this distribution have been defined, in order to637// keep compilers that support two-phase lookup happy.638#include <boost/math/distributions/detail/derived_accessors.hpp>639#include <boost/math/distributions/detail/generic_quantile.hpp>640 641#endif // BOOST_MATH_DISTRIBUTIONS_HYPEREXPONENTIAL642