1 /* boost random/hyperexponential_distribution.hpp header file
3 * Copyright Marco Guazzone 2014
4 * Distributed under the Boost Software License, Version 1.0. (See
5 * accompanying file LICENSE_1_0.txt or copy at
6 * http://www.boost.org/LICENSE_1_0.txt)
8 * See http://www.boost.org for most recent version including documentation.
10 * Much of the code here taken by boost::math::hyperexponential_distribution.
11 * To this end, we would like to thank Paul Bristow and John Maddock for their
14 * \author Marco Guazzone (marco.guazzone@gmail.com)
17 #ifndef BOOST_RANDOM_HYPEREXPONENTIAL_DISTRIBUTION_HPP
18 #define BOOST_RANDOM_HYPEREXPONENTIAL_DISTRIBUTION_HPP
21 #include <boost/config.hpp>
22 #include <boost/math/special_functions/fpclassify.hpp>
23 #include <boost/random/detail/operators.hpp>
24 #include <boost/random/detail/vector_io.hpp>
25 #include <boost/random/discrete_distribution.hpp>
26 #include <boost/random/exponential_distribution.hpp>
27 #include <boost/range/begin.hpp>
28 #include <boost/range/end.hpp>
29 #include <boost/range/size.hpp>
30 #include <boost/type_traits/has_pre_increment.hpp>
35 #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
36 # include <initializer_list>
37 #endif // BOOST_NO_CXX11_HDR_INITIALIZER_LIST
44 namespace boost { namespace random {
46 namespace hyperexp_detail {
49 std::vector<T>& normalize(std::vector<T>& v)
56 const T sum = std::accumulate(v.begin(), v.end(), static_cast<T>(0));
59 const typename std::vector<T>::iterator end = --v.end();
60 for (typename std::vector<T>::iterator it = v.begin();
67 *end = 1-final_sum; // avoids round off errors thus ensuring the probabilities really sum to 1
72 template <typename RealT>
73 bool check_probabilities(std::vector<RealT> const& probabilities)
75 const std::size_t n = probabilities.size();
77 for (std::size_t i = 0; i < n; ++i)
79 if (probabilities[i] < 0
80 || probabilities[i] > 1
81 || !(boost::math::isfinite)(probabilities[i]))
85 sum += probabilities[i];
88 //NOTE: the check below seems to fail on some architectures.
89 // So we commented it.
90 //// - We try to keep phase probabilities correctly normalized in the distribution constructors
91 //// - However in practice we have to allow for a very slight divergence from a sum of exactly 1:
92 ////if (std::abs(sum-1) > (std::numeric_limits<RealT>::epsilon()*2))
93 //// This is from Knuth "The Art of Computer Programming: Vol.2, 3rd Ed", and can be used to
94 //// check is two numbers are approximately equal
95 //const RealT one = 1;
96 //const RealT tol = std::numeric_limits<RealT>::epsilon()*2.0;
97 //if (std::abs(sum-one) > (std::max(std::abs(sum), std::abs(one))*tol))
105 template <typename RealT>
106 bool check_rates(std::vector<RealT> const& rates)
108 const std::size_t n = rates.size();
109 for (std::size_t i = 0; i < n; ++i)
112 || !(boost::math::isfinite)(rates[i]))
120 template <typename RealT>
121 bool check_params(std::vector<RealT> const& probabilities, std::vector<RealT> const& rates)
123 if (probabilities.size() != rates.size())
128 return check_probabilities(probabilities)
129 && check_rates(rates);
132 } // Namespace hyperexp_detail
136 * The hyperexponential distribution is a real-valued continuous distribution
137 * with two parameters, the <em>phase probability vector</em> \c probs and the
138 * <em>rate vector</em> \c rates.
140 * A \f$k\f$-phase hyperexponential distribution is a mixture of \f$k\f$
141 * exponential distributions.
142 * For this reason, it is also referred to as <em>mixed exponential
143 * distribution</em> or <em>parallel \f$k\f$-phase exponential
146 * A \f$k\f$-phase hyperexponential distribution is characterized by two
147 * parameters, namely a <em>phase probability vector</em> \f$\mathbf{\alpha}=(\alpha_1,\ldots,\alpha_k)\f$ and a <em>rate vector</em> \f$\mathbf{\lambda}=(\lambda_1,\ldots,\lambda_k)\f$.
149 * A \f$k\f$-phase hyperexponential distribution is frequently used in
150 * <em>queueing theory</em> to model the distribution of the superposition of
151 * \f$k\f$ independent events, like, for instance, the service time distribution
152 * of a queueing station with \f$k\f$ servers in parallel where the \f$i\f$-th
153 * server is chosen with probability \f$\alpha_i\f$ and its service time
154 * distribution is an exponential distribution with rate \f$\lambda_i\f$
155 * (Allen,1990; Papadopolous et al.,1993; Trivedi,2002).
157 * For instance, CPUs service-time distribution in a computing system has often
158 * been observed to possess such a distribution (Rosin,1965).
159 * Also, the arrival of different types of customer to a single queueing station
160 * is often modeled as a hyperexponential distribution (Papadopolous et al.,1993).
161 * Similarly, if a product manufactured in several parallel assemply lines and
162 * the outputs are merged, the failure density of the overall product is likely
163 * to be hyperexponential (Trivedi,2002).
165 * Finally, since the hyperexponential distribution exhibits a high Coefficient
166 * of Variation (CoV), that is a CoV > 1, it is especially suited to fit
167 * empirical data with large CoV (Feitelson,2014; Wolski et al.,2013) and to
168 * approximate <em>long-tail probability distributions</em> (Feldmann et al.,1998).
170 * See (Boost,2014) for more information and examples.
172 * A \f$k\f$-phase hyperexponential distribution has a probability density
175 * f(x) = \sum_{i=1}^k \alpha_i \lambda_i e^{-x\lambda_i}
178 * - \f$k\f$ is the <em>number of phases</em> and also the size of the input
180 * - \f$\mathbf{\alpha}=(\alpha_1,\ldots,\alpha_k)\f$ is the <em>phase probability
181 * vector</em> parameter, and
182 * - \f$\mathbf{\lambda}=(\lambda_1,\ldots,\lambda_k)\f$ is the <em>rate vector</em>
186 * Given a \f$k\f$-phase hyperexponential distribution with phase probability
187 * vector \f$\mathbf{\alpha}\f$ and rate vector \f$\mathbf{\lambda}\f$, the
188 * random variate generation algorithm consists of the following steps (Tyszer,1999):
189 * -# Generate a random variable \f$U\f$ uniformly distribution on the interval \f$(0,1)\f$.
190 * -# Use \f$U\f$ to select the appropriate \f$\lambda_i\f$ (e.g., the
191 * <em>alias method</em> can possibly be used for this step).
192 * -# Generate an exponentially distributed random variable \f$X\f$ with rate parameter \f$\lambda_i\f$.
197 * -# A.O. Allen, <em>Probability, Statistics, and Queuing Theory with Computer Science Applications, Second Edition</em>, Academic Press, 1990.
198 * -# Boost C++ Libraries, <em>Boost.Math / Statistical Distributions: Hyperexponential Distribution</em>, Online: http://www.boost.org/doc/libs/release/libs/math/doc/html/dist.html , 2014.
199 * -# D.G. Feitelson, <em>Workload Modeling for Computer Systems Performance Evaluation</em>, Cambridge University Press, 2014
200 * -# A. Feldmann and W. Whitt, <em>Fitting mixtures of exponentials to long-tail distributions to analyze network performance models</em>, Performance Evaluation 31(3-4):245, doi:10.1016/S0166-5316(97)00003-5, 1998.
201 * -# H.T. Papadopolous, C. Heavey and J. Browne, <em>Queueing Theory in Manufacturing Systems Analysis and Design</em>, Chapman & Hall/CRC, 1993, p. 35.
202 * -# R.F. Rosin, <em>Determining a computing center environment</em>, Communications of the ACM 8(7):463-468, 1965.
203 * -# K.S. Trivedi, <em>Probability and Statistics with Reliability, Queueing, and Computer Science Applications</em>, John Wiley & Sons, Inc., 2002.
204 * -# J. Tyszer, <em>Object-Oriented Computer Simulation of Discrete-Event Systems</em>, Springer, 1999.
205 * -# Wikipedia, <em>Hyperexponential Distribution</em>, Online: http://en.wikipedia.org/wiki/Hyperexponential_distribution , 2014.
206 * -# Wolfram Mathematica, <em>Hyperexponential Distribution</em>, Online: http://reference.wolfram.com/language/ref/HyperexponentialDistribution.html , 2014.
209 * \author Marco Guazzone (marco.guazzone@gmail.com)
211 template<class RealT = double>
212 class hyperexponential_distribution
214 public: typedef RealT result_type;
215 public: typedef RealT input_type;
219 * The parameters of a hyperexponential distribution.
221 * Stores the <em>phase probability vector</em> and the <em>rate vector</em>
222 * of the hyperexponential distribution.
224 * \author Marco Guazzone (marco.guazzone@gmail.com)
226 public: class param_type
228 public: typedef hyperexponential_distribution distribution_type;
231 * Constructs a \c param_type with the default parameters
232 * of the distribution.
241 * Constructs a \c param_type from the <em>phase probability vector</em>
242 * and <em>rate vector</em> parameters of the distribution.
244 * The <em>phase probability vector</em> parameter is given by the range
245 * defined by [\a prob_first, \a prob_last) iterator pair, and the
246 * <em>rate vector</em> parameter is given by the range defined by
247 * [\a rate_first, \a rate_last) iterator pair.
249 * \tparam ProbIterT Must meet the requirements of \c InputIterator concept (ISO,2014,sec. 24.2.3 [input.iterators]).
250 * \tparam RateIterT Must meet the requirements of \c InputIterator concept (ISO,2014,sec. 24.2.3 [input.iterators]).
252 * \param prob_first The iterator to the beginning of the range of non-negative real elements representing the phase probabilities; if elements don't sum to 1, they are normalized.
253 * \param prob_last The iterator to the ending of the range of non-negative real elements representing the phase probabilities; if elements don't sum to 1, they are normalized.
254 * \param rate_first The iterator to the beginning of the range of non-negative real elements representing the rates.
255 * \param rate_last The iterator to the ending of the range of non-negative real elements representing the rates.
258 * -# ISO, <em>ISO/IEC 14882-2014: Information technology - Programming languages - C++</em>, 2014
261 public: template <typename ProbIterT, typename RateIterT>
262 param_type(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)
267 hyperexp_detail::normalize(probs_);
269 assert( hyperexp_detail::check_params(probs_, rates_) );
273 * Constructs a \c param_type from the <em>phase probability vector</em>
274 * and <em>rate vector</em> parameters of the distribution.
276 * The <em>phase probability vector</em> parameter is given by the range
277 * defined by \a prob_range, and the <em>rate vector</em> parameter is
278 * given by the range defined by \a rate_range.
280 * \tparam ProbRangeT Must meet the requirements of <a href="boost:/libs/range/doc/html/range/concepts.html">Range</a> concept.
281 * \tparam RateRangeT Must meet the requirements of <a href="boost:/libs/range/doc/html/range/concepts.html">Range</a> concept.
283 * \param prob_range The range of non-negative real elements representing the phase probabilities; if elements don't sum to 1, they are normalized.
284 * \param rate_range The range of positive real elements representing the rates.
287 * The final \c disable_if parameter is an implementation detail that
288 * differentiates between this two argument constructor and the
289 * iterator-based two argument constructor described below.
291 // We SFINAE this out of existance if either argument type is
292 // incrementable as in that case the type is probably an iterator:
293 public: template <typename ProbRangeT, typename RateRangeT>
294 param_type(ProbRangeT const& prob_range,
295 RateRangeT const& rate_range,
296 typename boost::disable_if_c<boost::has_pre_increment<ProbRangeT>::value || boost::has_pre_increment<RateRangeT>::value>::type* = 0)
297 : probs_(boost::begin(prob_range), boost::end(prob_range)),
298 rates_(boost::begin(rate_range), boost::end(rate_range))
300 hyperexp_detail::normalize(probs_);
302 assert( hyperexp_detail::check_params(probs_, rates_) );
306 * Constructs a \c param_type from the <em>rate vector</em> parameter of
307 * the distribution and with equal phase probabilities.
309 * The <em>rate vector</em> parameter is given by the range defined by
310 * [\a rate_first, \a rate_last) iterator pair, and the <em>phase
311 * probability vector</em> parameter is set to the equal phase
312 * probabilities (i.e., to a vector of the same length \f$k\f$ of the
313 * <em>rate vector</em> and with each element set to \f$1.0/k\f$).
315 * \tparam RateIterT Must meet the requirements of \c InputIterator concept (ISO,2014,sec. 24.2.3 [input.iterators]).
316 * \tparam RateIterT2 Must meet the requirements of \c InputIterator concept (ISO,2014,sec. 24.2.3 [input.iterators]).
318 * \param rate_first The iterator to the beginning of the range of non-negative real elements representing the rates.
319 * \param rate_last The iterator to the ending of the range of non-negative real elements representing the rates.
322 * The final \c disable_if parameter is an implementation detail that
323 * differentiates between this two argument constructor and the
324 * range-based two argument constructor described above.
327 * -# ISO, <em>ISO/IEC 14882-2014: Information technology - Programming languages - C++</em>, 2014
330 // We SFINAE this out of existance if the argument type is
331 // incrementable as in that case the type is probably an iterator.
332 public: template <typename RateIterT>
333 param_type(RateIterT rate_first,
335 typename boost::enable_if_c<boost::has_pre_increment<RateIterT>::value>::type* = 0)
336 : probs_(std::distance(rate_first, rate_last), 1), // will be normalized below
337 rates_(rate_first, rate_last)
339 assert(probs_.size() == rates_.size());
343 * Constructs a @c param_type from the "rates" parameters
344 * of the distribution and with equal phase probabilities.
346 * The <em>rate vector</em> parameter is given by the range defined by
347 * \a rate_range, and the <em>phase probability vector</em> parameter is
348 * set to the equal phase probabilities (i.e., to a vector of the same
349 * length \f$k\f$ of the <em>rate vector</em> and with each element set
352 * \tparam RateRangeT Must meet the requirements of <a href="boost:/libs/range/doc/html/range/concepts.html">Range</a> concept.
354 * \param rate_range The range of positive real elements representing the rates.
356 public: template <typename RateRangeT>
357 param_type(RateRangeT const& rate_range)
358 : probs_(boost::size(rate_range), 1), // Will be normalized below
359 rates_(boost::begin(rate_range), boost::end(rate_range))
361 hyperexp_detail::normalize(probs_);
363 assert( hyperexp_detail::check_params(probs_, rates_) );
366 #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
368 * Constructs a \c param_type from the <em>phase probability vector</em>
369 * and <em>rate vector</em> parameters of the distribution.
371 * The <em>phase probability vector</em> parameter is given by the
372 * <em>brace-init-list</em> (ISO,2014,sec. 8.5.4 [dcl.init.list])
373 * defined by \a l1, and the <em>rate vector</em> parameter is given by the
374 * <em>brace-init-list</em> (ISO,2014,sec. 8.5.4 [dcl.init.list])
377 * \param l1 The initializer list for inizializing the phase probability vector.
378 * \param l2 The initializer list for inizializing the rate vector.
381 * -# ISO, <em>ISO/IEC 14882-2014: Information technology - Programming languages - C++</em>, 2014
384 public: param_type(std::initializer_list<RealT> l1, std::initializer_list<RealT> l2)
385 : probs_(l1.begin(), l1.end()),
386 rates_(l2.begin(), l2.end())
388 hyperexp_detail::normalize(probs_);
390 assert( hyperexp_detail::check_params(probs_, rates_) );
394 * Constructs a \c param_type from the <em>rate vector</em> parameter
395 * of the distribution and with equal phase probabilities.
397 * The <em>rate vector</em> parameter is given by the
398 * <em>brace-init-list</em> (ISO,2014,sec. 8.5.4 [dcl.init.list])
399 * defined by \a l1, and the <em>phase probability vector</em> parameter is
400 * set to the equal phase probabilities (i.e., to a vector of the same
401 * length \f$k\f$ of the <em>rate vector</em> and with each element set
404 * \param l1 The initializer list for inizializing the rate vector.
407 * -# ISO, <em>ISO/IEC 14882-2014: Information technology - Programming languages - C++</em>, 2014
410 public: param_type(std::initializer_list<RealT> l1)
411 : probs_(std::distance(l1.begin(), l1.end()), 1), // Will be normalized below
412 rates_(l1.begin(), l1.end())
414 hyperexp_detail::normalize(probs_);
416 assert( hyperexp_detail::check_params(probs_, rates_) );
418 #endif // BOOST_NO_CXX11_HDR_INITIALIZER_LIST
421 * Gets the <em>phase probability vector</em> parameter of the distribtuion.
423 * \return The <em>phase probability vector</em> parameter of the distribution.
426 * The returned probabilities are the normalized version of the ones
427 * passed at construction time.
429 public: std::vector<RealT> probabilities() const
435 * Gets the <em>rate vector</em> parameter of the distribtuion.
437 * \return The <em>rate vector</em> parameter of the distribution.
439 public: std::vector<RealT> rates() const
444 /** Writes a \c param_type to a \c std::ostream. */
445 public: BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, param)
447 detail::print_vector(os, param.probs_);
449 detail::print_vector(os, param.rates_);
454 /** Reads a \c param_type from a \c std::istream. */
455 public: BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, param)
457 // NOTE: if \c std::ios_base::exceptions is set, the code below may
458 // throw in case of a I/O failure.
459 // To prevent leaving the state of \c param inconsistent:
460 // - if an exception is thrown, the state of \c param is left
461 // unchanged (i.e., is the same as the one at the beginning
462 // of the function's execution), and
463 // - the state of \c param only after reading the whole input.
465 std::vector<RealT> probs;
466 std::vector<RealT> rates;
468 // Reads probability and rate vectors
469 detail::read_vector(is, probs);
475 detail::read_vector(is, rates);
481 // Update the state of the param_type object
482 if (probs.size() > 0)
484 param.probs_.swap(probs);
487 if (rates.size() > 0)
489 param.rates_.swap(rates);
495 // Adjust vector sizes (if needed)
496 if (param.probs_.size() != param.rates_.size()
497 || param.probs_.size() == 0)
501 const std::size_t np = param.probs_.size();
502 const std::size_t nr = param.rates_.size();
506 param.rates_.resize(np, 1);
510 param.probs_.resize(nr, 1);
514 param.probs_.resize(1, 1);
515 param.rates_.resize(1, 1);
519 // Normalize probabilities
520 // NOTE: this cannot be done earlier since the probability vector
521 // can be changed due to size conformance
522 hyperexp_detail::normalize(param.probs_);
524 // Set the error state in the underlying stream in case of invalid input
527 // This throws an exception if ios_base::exception(failbit) is enabled
528 is.setstate(std::ios_base::failbit);
531 //post: vector size conformance
532 assert(param.probs_.size() == param.rates_.size());
537 /** Returns true if the two sets of parameters are the same. */
538 public: BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
540 return lhs.probs_ == rhs.probs_
541 && lhs.rates_ == rhs.rates_;
544 /** Returns true if the two sets of parameters are the different. */
545 public: BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)
548 private: std::vector<RealT> probs_; ///< The <em>phase probability vector</em> parameter of the distribution
549 private: std::vector<RealT> rates_; ///< The <em>rate vector</em> parameter of the distribution
554 * Constructs a 1-phase \c hyperexponential_distribution (i.e., an
555 * exponential distribution) with rate 1.
557 public: hyperexponential_distribution()
558 : dd_(std::vector<RealT>(1, 1)),
565 * Constructs a \c hyperexponential_distribution from the <em>phase
566 * probability vector</em> and <em>rate vector</em> parameters of the
569 * The <em>phase probability vector</em> parameter is given by the range
570 * defined by [\a prob_first, \a prob_last) iterator pair, and the
571 * <em>rate vector</em> parameter is given by the range defined by
572 * [\a rate_first, \a rate_last) iterator pair.
574 * \tparam ProbIterT Must meet the requirements of \c InputIterator concept (ISO,2014,sec. 24.2.3 [input.iterators]).
575 * \tparam RateIterT Must meet the requirements of \c InputIterator concept (ISO,2014,sec. 24.2.3 [input.iterators]).
577 * \param prob_first The iterator to the beginning of the range of non-negative real elements representing the phase probabilities; if elements don't sum to 1, they are normalized.
578 * \param prob_last The iterator to the ending of the range of non-negative real elements representing the phase probabilities; if elements don't sum to 1, they are normalized.
579 * \param rate_first The iterator to the beginning of the range of non-negative real elements representing the rates.
580 * \param rate_last The iterator to the ending of the range of non-negative real elements representing the rates.
583 * -# ISO, <em>ISO/IEC 14882-2014: Information technology - Programming languages - C++</em>, 2014
586 public: template <typename ProbIterT, typename RateIterT>
587 hyperexponential_distribution(ProbIterT prob_first, ProbIterT prob_last,
588 RateIterT rate_first, RateIterT rate_last)
589 : dd_(prob_first, prob_last),
590 rates_(rate_first, rate_last)
592 assert( hyperexp_detail::check_params(dd_.probabilities(), rates_) );
596 * Constructs a \c hyperexponential_distribution from the <em>phase
597 * probability vector</em> and <em>rate vector</em> parameters of the
600 * The <em>phase probability vector</em> parameter is given by the range
601 * defined by \a prob_range, and the <em>rate vector</em> parameter is
602 * given by the range defined by \a rate_range.
604 * \tparam ProbRangeT Must meet the requirements of <a href="boost:/libs/range/doc/html/range/concepts.html">Range</a> concept.
605 * \tparam RateRangeT Must meet the requirements of <a href="boost:/libs/range/doc/html/range/concepts.html">Range</a> concept.
607 * \param prob_range The range of non-negative real elements representing the phase probabilities; if elements don't sum to 1, they are normalized.
608 * \param rate_range The range of positive real elements representing the rates.
611 * The final \c disable_if parameter is an implementation detail that
612 * differentiates between this two argument constructor and the
613 * iterator-based two argument constructor described below.
615 // We SFINAE this out of existance if either argument type is
616 // incrementable as in that case the type is probably an iterator:
617 public: template <typename ProbRangeT, typename RateRangeT>
618 hyperexponential_distribution(ProbRangeT const& prob_range,
619 RateRangeT const& rate_range,
620 typename boost::disable_if_c<boost::has_pre_increment<ProbRangeT>::value || boost::has_pre_increment<RateRangeT>::value>::type* = 0)
622 rates_(boost::begin(rate_range), boost::end(rate_range))
624 assert( hyperexp_detail::check_params(dd_.probabilities(), rates_) );
628 * Constructs a \c hyperexponential_distribution from the <em>rate
629 * vector</em> parameter of the distribution and with equal phase
632 * The <em>rate vector</em> parameter is given by the range defined by
633 * [\a rate_first, \a rate_last) iterator pair, and the <em>phase
634 * probability vector</em> parameter is set to the equal phase
635 * probabilities (i.e., to a vector of the same length \f$k\f$ of the
636 * <em>rate vector</em> and with each element set to \f$1.0/k\f$).
638 * \tparam RateIterT Must meet the requirements of \c InputIterator concept (ISO,2014,sec. 24.2.3 [input.iterators]).
639 * \tparam RateIterT2 Must meet the requirements of \c InputIterator concept (ISO,2014,sec. 24.2.3 [input.iterators]).
641 * \param rate_first The iterator to the beginning of the range of non-negative real elements representing the rates.
642 * \param rate_last The iterator to the ending of the range of non-negative real elements representing the rates.
645 * The final \c disable_if parameter is an implementation detail that
646 * differentiates between this two argument constructor and the
647 * range-based two argument constructor described above.
650 * -# ISO, <em>ISO/IEC 14882-2014: Information technology - Programming languages - C++</em>, 2014
653 // We SFINAE this out of existance if the argument type is
654 // incrementable as in that case the type is probably an iterator.
655 public: template <typename RateIterT>
656 hyperexponential_distribution(RateIterT rate_first,
658 typename boost::enable_if_c<boost::has_pre_increment<RateIterT>::value>::type* = 0)
659 : dd_(std::vector<RealT>(std::distance(rate_first, rate_last), 1)),
660 rates_(rate_first, rate_last)
662 assert( hyperexp_detail::check_params(dd_.probabilities(), rates_) );
666 * Constructs a @c param_type from the "rates" parameters
667 * of the distribution and with equal phase probabilities.
669 * The <em>rate vector</em> parameter is given by the range defined by
670 * \a rate_range, and the <em>phase probability vector</em> parameter is
671 * set to the equal phase probabilities (i.e., to a vector of the same
672 * length \f$k\f$ of the <em>rate vector</em> and with each element set
675 * \tparam RateRangeT Must meet the requirements of <a href="boost:/libs/range/doc/html/range/concepts.html">Range</a> concept.
677 * \param rate_range The range of positive real elements representing the rates.
679 public: template <typename RateRangeT>
680 hyperexponential_distribution(RateRangeT const& rate_range)
681 : dd_(std::vector<RealT>(boost::size(rate_range), 1)),
682 rates_(boost::begin(rate_range), boost::end(rate_range))
684 assert( hyperexp_detail::check_params(dd_.probabilities(), rates_) );
688 * Constructs a \c hyperexponential_distribution from its parameters.
690 * \param param The parameters of the distribution.
692 public: explicit hyperexponential_distribution(param_type const& param)
693 : dd_(param.probabilities()),
694 rates_(param.rates())
696 assert( hyperexp_detail::check_params(dd_.probabilities(), rates_) );
699 #ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
701 * Constructs a \c hyperexponential_distribution from the <em>phase
702 * probability vector</em> and <em>rate vector</em> parameters of the
705 * The <em>phase probability vector</em> parameter is given by the
706 * <em>brace-init-list</em> (ISO,2014,sec. 8.5.4 [dcl.init.list])
707 * defined by \a l1, and the <em>rate vector</em> parameter is given by the
708 * <em>brace-init-list</em> (ISO,2014,sec. 8.5.4 [dcl.init.list])
711 * \param l1 The initializer list for inizializing the phase probability vector.
712 * \param l2 The initializer list for inizializing the rate vector.
715 * -# ISO, <em>ISO/IEC 14882-2014: Information technology - Programming languages - C++</em>, 2014
718 public: hyperexponential_distribution(std::initializer_list<RealT> const& l1, std::initializer_list<RealT> const& l2)
719 : dd_(l1.begin(), l1.end()),
720 rates_(l2.begin(), l2.end())
722 assert( hyperexp_detail::check_params(dd_.probabilities(), rates_) );
726 * Constructs a \c hyperexponential_distribution from the <em>rate
727 * vector</em> parameter of the distribution and with equal phase
730 * The <em>rate vector</em> parameter is given by the
731 * <em>brace-init-list</em> (ISO,2014,sec. 8.5.4 [dcl.init.list])
732 * defined by \a l1, and the <em>phase probability vector</em> parameter is
733 * set to the equal phase probabilities (i.e., to a vector of the same
734 * length \f$k\f$ of the <em>rate vector</em> and with each element set
737 * \param l1 The initializer list for inizializing the rate vector.
740 * -# ISO, <em>ISO/IEC 14882-2014: Information technology - Programming languages - C++</em>, 2014
743 public: hyperexponential_distribution(std::initializer_list<RealT> const& l1)
744 : dd_(std::vector<RealT>(std::distance(l1.begin(), l1.end()), 1)),
745 rates_(l1.begin(), l1.end())
747 assert( hyperexp_detail::check_params(dd_.probabilities(), rates_) );
752 * Gets a random variate distributed according to the
753 * hyperexponential distribution.
755 * \tparam URNG Must meet the requirements of \uniform_random_number_generator.
757 * \param urng A uniform random number generator object.
759 * \return A random variate distributed according to the hyperexponential distribution.
761 public: template<class URNG>\
762 RealT operator()(URNG& urng) const
764 const int i = dd_(urng);
766 return boost::random::exponential_distribution<RealT>(rates_[i])(urng);
770 * Gets a random variate distributed according to the hyperexponential
771 * distribution with parameters specified by \c param.
773 * \tparam URNG Must meet the requirements of \uniform_random_number_generator.
775 * \param urng A uniform random number generator object.
776 * \param param A distribution parameter object.
778 * \return A random variate distributed according to the hyperexponential distribution.
779 * distribution with parameters specified by \c param.
781 public: template<class URNG>
782 RealT operator()(URNG& urng, const param_type& param) const
784 return hyperexponential_distribution(param)(urng);
787 /** Returns the number of phases of the distribution. */
788 public: std::size_t num_phases() const
790 return rates_.size();
793 /** Returns the <em>phase probability vector</em> parameter of the distribution. */
794 public: std::vector<RealT> probabilities() const
796 return dd_.probabilities();
799 /** Returns the <em>rate vector</em> parameter of the distribution. */
800 public: std::vector<RealT> rates() const
805 /** Returns the smallest value that the distribution can produce. */
806 public: RealT min BOOST_PREVENT_MACRO_SUBSTITUTION () const
811 /** Returns the largest value that the distribution can produce. */
812 public: RealT max BOOST_PREVENT_MACRO_SUBSTITUTION () const
814 return std::numeric_limits<RealT>::infinity();
817 /** Returns the parameters of the distribution. */
818 public: param_type param() const
820 std::vector<RealT> probs = dd_.probabilities();
822 return param_type(probs.begin(), probs.end(), rates_.begin(), rates_.end());
825 /** Sets the parameters of the distribution. */
826 public: void param(param_type const& param)
828 dd_.param(typename boost::random::discrete_distribution<int,RealT>::param_type(param.probabilities()));
829 rates_ = param.rates();
833 * Effects: Subsequent uses of the distribution do not depend
834 * on values produced by any engine prior to invoking reset.
841 /** Writes an @c hyperexponential_distribution to a @c std::ostream. */
842 public: BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, hyperexponential_distribution, hd)
848 /** Reads an @c hyperexponential_distribution from a @c std::istream. */
849 public: BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, hyperexponential_distribution, hd)
860 * Returns true if the two instances of @c hyperexponential_distribution will
861 * return identical sequences of values given equal generators.
863 public: BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(hyperexponential_distribution, lhs, rhs)
865 return lhs.dd_ == rhs.dd_
866 && lhs.rates_ == rhs.rates_;
870 * Returns true if the two instances of @c hyperexponential_distribution will
871 * return different sequences of values given equal generators.
873 public: BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(hyperexponential_distribution)
876 private: boost::random::discrete_distribution<int,RealT> dd_; ///< The \c discrete_distribution used to sample the phase probability and choose the rate
877 private: std::vector<RealT> rates_; ///< The <em>rate vector</em> parameter of the distribution
878 }; // hyperexponential_distribution
880 }} // namespace boost::random
883 #endif // BOOST_RANDOM_HYPEREXPONENTIAL_DISTRIBUTION_HPP