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.
5 // (See accompanying file LICENSE_1_0.txt
6 // or copy at http://www.boost.org/LICENSE_1_0.txt)
7 //
8
9 #include <algorithm>
10 #include <boost/math/concepts/real_concept.hpp>
11 #include <boost/math/distributions/complement.hpp>
12 #include <boost/math/distributions/hyperexponential.hpp>
13 #include <boost/math/tools/precision.hpp>
14
15 #define BOOST_TEST_MAIN
16 #include <boost/test/unit_test.hpp>
17 #include <boost/test/floating_point_comparison.hpp>
18
19 #include <cstddef>
20 #include <iostream>
21 #include <vector>
22
23 #define BOOST_MATH_HYPEREXP_CHECK_CLOSE_COLLECTIONS(T, actual, expected, tol) \
24 do { \
25 std::vector<T> x = (actual); \
26 std::vector<T> y = (expected); \
27 BOOST_CHECK_EQUAL( x.size(), y.size() ); \
28 const std::size_t n = x.size(); \
29 for (std::size_t i = 0; i < n; ++i) \
30 { \
31 BOOST_CHECK_CLOSE( x[i], y[i], tol ); \
32 } \
33 } while(false)
34
35 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
36 typedef boost::mpl::list<float, double, long double, boost::math::concepts::real_concept> test_types;
37 #else
38 typedef boost::mpl::list<float, double> test_types;
39 #endif
40
41 template <typename RealT>
make_tolerance()42 RealT make_tolerance()
43 {
44 // Tolerance is 100eps expressed as a persentage (as required by Boost.Build):
45 return boost::math::tools::epsilon<RealT>() * 100 * 100;
46 }
47
BOOST_AUTO_TEST_CASE_TEMPLATE(klass,RealT,test_types)48 BOOST_AUTO_TEST_CASE_TEMPLATE(klass, RealT, test_types)
49 {
50 const RealT tol = make_tolerance<RealT>();
51
52 boost::math::hyperexponential_distribution<RealT> dist;
53 BOOST_CHECK_EQUAL(dist.num_phases(), 1);
54 BOOST_CHECK_CLOSE(dist.probabilities()[0], static_cast<RealT>(1L), tol);
55 BOOST_CHECK_CLOSE(dist.rates()[0], static_cast<RealT>(1L), tol);
56
57 const RealT probs[] = { static_cast<RealT>(0.2L), static_cast<RealT>(0.3L), static_cast<RealT>(0.5L) };
58 const RealT rates[] = { static_cast<RealT>(0.5L), static_cast<RealT>(1.0L), static_cast<RealT>(1.5L) };
59 const std::size_t n = sizeof(probs) / sizeof(RealT);
60
61 boost::math::hyperexponential_distribution<RealT> dist_it(probs, probs+n, rates, rates+n);
62 BOOST_CHECK_EQUAL(dist_it.num_phases(), n);
63 BOOST_MATH_HYPEREXP_CHECK_CLOSE_COLLECTIONS(RealT, dist_it.probabilities(), std::vector<RealT>(probs, probs+n), tol);
64 BOOST_MATH_HYPEREXP_CHECK_CLOSE_COLLECTIONS(RealT, dist_it.rates(), std::vector<RealT>(rates, rates+n), tol);
65
66 boost::math::hyperexponential_distribution<RealT> dist_r(probs, rates);
67 BOOST_CHECK_EQUAL(dist_r.num_phases(), n);
68 BOOST_MATH_HYPEREXP_CHECK_CLOSE_COLLECTIONS(RealT, dist_r.probabilities(), std::vector<RealT>(probs, probs+n), tol);
69 BOOST_MATH_HYPEREXP_CHECK_CLOSE_COLLECTIONS(RealT, dist_r.rates(), std::vector<RealT>(rates, rates+n), tol);
70
71 #if !defined(BOOST_NO_CXX11_HDR_INITIALIZER_LIST) && !(defined(BOOST_GCC_VERSION) && (BOOST_GCC_VERSION < 40500))
72 boost::math::hyperexponential_distribution<RealT> dist_il = {{static_cast<RealT>(0.2L), static_cast<RealT>(0.3L), static_cast<RealT>(0.5L)}, {static_cast<RealT>(0.5L), static_cast<RealT>(1.0L), static_cast<RealT>(1.5L)}};
73 BOOST_CHECK_EQUAL(dist_il.num_phases(), n);
74 BOOST_MATH_HYPEREXP_CHECK_CLOSE_COLLECTIONS(RealT, dist_il.probabilities(), std::vector<RealT>(probs, probs+n), tol);
75 BOOST_MATH_HYPEREXP_CHECK_CLOSE_COLLECTIONS(RealT, dist_il.rates(), std::vector<RealT>(rates, rates+n), tol);
76
77 boost::math::hyperexponential_distribution<RealT> dist_n_r = { static_cast<RealT>(0.5L), static_cast<RealT>(1.0L), static_cast<RealT>(1.5L) };
78 BOOST_CHECK_EQUAL(dist_n_r.num_phases(), n);
79 BOOST_MATH_HYPEREXP_CHECK_CLOSE_COLLECTIONS(RealT, dist_n_r.probabilities(), std::vector<RealT>(n, static_cast<RealT>(1.0L / 3.0L)), tol);
80 BOOST_MATH_HYPEREXP_CHECK_CLOSE_COLLECTIONS(RealT, dist_n_r.rates(), std::vector<RealT>(rates, rates + n), tol);
81 #endif // BOOST_NO_CXX11_HDR_INITIALIZER_LIST
82
83 boost::math::hyperexponential_distribution<RealT> dist_n_it(rates, rates+n);
84 BOOST_CHECK_EQUAL(dist_n_it.num_phases(), n);
85 BOOST_MATH_HYPEREXP_CHECK_CLOSE_COLLECTIONS(RealT, dist_n_it.probabilities(), std::vector<RealT>(n, static_cast<RealT>(1.0L/3.0L)), tol);
86 BOOST_MATH_HYPEREXP_CHECK_CLOSE_COLLECTIONS(RealT, dist_n_it.rates(), std::vector<RealT>(rates, rates+n), tol);
87
88 boost::math::hyperexponential_distribution<RealT> dist_n_r2(rates);
89 BOOST_CHECK_EQUAL(dist_n_r2.num_phases(), n);
90 BOOST_MATH_HYPEREXP_CHECK_CLOSE_COLLECTIONS(RealT, dist_n_r2.probabilities(), std::vector<RealT>(n, static_cast<RealT>(1.0L/3.0L)), tol);
91 BOOST_MATH_HYPEREXP_CHECK_CLOSE_COLLECTIONS(RealT, dist_n_r2.rates(), std::vector<RealT>(rates, rates+n), tol);
92 }
93
BOOST_AUTO_TEST_CASE_TEMPLATE(range,RealT,test_types)94 BOOST_AUTO_TEST_CASE_TEMPLATE(range, RealT, test_types)
95 {
96 const RealT tol = make_tolerance<RealT>();
97
98 const RealT probs[] = { static_cast<RealT>(0.2L), static_cast<RealT>(0.3L), static_cast<RealT>(0.5L) };
99 const RealT rates[] = { static_cast<RealT>(0.5L), static_cast<RealT>(1.0L), static_cast<RealT>(1.5L) };
100 const std::size_t n = sizeof(probs) / sizeof(RealT);
101
102 boost::math::hyperexponential_distribution<RealT> dist(probs, probs+n, rates, rates+n);
103
104 std::pair<RealT,RealT> res;
105 res = boost::math::range(dist);
106
107 BOOST_CHECK_CLOSE( res.first, static_cast<RealT>(0), tol );
108 if(std::numeric_limits<RealT>::has_infinity)
109 {
110 BOOST_CHECK_EQUAL(res.second, std::numeric_limits<RealT>::infinity());
111 }
112 else
113 {
114 BOOST_CHECK_EQUAL(res.second, boost::math::tools::max_value<RealT>());
115 }
116 }
117
BOOST_AUTO_TEST_CASE_TEMPLATE(support,RealT,test_types)118 BOOST_AUTO_TEST_CASE_TEMPLATE(support, RealT, test_types)
119 {
120 const RealT tol = make_tolerance<RealT>();
121
122 const RealT probs[] = { static_cast<RealT>(0.2L), static_cast<RealT>(0.3L), static_cast<RealT>(0.5L) };
123 const RealT rates[] = { static_cast<RealT>(0.5L), static_cast<RealT>(1), static_cast<RealT>(1.5L) };
124 const std::size_t n = sizeof(probs)/sizeof(RealT);
125
126 boost::math::hyperexponential_distribution<RealT> dist(probs, probs+n, rates, rates+n);
127
128 std::pair<RealT,RealT> res;
129 res = boost::math::support(dist);
130
131 BOOST_CHECK_CLOSE( res.first, boost::math::tools::min_value<RealT>(), tol );
132 BOOST_CHECK_CLOSE( res.second, boost::math::tools::max_value<RealT>(), tol );
133 }
134
BOOST_AUTO_TEST_CASE_TEMPLATE(pdf,RealT,test_types)135 BOOST_AUTO_TEST_CASE_TEMPLATE(pdf, RealT, test_types)
136 {
137 const RealT tol = make_tolerance<RealT>();
138
139 const RealT probs[] = { static_cast<RealT>(0.2L), static_cast<RealT>(0.3L), static_cast<RealT>(0.5L) };
140 const RealT rates[] = { static_cast<RealT>(0.5L), static_cast<RealT>(1), static_cast<RealT>(1.5) };
141 const std::size_t n = sizeof(probs)/sizeof(RealT);
142
143 boost::math::hyperexponential_distribution<RealT> dist(probs, probs+n, rates, rates+n);
144
145 // Mathematica: Table[N[PDF[HyperexponentialDistribution[{1/5, 3/10, 1/2}, {1/2, 1, 3/2}], x], 35], {x, 0, 4}]
146 BOOST_CHECK_CLOSE( boost::math::pdf(dist, static_cast<RealT>(0)), static_cast<RealT>(1.15L), tol );
147 BOOST_CHECK_CLOSE( boost::math::pdf(dist, static_cast<RealT>(1)), static_cast<RealT>(0.33836451843401841053899743762056570L), tol );
148 BOOST_CHECK_CLOSE( boost::math::pdf(dist, static_cast<RealT>(2)), static_cast<RealT>(0.11472883036402599696225903724543774L), tol );
149 BOOST_CHECK_CLOSE( boost::math::pdf(dist, static_cast<RealT>(3)), static_cast<RealT>(0.045580883928883895659238122486617681L), tol );
150 BOOST_CHECK_CLOSE( boost::math::pdf(dist, static_cast<RealT>(4)), static_cast<RealT>(0.020887284122781292094799231452333314L), tol );
151 }
152
BOOST_AUTO_TEST_CASE_TEMPLATE(cdf,RealT,test_types)153 BOOST_AUTO_TEST_CASE_TEMPLATE(cdf, RealT, test_types)
154 {
155 const RealT tol = make_tolerance<RealT>();
156
157 const RealT probs[] = { static_cast<RealT>(0.2L), static_cast<RealT>(0.3L), static_cast<RealT>(0.5L) };
158 const RealT rates[] = { static_cast<RealT>(0.5L), static_cast<RealT>(1.0L), static_cast<RealT>(1.5L) };
159 const std::size_t n = sizeof(probs)/sizeof(RealT);
160
161 boost::math::hyperexponential_distribution<RealT> dist(probs, probs+n, rates, rates+n);
162
163 // Mathematica: Table[N[CDF[HyperexponentialDistribution[{1/5, 3/10, 1/2}, {1/2, 1, 3/2}], x], 35], {x, 0, 4}]
164 BOOST_CHECK_CLOSE( boost::math::cdf(dist, static_cast<RealT>(0)), static_cast<RealT>(0), tol );
165 BOOST_CHECK_CLOSE( boost::math::cdf(dist, static_cast<RealT>(1)), static_cast<RealT>(0.65676495563182570433394272657131939L), tol );
166 BOOST_CHECK_CLOSE( boost::math::cdf(dist, static_cast<RealT>(2)), static_cast<RealT>(0.86092999261079575662302418965093162L), tol );
167 BOOST_CHECK_CLOSE( boost::math::cdf(dist, static_cast<RealT>(3)), static_cast<RealT>(0.93488334919083369807146961400871370L), tol );
168 BOOST_CHECK_CLOSE( boost::math::cdf(dist, static_cast<RealT>(4)), static_cast<RealT>(0.96619887559772402832156211090812241L), tol );
169 }
170
171
BOOST_AUTO_TEST_CASE_TEMPLATE(quantile,RealT,test_types)172 BOOST_AUTO_TEST_CASE_TEMPLATE(quantile, RealT, test_types)
173 {
174 const RealT tol = make_tolerance<RealT>();
175
176 const RealT probs[] = { static_cast<RealT>(0.2L), static_cast<RealT>(0.3L), static_cast<RealT>(0.5L) };
177 const RealT rates[] = { static_cast<RealT>(0.5L), static_cast<RealT>(1.0L), static_cast<RealT>(1.5L) };
178 const std::size_t n = sizeof(probs)/sizeof(RealT);
179
180 boost::math::hyperexponential_distribution<RealT> dist(probs, probs+n, rates, rates+n);
181
182 // Mathematica: Table[N[Quantile[HyperexponentialDistribution[{1/5, 3/10, 1/2}, {1/2, 1, 3/2}], p], 35], {p, {0.`35, 0.6567649556318257043339427265713193884067872189124925936717`35, 0.8609299926107957566230241896509316171726985139265620607067`35, 0.9348833491908336980714696140087136988562861627183715044229`35, 0.9661988755977240283215621109081224127091468307592751727719`35}}]
183 BOOST_CHECK_CLOSE( boost::math::quantile(dist, static_cast<RealT>(0)), static_cast<RealT>(0), tol );
184 BOOST_CHECK_CLOSE( boost::math::quantile(dist, static_cast<RealT>(0.65676495563182570433394272657131939L)), static_cast<RealT>(1), tol );
185 BOOST_CHECK_CLOSE( boost::math::quantile(dist, static_cast<RealT>(0.86092999261079575662302418965093162L)), static_cast<RealT>(2), tol );
186 BOOST_CHECK_CLOSE( boost::math::quantile(dist, static_cast<RealT>(0.93488334919083369807146961400871370L)), static_cast<RealT>(3), tol );
187 BOOST_CHECK_CLOSE( boost::math::quantile(dist, static_cast<RealT>(0.96619887559772402832156211090812241L)), static_cast<RealT>(4), tol );
188 }
189
BOOST_AUTO_TEST_CASE_TEMPLATE(ccdf,RealT,test_types)190 BOOST_AUTO_TEST_CASE_TEMPLATE(ccdf, RealT, test_types)
191 {
192 const RealT tol = make_tolerance<RealT>();
193
194 const RealT probs[] = { static_cast<RealT>(0.2L), static_cast<RealT>(0.3L), static_cast<RealT>(0.5L) };
195 const RealT rates[] = { static_cast<RealT>(0.5L), static_cast<RealT>(1.0L), static_cast<RealT>(1.5L) };
196 const std::size_t n = sizeof(probs)/sizeof(RealT);
197
198 boost::math::hyperexponential_distribution<RealT> dist(probs, probs+n, rates, rates+n);
199
200 // Mathematica: Table[N[SurvivalFunction[HyperexponentialDistribution[{1/5, 3/10, 1/2}, {1/2, 1, 3/2}], x], 35], {x, 0, 4}]
201 BOOST_CHECK_CLOSE( boost::math::cdf(boost::math::complement(dist, static_cast<RealT>(0))), static_cast<RealT>(1), tol );
202 BOOST_CHECK_CLOSE( boost::math::cdf(boost::math::complement(dist, static_cast<RealT>(1))), static_cast<RealT>(0.34323504436817429566605727342868061L), tol );
203 BOOST_CHECK_CLOSE( boost::math::cdf(boost::math::complement(dist, static_cast<RealT>(2))), static_cast<RealT>(0.13907000738920424337697581034906838L), tol );
204 BOOST_CHECK_CLOSE( boost::math::cdf(boost::math::complement(dist, static_cast<RealT>(3))), static_cast<RealT>(0.065116650809166301928530385991286301L), tol );
205 BOOST_CHECK_CLOSE( boost::math::cdf(boost::math::complement(dist, static_cast<RealT>(4))), static_cast<RealT>(0.033801124402275971678437889091877587L), tol );
206 }
207
208
BOOST_AUTO_TEST_CASE_TEMPLATE(cquantile,RealT,test_types)209 BOOST_AUTO_TEST_CASE_TEMPLATE(cquantile, RealT, test_types)
210 {
211 const RealT tol = make_tolerance<RealT>();
212
213 const RealT probs[] = { static_cast<RealT>(0.2L), static_cast<RealT>(0.3L), static_cast<RealT>(0.5L) };
214 const RealT rates[] = { static_cast<RealT>(0.5L), static_cast<RealT>(1.0L), static_cast<RealT>(1.5L) };
215 const std::size_t n = sizeof(probs) / sizeof(RealT);
216
217 boost::math::hyperexponential_distribution<RealT> dist(probs, probs+n, rates, rates+n);
218
219 // Mathematica: Table[N[InverseSurvivalFunction[HyperexponentialDistribution[{1/5, 3/10, 1/2}, {1/2, 1, 3/2}], p], 35], {p, {1.`35, 0.3432350443681742956660572734286806115932127810875074063283`35, 0.1390700073892042433769758103490683828273014860734379392933`35, 0.0651166508091663019285303859912863011437138372816284955771`35, 0.0338011244022759716784378890918775872908531692407248272281`35}}]
220 BOOST_CHECK_CLOSE( boost::math::quantile(boost::math::complement(dist, static_cast<RealT>(1))), static_cast<RealT>(0), tol );
221 BOOST_CHECK_CLOSE( boost::math::quantile(boost::math::complement(dist, static_cast<RealT>(0.34323504436817429566605727342868061L))), static_cast<RealT>(1), tol );
222 BOOST_CHECK_CLOSE( boost::math::quantile(boost::math::complement(dist, static_cast<RealT>(0.13907000738920424337697581034906838L))), static_cast<RealT>(2), tol );
223 BOOST_CHECK_CLOSE( boost::math::quantile(boost::math::complement(dist, static_cast<RealT>(0.065116650809166301928530385991286301L))), static_cast<RealT>(3), tol );
224 BOOST_CHECK_CLOSE( boost::math::quantile(boost::math::complement(dist, static_cast<RealT>(0.033801124402275971678437889091877587L))), static_cast<RealT>(4), tol );
225 }
226
BOOST_AUTO_TEST_CASE_TEMPLATE(mean,RealT,test_types)227 BOOST_AUTO_TEST_CASE_TEMPLATE(mean, RealT, test_types)
228 {
229 const RealT tol = make_tolerance<RealT>();
230
231 const RealT probs[] = { static_cast<RealT>(0.2L), static_cast<RealT>(0.3L), static_cast<RealT>(0.5L) };
232 const RealT rates[] = { static_cast<RealT>(0.5L), static_cast<RealT>(1.0L), static_cast<RealT>(1.5L) };
233 const std::size_t n = sizeof(probs) / sizeof(RealT);
234
235 boost::math::hyperexponential_distribution<RealT> dist(probs, probs+n, rates, rates+n);
236
237 // Mathematica: N[Mean[HyperexponentialDistribution[{1/5, 3/10, 1/2}, {1/2, 1, 3/2}]], 35]
238 BOOST_CHECK_CLOSE( boost::math::mean(dist), static_cast<RealT>(1.0333333333333333333333333333333333L), tol );
239 }
240
BOOST_AUTO_TEST_CASE_TEMPLATE(variance,RealT,test_types)241 BOOST_AUTO_TEST_CASE_TEMPLATE(variance, RealT, test_types)
242 {
243 const RealT tol = make_tolerance<RealT>();
244
245 const RealT probs[] = { static_cast<RealT>(0.2L), static_cast<RealT>(0.3L), static_cast<RealT>(0.5L) };
246 const RealT rates[] = { static_cast<RealT>(0.5L), static_cast<RealT>(1.0L), static_cast<RealT>(1.5L) };
247 const std::size_t n = sizeof(probs) / sizeof(RealT);
248
249 boost::math::hyperexponential_distribution<RealT> dist(probs, probs+n, rates, rates+n);
250
251 // Mathematica: N[Variance[HyperexponentialDistribution[{1/5, 3/10, 1/2}, {1/2, 1, 3/2}]], 35]
252 BOOST_CHECK_CLOSE( boost::math::variance(dist), static_cast<RealT>(1.5766666666666666666666666666666667L), tol );
253 }
254
BOOST_AUTO_TEST_CASE_TEMPLATE(kurtosis,RealT,test_types)255 BOOST_AUTO_TEST_CASE_TEMPLATE(kurtosis, RealT, test_types)
256 {
257 const RealT tol = make_tolerance<RealT>();
258
259 const RealT probs[] = { static_cast<RealT>(0.2L), static_cast<RealT>(0.3L), static_cast<RealT>(0.5L) };
260 const RealT rates[] = { static_cast<RealT>(0.5L), static_cast<RealT>(1.0L), static_cast<RealT>(1.5L) };
261 const std::size_t n = sizeof(probs) / sizeof(RealT);
262
263 boost::math::hyperexponential_distribution<RealT> dist(probs, probs+n, rates, rates+n);
264
265 // Mathematica: N[Kurtosis[HyperexponentialDistribution[{1/5, 3/10, 1/2}, {1/2, 1, 3/2}]], 35]
266 BOOST_CHECK_CLOSE( boost::math::kurtosis(dist), static_cast<RealT>(19.750738616808728416968743435138046L), tol );
267 // Mathematica: N[Kurtosis[HyperexponentialDistribution[{1/5, 3/10, 1/2}, {1/2, 1, 3/2}] - 3.`35], 35]
268 BOOST_CHECK_CLOSE( boost::math::kurtosis_excess(dist), static_cast<RealT>(16.750738616808728416968743435138046L), tol );
269 }
270
BOOST_AUTO_TEST_CASE_TEMPLATE(skewness,RealT,test_types)271 BOOST_AUTO_TEST_CASE_TEMPLATE(skewness, RealT, test_types)
272 {
273 const RealT tol = make_tolerance<RealT>();
274
275 const RealT probs[] = { static_cast<RealT>(0.2L), static_cast<RealT>(0.3L), static_cast<RealT>(0.5L) };
276 const RealT rates[] = { static_cast<RealT>(0.5L), static_cast<RealT>(1.0L), static_cast<RealT>(1.5L) };
277 const std::size_t n = sizeof(probs) / sizeof(RealT);
278
279 boost::math::hyperexponential_distribution<RealT> dist(probs, probs+n, rates, rates+n);
280
281 // Mathematica: N[Skewness[HyperexponentialDistribution[{1/5, 3/10, 1/2}, {1/2, 1, 3/2}]], 35]
282 BOOST_CHECK_CLOSE( boost::math::skewness(dist), static_cast<RealT>(3.1811387449963809211146099116375685L), tol );
283 }
284
BOOST_AUTO_TEST_CASE_TEMPLATE(mode,RealT,test_types)285 BOOST_AUTO_TEST_CASE_TEMPLATE(mode, RealT, test_types)
286 {
287 const RealT tol = make_tolerance<RealT>();
288
289 const RealT probs[] = { static_cast<RealT>(0.2L), static_cast<RealT>(0.3L), static_cast<RealT>(0.5L) };
290 const RealT rates[] = { static_cast<RealT>(0.5L), static_cast<RealT>(1.0L), static_cast<RealT>(1.5L) };
291 const std::size_t n = sizeof(probs) / sizeof(RealT);
292
293 boost::math::hyperexponential_distribution<RealT> dist(probs, probs+n, rates, rates+n);
294
295 BOOST_CHECK_CLOSE( boost::math::mode(dist), static_cast<RealT>(0), tol );
296 }
297
298 template <class T>
f(T t)299 void f(T t)
300 {
301 std::cout << typeid(t).name() << std::endl;
302 }
303
BOOST_AUTO_TEST_CASE(construct)304 BOOST_AUTO_TEST_CASE(construct)
305 {
306 boost::array<double, 3> da1 = { { 0.5, 1, 1.5 } };
307 boost::array<double, 3> da2 = { { 0.25, 0.5, 0.25 } };
308 std::vector<double> v1(da1.begin(), da1.end());
309 std::vector<double> v2(da2.begin(), da2.end());
310
311 std::vector<double> result_v;
312 boost::math::hyperexponential he1(v2, v1);
313 result_v = he1.rates();
314 BOOST_CHECK_EQUAL_COLLECTIONS(v1.begin(), v1.end(), result_v.begin(), result_v.end());
315 result_v = he1.probabilities();
316 BOOST_CHECK_EQUAL_COLLECTIONS(v2.begin(), v2.end(), result_v.begin(), result_v.end());
317
318 boost::math::hyperexponential he2(da2, da1);
319 result_v = he2.rates();
320 BOOST_CHECK_EQUAL_COLLECTIONS(v1.begin(), v1.end(), result_v.begin(), result_v.end());
321 result_v = he2.probabilities();
322 BOOST_CHECK_EQUAL_COLLECTIONS(v2.begin(), v2.end(), result_v.begin(), result_v.end());
323
324 #if !defined(BOOST_NO_CXX11_HDR_INITIALIZER_LIST) && !(defined(BOOST_GCC_VERSION) && (BOOST_GCC_VERSION < 40500))
325 std::initializer_list<double> il = { 0.25, 0.5, 0.25 };
326 std::initializer_list<double> il2 = { 0.5, 1, 1.5 };
327 boost::math::hyperexponential he3(il, il2);
328 result_v = he3.rates();
329 BOOST_CHECK_EQUAL_COLLECTIONS(v1.begin(), v1.end(), result_v.begin(), result_v.end());
330 result_v = he3.probabilities();
331 BOOST_CHECK_EQUAL_COLLECTIONS(v2.begin(), v2.end(), result_v.begin(), result_v.end());
332
333 boost::math::hyperexponential he4({ 0.25, 0.5, 0.25 }, { 0.5, 1.0, 1.5 });
334 result_v = he4.rates();
335 BOOST_CHECK_EQUAL_COLLECTIONS(v1.begin(), v1.end(), result_v.begin(), result_v.end());
336 result_v = he4.probabilities();
337 BOOST_CHECK_EQUAL_COLLECTIONS(v2.begin(), v2.end(), result_v.begin(), result_v.end());
338 #endif
339 }
340
BOOST_AUTO_TEST_CASE_TEMPLATE(special_cases,RealT,test_types)341 BOOST_AUTO_TEST_CASE_TEMPLATE(special_cases, RealT, test_types)
342 {
343 const RealT tol = make_tolerance<RealT>();
344
345 // When the number of phases is 1, the hyperexponential distribution is an exponential distribution
346 const RealT rates1[] = { static_cast<RealT>(0.5L) };
347 boost::math::hyperexponential_distribution<RealT> hexp1(rates1);
348 boost::math::exponential_distribution<RealT> exp1(rates1[0]);
349 BOOST_CHECK_CLOSE(boost::math::pdf(hexp1, static_cast<RealT>(1L)), boost::math::pdf(exp1, static_cast<RealT>(1L)), tol);
350 BOOST_CHECK_CLOSE(boost::math::cdf(hexp1, static_cast<RealT>(1L)), boost::math::cdf(exp1, static_cast<RealT>(1L)), tol);
351 BOOST_CHECK_CLOSE(boost::math::mean(hexp1), boost::math::mean(exp1), tol);
352 BOOST_CHECK_CLOSE(boost::math::variance(hexp1), boost::math::variance(exp1), tol);
353 BOOST_CHECK_CLOSE(boost::math::quantile(hexp1, static_cast<RealT>(0.25L)), boost::math::quantile(exp1, static_cast<RealT>(0.25L)), tol);
354 BOOST_CHECK_CLOSE(boost::math::median(hexp1), boost::math::median(exp1), tol);
355 BOOST_CHECK_CLOSE(boost::math::quantile(hexp1, static_cast<RealT>(0.75L)), boost::math::quantile(exp1, static_cast<RealT>(0.75L)), tol);
356 BOOST_CHECK_CLOSE(boost::math::mode(hexp1), boost::math::mode(exp1), tol);
357
358 // When a k-phase hyperexponential distribution has all rates equal to r, the distribution is an exponential distribution with rate r
359 const RealT rate2 = static_cast<RealT>(0.5L);
360 const RealT rates2[] = { rate2, rate2, rate2 };
361 boost::math::hyperexponential_distribution<RealT> hexp2(rates2);
362 boost::math::exponential_distribution<RealT> exp2(rate2);
363 BOOST_CHECK_CLOSE(boost::math::pdf(hexp2, static_cast<RealT>(1L)), boost::math::pdf(exp2, static_cast<RealT>(1L)), tol);
364 BOOST_CHECK_CLOSE(boost::math::cdf(hexp2, static_cast<RealT>(1L)), boost::math::cdf(exp2, static_cast<RealT>(1L)), tol);
365 BOOST_CHECK_CLOSE(boost::math::mean(hexp2), boost::math::mean(exp2), tol);
366 BOOST_CHECK_CLOSE(boost::math::variance(hexp2), boost::math::variance(exp2), tol);
367 BOOST_CHECK_CLOSE(boost::math::quantile(hexp2, static_cast<RealT>(0.25L)), boost::math::quantile(exp2, static_cast<RealT>(0.25L)), tol);
368 BOOST_CHECK_CLOSE(boost::math::median(hexp2), boost::math::median(exp2), tol);
369 BOOST_CHECK_CLOSE(boost::math::quantile(hexp2, static_cast<RealT>(0.75L)), boost::math::quantile(exp2, static_cast<RealT>(0.75L)), tol);
370 BOOST_CHECK_CLOSE(boost::math::mode(hexp2), boost::math::mode(exp2), tol);
371 }
372
BOOST_AUTO_TEST_CASE_TEMPLATE(error_cases,RealT,test_types)373 BOOST_AUTO_TEST_CASE_TEMPLATE(error_cases, RealT, test_types)
374 {
375 typedef boost::math::hyperexponential_distribution<RealT> dist_t;
376 boost::array<RealT, 2> probs = { { 1, 2 } };
377 boost::array<RealT, 3> probs2 = { { 1, 2, 3 } };
378 boost::array<RealT, 3> rates = { { 1, 2, 3 } };
379 BOOST_CHECK_THROW(dist_t(probs.begin(), probs.end(), rates.begin(), rates.end()), std::domain_error);
380 BOOST_CHECK_THROW(dist_t(probs, rates), std::domain_error);
381 rates[1] = 0;
382 BOOST_CHECK_THROW(dist_t(probs2, rates), std::domain_error);
383 rates[1] = -1;
384 BOOST_CHECK_THROW(dist_t(probs2, rates), std::domain_error);
385 BOOST_CHECK_THROW(dist_t(probs.begin(), probs.begin(), rates.begin(), rates.begin()), std::domain_error);
386 BOOST_CHECK_THROW(dist_t(rates.begin(), rates.begin()), std::domain_error);
387 }