libs/math/test/test_hyperexponential_dist.cpp

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