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