boost/math
2
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2026-01-19 04:22:09 +00:00
Code Issues Projects Releases Wiki Activity
Files
4fbfc1f2cedb1da21c0818e056fe51609728d670
math/test/test_hyperexponential_dist.cpp

348 lines
12 KiB
C++
Raw Blame History

// Copyright 2014 Marco Guazzone (marco.guazzone@gmail.com).
//
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt
// or copy at http://www.boost.org/LICENSE_1_0.txt)
//
#include <algorithm>
#include <boost/math/concepts/real_concept.hpp>
#include <boost/math/distributions/complement.hpp>
#include <boost/math/distributions/hyperexponential.hpp>
#include <boost/math/tools/precision.hpp>
#define BOOST_TEST_MAIN
#include <boost/test/unit_test.hpp>
#include <boost/test/floating_point_comparison.hpp>
#include <cstddef>
#include <iostream>
#include <vector>
template <typename RealT>
RealT make_tolerance()
{
// We must use a low precision since it seems the involved computations are very challenging from the numerical point of view.
// Indeed, both Octave 3.6.4, MATLAB 2012a and Mathematica 10 provides different results.
// E.g.:
// x = [0 1 2 3 4]
// p = [0.2 0.3 0.5]
// r = [0.5 1.0 1.5]
// PDF(x)
// - MATLAB: 1.033333333333333, 0.335636985323608, 0.135792553231720, 0.061039382459897, 0.028790027125382
// - Octave: 1.0333333333333332, 0.3356369853236084, 0.1357925532317197, 0.0610393824598966, 0.0287900271253818
// - Mathematica: 1.15, 0.3383645184340184, 0.11472883036402601, 0.04558088392888389, 0.02088728412278129
//
// (Tested under Fedora Linux 20 x86_64 running on Intel(R) Core(TM) i7-3540M)
//
/*
RealT tol = std::max(boost::math::tools::epsilon<RealT>(),
static_cast<RealT>(boost::math::tools::epsilon<double>()*5)*150);
// At float precision we need to up the tolerance, since
// the input values are rounded off to inexact quantities
// the results get thrown off by a noticeable amount.
if (boost::math::tools::digits<RealT>() < 50)
{
tol *= 50;
}
if (boost::is_floating_point<RealT>::value != 1)
{
tol *= 20; // real_concept special functions are less accurate
}
*/
const RealT tol = 1e-4;
//std::cout << "[" << __func__ << "] Tolerance: " << tol << "%" << std::endl;
return tol;
}
template <typename RealT>
void test_range()
{
const RealT tol = make_tolerance<RealT>();
const RealT probs[] = {0.2, 0.3, 0.5};
const RealT rates[] = {0.5, 1.0, 1.5};
const std::size_t n = sizeof(probs)/sizeof(RealT);
boost::math::hyperexponential_distribution<RealT> dist(probs, probs+n, rates, rates+n);
std::pair<RealT,RealT> res;
res = boost::math::range(dist);
BOOST_CHECK_CLOSE( res.first, 0, tol );
BOOST_CHECK_EQUAL( res.second, std::numeric_limits<RealT>::infinity() );
}
template <typename RealT>
void test_support()
{
const RealT tol = make_tolerance<RealT>();
const RealT probs[] = {0.2, 0.3, 0.5};
const RealT rates[] = {0.5, 1.0, 1.5};
const std::size_t n = sizeof(probs)/sizeof(RealT);
boost::math::hyperexponential_distribution<RealT> dist(probs, probs+n, rates, rates+n);
std::pair<RealT,RealT> res;
res = boost::math::support(dist);
BOOST_CHECK_CLOSE( res.first, boost::math::tools::min_value<RealT>(), tol );
BOOST_CHECK_CLOSE( res.second, boost::math::tools::max_value<RealT>(), tol );
}
template <typename RealT>
void test_pdf()
{
const RealT tol = make_tolerance<RealT>();
const RealT probs[] = {0.2, 0.3, 0.5};
const RealT rates[] = {0.5, 1.0, 1.5};
const std::size_t n = sizeof(probs)/sizeof(RealT);
boost::math::hyperexponential_distribution<RealT> dist(probs, probs+n, rates, rates+n);
// Mathematica: Table[PDF[HyperexponentialDistribution[{0.2, 0.3, 0.5}, {.5, 1.0, 1.5}], x], {x, 0, 4}]
BOOST_CHECK_CLOSE( boost::math::pdf(dist, static_cast<RealT>(0)), static_cast<RealT>(1.15), tol );
BOOST_CHECK_CLOSE( boost::math::pdf(dist, static_cast<RealT>(1)), static_cast<RealT>(0.3383645184340184), tol );
BOOST_CHECK_CLOSE( boost::math::pdf(dist, static_cast<RealT>(2)), static_cast<RealT>(0.11472883036402601), tol );
BOOST_CHECK_CLOSE( boost::math::pdf(dist, static_cast<RealT>(3)), static_cast<RealT>(0.04558088392888389), tol );
BOOST_CHECK_CLOSE( boost::math::pdf(dist, static_cast<RealT>(4)), static_cast<RealT>(0.02088728412278129), tol );
}
template <typename RealT>
void test_cdf()
{
const RealT tol = make_tolerance<RealT>();
const RealT probs[] = {0.2, 0.3, 0.5};
const RealT rates[] = {0.5, 1.0, 1.5};
const std::size_t n = sizeof(probs)/sizeof(RealT);
boost::math::hyperexponential_distribution<RealT> dist(probs, probs+n, rates, rates+n);
// Mathematica: Table[CDF[HyperexponentialDistribution[{0.2, 0.3, 0.5}, {.5, 1.0, 1.5}], x], {x, 0, 4}]
BOOST_CHECK_CLOSE( boost::math::cdf(dist, static_cast<RealT>(0)), static_cast<RealT>(0), tol );
BOOST_CHECK_CLOSE( boost::math::cdf(dist, static_cast<RealT>(1)), static_cast<RealT>(0.6567649556318257), tol );
BOOST_CHECK_CLOSE( boost::math::cdf(dist, static_cast<RealT>(2)), static_cast<RealT>(0.8609299926107957), tol );
BOOST_CHECK_CLOSE( boost::math::cdf(dist, static_cast<RealT>(3)), static_cast<RealT>(0.9348833491908337), tol );
BOOST_CHECK_CLOSE( boost::math::cdf(dist, static_cast<RealT>(4)), static_cast<RealT>(0.966198875597724), tol );
}
template <typename RealT>
void test_quantile()
{
const RealT tol = make_tolerance<RealT>();
const RealT probs[] = {0.2, 0.3, 0.5};
const RealT rates[] = {0.5, 1.0, 1.5};
const std::size_t n = sizeof(probs)/sizeof(RealT);
boost::math::hyperexponential_distribution<RealT> dist(probs, probs+n, rates, rates+n);
// Mathematica: Table[Quantile[HyperexponentialDistribution[{0.2, 0.3, 0.5}, {.5, 1.0, 1.5}], p], {p, {0, 0.6567649556318257, 0.8609299926107957, 0.9348833491908337, 0.966198875597724}}]
BOOST_CHECK_CLOSE( boost::math::quantile(dist, static_cast<RealT>(0)), static_cast<RealT>(0), tol );
BOOST_CHECK_CLOSE( boost::math::quantile(dist, static_cast<RealT>(0.6567649556318257)), static_cast<RealT>(1.0000000000000036), tol );
BOOST_CHECK_CLOSE( boost::math::quantile(dist, static_cast<RealT>(0.8609299926107957)), static_cast<RealT>(1.9999999999999947), tol );
BOOST_CHECK_CLOSE( boost::math::quantile(dist, static_cast<RealT>(0.9348833491908337)), static_cast<RealT>(3), tol );
BOOST_CHECK_CLOSE( boost::math::quantile(dist, static_cast<RealT>(0.966198875597724)), static_cast<RealT>(3.9999999999999964), tol );
}
template <typename RealT>
void test_ccdf()
{
const RealT tol = make_tolerance<RealT>();
const RealT probs[] = {0.2, 0.3, 0.5};
const RealT rates[] = {0.5, 1.0, 1.5};
const std::size_t n = sizeof(probs)/sizeof(RealT);
boost::math::hyperexponential_distribution<RealT> dist(probs, probs+n, rates, rates+n);
// Mathematica: Table[SurvivalFunction[HyperexponentialDistribution[{0.2, 0.3, 0.5}, {.5, 1.0, 1.5}], x], {x, 0, 4}]
BOOST_CHECK_CLOSE( boost::math::cdf(boost::math::complement(dist, static_cast<RealT>(0))), static_cast<RealT>(1), tol );
BOOST_CHECK_CLOSE( boost::math::cdf(boost::math::complement(dist, static_cast<RealT>(1))), static_cast<RealT>(0.3432350443681743), tol );
BOOST_CHECK_CLOSE( boost::math::cdf(boost::math::complement(dist, static_cast<RealT>(2))), static_cast<RealT>(0.13907000738920425), tol );
BOOST_CHECK_CLOSE( boost::math::cdf(boost::math::complement(dist, static_cast<RealT>(3))), static_cast<RealT>(0.0651166508091663), tol );
BOOST_CHECK_CLOSE( boost::math::cdf(boost::math::complement(dist, static_cast<RealT>(4))), static_cast<RealT>(0.03380112440227598), tol );
}
template <typename RealT>
void test_cquantile()
{
const RealT tol = make_tolerance<RealT>();
const RealT probs[] = {0.2, 0.3, 0.5};
const RealT rates[] = {0.5, 1.0, 1.5};
const std::size_t n = sizeof(probs)/sizeof(RealT);
boost::math::hyperexponential_distribution<RealT> dist(probs, probs+n, rates, rates+n);
// Mathematica: Table[SurvivalFunction[HyperexponentialDistribution[{0.2, 0.3, 0.5}, {.5, 1.0, 1.5}], p], {p, {1., 0.3432350443681743, 0.13907000738920425, 0.0651166508091663, 0.03380112440227598}}]
BOOST_CHECK_CLOSE( boost::math::quantile(boost::math::complement(dist, static_cast<RealT>(1))), static_cast<RealT>(0), tol );
BOOST_CHECK_CLOSE( boost::math::quantile(boost::math::complement(dist, static_cast<RealT>(0.3432350443681743))), static_cast<RealT>(1.0000000000000036), tol );
BOOST_CHECK_CLOSE( boost::math::quantile(boost::math::complement(dist, static_cast<RealT>(0.13907000738920425))), static_cast<RealT>(1.9999999999999947), tol );
BOOST_CHECK_CLOSE( boost::math::quantile(boost::math::complement(dist, static_cast<RealT>(0.0651166508091663))), static_cast<RealT>(3), tol );
BOOST_CHECK_CLOSE( boost::math::quantile(boost::math::complement(dist, static_cast<RealT>(0.03380112440227598))), static_cast<RealT>(3.9999999999999964), tol );
}
template <typename RealT>
void test_mean()
{
const RealT tol = make_tolerance<RealT>();
const RealT probs[] = {0.2, 0.3, 0.5};
const RealT rates[] = {0.5, 1.0, 1.5};
const std::size_t n = sizeof(probs)/sizeof(RealT);
boost::math::hyperexponential_distribution<RealT> dist(probs, probs+n, rates, rates+n);
// Mathematica: Mean[HyperexponentialDistribution[{0.2, 0.3, 0.5}, {.5, 1.0, 1.5}]]
BOOST_CHECK_CLOSE( boost::math::mean(dist), static_cast<RealT>(1.0333333333333332), tol );
}
template <typename RealT>
void test_variance()
{
const RealT tol = make_tolerance<RealT>();
const RealT probs[] = {0.2, 0.3, 0.5};
const RealT rates[] = {0.5, 1.0, 1.5};
const std::size_t n = sizeof(probs)/sizeof(RealT);
boost::math::hyperexponential_distribution<RealT> dist(probs, probs+n, rates, rates+n);
// Mathematica: Mean[HyperexponentialDistribution[{0.2, 0.3, 0.5}, {.5, 1.0, 1.5}]]
BOOST_CHECK_CLOSE( boost::math::variance(dist), static_cast<RealT>(1.5766666666666673), tol );
}
template <typename RealT>
void test_kurtosis()
{
const RealT tol = make_tolerance<RealT>();
const RealT probs[] = {0.2, 0.3, 0.5};
const RealT rates[] = {0.5, 1.0, 1.5};
const std::size_t n = sizeof(probs)/sizeof(RealT);
boost::math::hyperexponential_distribution<RealT> dist(probs, probs+n, rates, rates+n);
// Mathematica: Kurtosis[HyperexponentialDistribution[{0.2, 0.3, 0.5}, {.5, 1.0, 1.5}]]
BOOST_CHECK_CLOSE( boost::math::kurtosis(dist), static_cast<RealT>(19.75073861680871), tol );
BOOST_CHECK_CLOSE( boost::math::kurtosis_excess(dist), static_cast<RealT>(19.75073861680871)-static_cast<RealT>(3), tol );
}
template <typename RealT>
void test_skewness()
{
const RealT tol = make_tolerance<RealT>();
const RealT probs[] = {0.2, 0.3, 0.5};
const RealT rates[] = {0.5, 1.0, 1.5};
const std::size_t n = sizeof(probs)/sizeof(RealT);
boost::math::hyperexponential_distribution<RealT> dist(probs, probs+n, rates, rates+n);
// Mathematica: Skewness[HyperexponentialDistribution[{0.2, 0.3, 0.5}, {.5, 1.0, 1.5}]]
BOOST_CHECK_CLOSE( boost::math::skewness(dist), static_cast<RealT>(3.181138744996378), tol );
}
template <typename RealT>
void test_mode()
{
const RealT tol = make_tolerance<RealT>();
const RealT probs[] = {0.2, 0.3, 0.5};
const RealT rates[] = {0.5, 1.0, 1.5};
const std::size_t n = sizeof(probs)/sizeof(RealT);
boost::math::hyperexponential_distribution<RealT> dist(probs, probs+n, rates, rates+n);
BOOST_CHECK_CLOSE( boost::math::mode(dist), static_cast<RealT>(0), tol );
}
BOOST_AUTO_TEST_CASE( range )
{
test_range<float>();
test_range<double>();
//test_range<boost::math::concepts::real_concept>();
}
BOOST_AUTO_TEST_CASE( support )
{
test_support<float>();
test_support<double>();
}
BOOST_AUTO_TEST_CASE( pdf )
{
test_pdf<float>();
test_pdf<double>();
}
BOOST_AUTO_TEST_CASE( cdf )
{
test_cdf<float>();
test_cdf<double>();
}
BOOST_AUTO_TEST_CASE( quantile )
{
test_quantile<float>();
test_quantile<double>();
}
BOOST_AUTO_TEST_CASE( ccdf )
{
test_ccdf<float>();
test_ccdf<double>();
}
BOOST_AUTO_TEST_CASE( cquantile )
{
test_cquantile<float>();
test_cquantile<double>();
}
BOOST_AUTO_TEST_CASE( mean )
{
test_mean<float>();
test_mean<double>();
}
BOOST_AUTO_TEST_CASE( variance )
{
test_variance<float>();
test_variance<double>();
}
BOOST_AUTO_TEST_CASE( kurtosis )
{
test_kurtosis<float>();
test_kurtosis<double>();
}
BOOST_AUTO_TEST_CASE( skewness )
{
test_skewness<float>();
test_skewness<double>();
}
BOOST_AUTO_TEST_CASE( mode )
{
test_mode<float>();
test_mode<double>();
}
Reference in New Issue View Git Blame Copy Permalink