boost_mirror/math
1
0
Fork 0
mirror of https://github.com/boostorg/math.git synced 2025-05-11 21:33:52 +00:00
Code Issues Packages Projects Releases Wiki Activity
math/test/test_geometric.cpp
Matt Borland d01893d215
Add GPU markers to fisher f dist 
Add SYCL testing of fisher f dist

Add CUDA fisher f dist testing

Add NVRTC fisher f dist testing

Add GPU support to gamma dist

Add SYCL testing of gamma dist

Add CUDA gamma dist testing

Add NVRTC gamma dist testing

Reduce number of threads per block since it can crash CI

Add GPU support to the geometric dist

Add SYCL testing of geometric dist

Add cuda::std::tie

Add GPU support to inv_discrete_quantile

Add CUDA testing of geometric dist

Add NVRTC testing of geometric dist

Add SYCL testing of inverse_chi_squared dist

Adjust tol

Add NVRTC inverse chi squared dist testing

Add CUDA inverse chi squared dist testing

Add GPU support to inverse gamma dist

Add SYCL testing to inverse gamma dist

Add NVRTC testing of inverse gamma dist

Add CUDA testing of inverse gamma dist
2024-09-04 11:07:17 -04:00

851 lines
33 KiB
C++
Raw Permalink Blame History

// test_geometric.cpp
// Copyright Paul A. Bristow 2010.
// Copyright John Maddock 2010.
// 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)
// Tests for Geometric Distribution.
// Note that these defines must be placed BEFORE #includes.
#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
// because several tests overflow & underflow by design.
#define BOOST_MATH_DISCRETE_QUANTILE_POLICY real
#ifdef _MSC_VER
# pragma warning(disable: 4127) // conditional expression is constant.
#endif
#if !defined(TEST_FLOAT) && !defined(TEST_DOUBLE) && !defined(TEST_LDOUBLE) && !defined(TEST_REAL_CONCEPT)
# define TEST_FLOAT
# define TEST_DOUBLE
# define TEST_LDOUBLE
# define TEST_REAL_CONCEPT
#endif
#include <boost/math/tools/config.hpp>
#include "../include_private/boost/math/tools/test.hpp"
#ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
#include <boost/math/concepts/real_concept.hpp> // for real_concept
using ::boost::math::concepts::real_concept;
#endif
#include <boost/math/distributions/geometric.hpp> // for geometric_distribution
using boost::math::geometric_distribution;
using boost::math::geometric; // using typedef for geometric_distribution<double>
#include <boost/math/distributions/negative_binomial.hpp> // for some comparisons.
#define BOOST_TEST_MAIN
#include <boost/test/unit_test.hpp> // for test_main
#include <boost/test/tools/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE_FRACTION
#include "test_out_of_range.hpp"
#include <iostream>
using std::cout;
using std::endl;
using std::setprecision;
using std::showpoint;
#include <limits>
using std::numeric_limits;
#include <cmath>
using std::log;
using std::abs;
#include <type_traits>
template <class RealType>
void test_spot( // Test a single spot value against 'known good' values.
RealType k, // Number of failures.
RealType p, // Probability of success_fraction.
RealType P, // CDF probability.
RealType Q, // Complement of CDF.
RealType logP, // Logcdf probability
RealType logQ, // Complement of logcdf
RealType tol, // Test tolerance
RealType logtol) // Logcdf Test tolerance.
{
BOOST_IF_CONSTEXPR (std::is_same<RealType, long double>::value
#ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
|| std::is_same<RealType, real_concept>::value
#endif
)
{
logtol *= 100;
}
boost::math::geometric_distribution<RealType> g(p);
BOOST_CHECK_EQUAL(p, g.success_fraction());
BOOST_CHECK_CLOSE_FRACTION(cdf(g, k), P, tol);
BOOST_CHECK_CLOSE_FRACTION(logcdf(g, k), logP, logtol);
if((P < 0.99) && (Q < 0.99))
{
// We can only check this if P is not too close to 1,
// so that we can guarantee that Q is free of error:
//
BOOST_CHECK_CLOSE_FRACTION(
cdf(complement(g, k)), Q, tol);
BOOST_CHECK_CLOSE_FRACTION(
logcdf(complement(g, k)), logQ, logtol);
if(k != 0)
{
BOOST_CHECK_CLOSE_FRACTION(
quantile(g, P), k, tol);
}
else
{
// Just check quantile is very small:
if((std::numeric_limits<RealType>::max_exponent <= std::numeric_limits<double>::max_exponent)
&& (boost::is_floating_point<RealType>::value))
{
// Limit where this is checked: if exponent range is very large we may
// run out of iterations in our root finding algorithm.
BOOST_CHECK(quantile(g, P) < boost::math::tools::epsilon<RealType>() * 10);
}
}
if(k != 0)
{
BOOST_CHECK_CLOSE_FRACTION(
quantile(complement(g, Q)), k, tol);
}
else
{
// Just check quantile is very small:
if((std::numeric_limits<RealType>::max_exponent <= std::numeric_limits<double>::max_exponent)
&& (boost::is_floating_point<RealType>::value))
{
// Limit where this is checked: if exponent range is very large we may
// run out of iterations in our root finding algorithm.
BOOST_CHECK(quantile(complement(g, Q)) < boost::math::tools::epsilon<RealType>() * 10);
}
}
} // if((P < 0.99) && (Q < 0.99))
// Parameter estimation test: estimate success ratio:
BOOST_CHECK_CLOSE_FRACTION(
geometric_distribution<RealType>::find_lower_bound_on_p(
1+k, P),
p, 0.02); // Wide tolerance needed for some tests.
// Note we bump up the sample size here, purely for the sake of the test,
// internally the function has to adjust the sample size so that we get
// the right upper bound, our test undoes this, so we can verify the result.
BOOST_CHECK_CLOSE_FRACTION(
geometric_distribution<RealType>::find_upper_bound_on_p(
1+k+1, Q),
p, 0.02);
if(Q < P)
{
//
// We check two things here, that the upper and lower bounds
// are the right way around, and that they do actually bracket
// the naive estimate of p = successes / (sample size)
//
BOOST_CHECK(
geometric_distribution<RealType>::find_lower_bound_on_p(
1+k, Q)
<=
geometric_distribution<RealType>::find_upper_bound_on_p(
1+k, Q)
);
BOOST_CHECK(
geometric_distribution<RealType>::find_lower_bound_on_p(
1+k, Q)
<=
1 / (1+k)
);
BOOST_CHECK(
1 / (1+k)
<=
geometric_distribution<RealType>::find_upper_bound_on_p(
1+k, Q)
);
}
else
{
// As above but when P is small.
BOOST_CHECK(
geometric_distribution<RealType>::find_lower_bound_on_p(
1+k, P)
<=
geometric_distribution<RealType>::find_upper_bound_on_p(
1+k, P)
);
BOOST_CHECK(
geometric_distribution<RealType>::find_lower_bound_on_p(
1+k, P)
<=
1 / (1+k)
);
BOOST_CHECK(
1 / (1+k)
<=
geometric_distribution<RealType>::find_upper_bound_on_p(
1+k, P)
);
}
// Estimate sample size:
BOOST_CHECK_CLOSE_FRACTION(
geometric_distribution<RealType>::find_minimum_number_of_trials(
k, p, P),
1+k, 0.02); // Can differ 50 to 51 for small p
BOOST_CHECK_CLOSE_FRACTION(
geometric_distribution<RealType>::find_maximum_number_of_trials(
k, p, Q),
1+k, 0.02);
} // test_spot
template <class RealType> // Any floating-point type RealType.
void test_spots(RealType)
{
// Basic sanity checks.
// Most test data is to double precision (17 decimal digits) only,
cout << "Floating point Type is " << typeid(RealType).name() << endl;
// so set tolerance to 1000 eps expressed as a fraction,
// or 1000 eps of type double expressed as a fraction,
// whichever is the larger.
RealType tolerance = (std::max)
(boost::math::tools::epsilon<RealType>(),
static_cast<RealType>(std::numeric_limits<double>::epsilon()));
tolerance *= 10; // 10 eps
cout << "Tolerance = " << tolerance << "." << endl;
RealType tol1eps = boost::math::tools::epsilon<RealType>(); // Very tight, suit exact values.
//RealType tol2eps = boost::math::tools::epsilon<RealType>() * 2; // Tight, values.
RealType tol5eps = boost::math::tools::epsilon<RealType>() * 5; // Wider 5 epsilon.
cout << "Tolerance 5 eps = " << tol5eps << "." << endl;
// Sources of spot test values are mainly R.
using boost::math::geometric_distribution;
using boost::math::geometric;
using boost::math::cdf;
using boost::math::pdf;
using boost::math::quantile;
using boost::math::complement;
BOOST_MATH_STD_USING // for std math functions
// Test geometric using cdf spot values R
// These test quantiles and complements as well.
test_spot( //
static_cast<RealType>(2), // Number of failures, k
static_cast<RealType>(0.5), // Probability of success as fraction, p
static_cast<RealType>(0.875L), // Probability of result (CDF), P
static_cast<RealType>(0.125L), // complement CCDF Q = 1 - P
static_cast<RealType>(-0.1335313926245226231463436209313499745894L),
static_cast<RealType>(-2.079441541679835928251696364374529704227L),
tolerance,
tolerance);
test_spot( //
static_cast<RealType>(0), // Number of failures, k
static_cast<RealType>(0.25), // Probability of success as fraction, p
static_cast<RealType>(0.25), // Probability of result (CDF), P
static_cast<RealType>(0.75), // Q = 1 - P
static_cast<RealType>(-1.386294361119890618834464242916353136151L),
static_cast<RealType>(-0.2876820724517809274392190059938274315035L),
tolerance,
tolerance);
test_spot(
// R formatC(pgeom(10,0.25), digits=17) [1] "0.95776486396789551"
// formatC(pgeom(10,0.25, FALSE), digits=17) [1] "0.042235136032104499"
static_cast<RealType>(10), // Number of failures, k
static_cast<RealType>(0.25), // Probability of success, p
static_cast<RealType>(0.95776486396789551L), // Probability of result (CDF), P
static_cast<RealType>(0.042235136032104499L), // Q = 1 - P
static_cast<RealType>(-0.04315297584768019483875419429616349387993L),
static_cast<RealType>(-3.164502796969590201831409065932101746539L),
tolerance,
tolerance);
test_spot( //
// > R formatC(pgeom(50,0.25, TRUE), digits=17) [1] "0.99999957525875771"
// > R formatC(pgeom(50,0.25, FALSE), digits=17) [1] "4.2474124232020353e-07"
static_cast<RealType>(50), // Number of failures, k
static_cast<RealType>(0.25), // Probability of success, p
static_cast<RealType>(0.99999957525875771), // Probability of result (CDF), P
static_cast<RealType>(4.2474124232020353e-07), // Q = 1 - P
static_cast<RealType>(-4.247413325227902241937783772756893037512e-7L),
static_cast<RealType>(-14.67178569504082729940016930568519898711L),
tolerance,
tolerance);
/*
// This causes failures in find_upper_bound_on_p p is small branch.
test_spot( // formatC(pgeom(50,0.01, TRUE), digits=17)[1] "0.40104399353383874"
// > formatC(pgeom(50,0.01, FALSE), digits=17) [1] "0.59895600646616121"
static_cast<RealType>(50), // Number of failures, k
static_cast<RealType>(0.01), // Probability of success, p
static_cast<RealType>(0.40104399353383874), // Probability of result (CDF), P
static_cast<RealType>(0.59895600646616121), // Q = 1 - P
tolerance);
*/
test_spot( // > formatC(pgeom(50,0.99, TRUE), digits=17) [1] " 1"
// formatC(pgeom(50,0.99, FALSE), digits=17) [1] "1.0000000000000364e-102"
static_cast<RealType>(50), // Number of failures, k
static_cast<RealType>(0.99), // Probability of success, p
static_cast<RealType>(1), // Probability of result (CDF), P
static_cast<RealType>(1.0000000000000364e-102), // Q = 1 - P
static_cast<RealType>(-1.0000000000000364e-102L),
static_cast<RealType>(-std::numeric_limits<RealType>::infinity()),
tolerance,
tolerance * 100);
test_spot( // > formatC(pgeom(1,0.99, TRUE), digits=17) [1] "0.99990000000000001"
// > formatC(pgeom(1,0.99, FALSE), digits=17) [1] "0.00010000000000000009"
static_cast<RealType>(1), // Number of failures, k
static_cast<RealType>(0.99), // Probability of success, p
static_cast<RealType>(0.9999), // Probability of result (CDF), P
static_cast<RealType>(0.0001), // Q = 1 - P
static_cast<RealType>(-0.0001000050003333583353335000142869643968354L),
static_cast<RealType>(-9.210340371976182736071965818737456830404L),
tolerance,
tolerance * 100);
if(std::numeric_limits<RealType>::is_specialized)
{ // An extreme value test that is more accurate than using negative binomial.
// Since geometric only uses exp and log functions.
test_spot( // > formatC(pgeom(10000, 0.001, TRUE), digits=17) [1] "0.99995487182736897"
// > formatC(pgeom(10000,0.001, FALSE), digits=17) [1] "4.5128172631071587e-05"
static_cast<RealType>(10000L), // Number of failures, k
static_cast<RealType>(0.001L), // Probability of success, p
static_cast<RealType>(0.99995487182736897L), // Probability of result (CDF), P
static_cast<RealType>(4.5128172631071587e-05L), // Q = 1 - P,
static_cast<RealType>(-0.00004512919093769043386238651458397312570531L),
static_cast<RealType>(-10.00600383616891853492996552293751795172L),
tolerance,
tolerance * 100); //
} // numeric_limit is specialized
// End of single spot tests using RealType
// Tests on PDF:
BOOST_CHECK_CLOSE_FRACTION( //> formatC(dgeom(0,0.5), digits=17)[1] " 0.5"
pdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)),
static_cast<RealType>(0.0) ), // Number of failures, k is very small but not integral,
static_cast<RealType>(0.5), // nearly success probability.
tolerance);
BOOST_CHECK_CLOSE_FRACTION( //> formatC(dgeom(0,0.5), digits=17)[1] " 0.5"
// R treats geom as a discrete distribution.
// > formatC(dgeom(1.999999,0.5, FALSE), digits=17) [1] " 0"
// Warning message:
// In dgeom(1.999999, 0.5, FALSE) : non-integer x = 1.999999
pdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)),
static_cast<RealType>(0.0001L) ), // Number of failures, k is very small but not integral,
static_cast<RealType>(0.4999653438420768L), // nearly success probability.
tolerance);
BOOST_CHECK_CLOSE_FRACTION( // > formatC(pgeom(0.0001,0.5, TRUE), digits=17)[1] " 0.5"
// > formatC(pgeom(0.0001,0.5, FALSE), digits=17) [1] " 0.5"
// R treats geom as a discrete distribution.
pdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)),
static_cast<RealType>(0.0001L) ), // Number of failures, k is very small but not integral,
static_cast<RealType>(0.4999653438420768L), // nearly success probability.
tolerance);
BOOST_CHECK_CLOSE_FRACTION( // formatC(dgeom(1,0.01), digits=17)[1] "0.0099000000000000008"
pdf(geometric_distribution<RealType>(static_cast<RealType>(0.01L)),
static_cast<RealType>(1) ), // Number of failures, k
static_cast<RealType>(0.0099000000000000008), //
tolerance);
BOOST_CHECK_CLOSE_FRACTION( //> formatC(dgeom(1,0.99), digits=17)[1] "0.0099000000000000043"
pdf(geometric_distribution<RealType>(static_cast<RealType>(0.99L)),
static_cast<RealType>(1) ), // Number of failures, k
static_cast<RealType>(0.00990000000000000043L), //
tolerance);
BOOST_CHECK_CLOSE_FRACTION( //> > formatC(dgeom(0,0.99), digits=17)[1] "0.98999999999999999"
pdf(geometric_distribution<RealType>(static_cast<RealType>(0.99L)),
static_cast<RealType>(0) ), // Number of failures, k
static_cast<RealType>(0.98999999999999999L), //
tolerance);
// p near unity.
BOOST_CHECK_CLOSE_FRACTION( // > formatC(dgeom(100,0.99), digits=17)[1] "9.9000000000003448e-201"
pdf(geometric_distribution<RealType>(static_cast<RealType>(0.99L)),
static_cast<RealType>(100) ), // Number of failures, k
static_cast<RealType>(9.9000000000003448e-201L), //
100 * tolerance); // Note difference
// p nearer unity.
// On GPU this gets flushed to 0 which has an eps difference of 3.4e+38
#ifndef BOOST_MATH_HAS_GPU_SUPPORT
BOOST_CHECK_CLOSE_FRACTION( //
pdf(geometric_distribution<RealType>(static_cast<RealType>(0.9999)),
static_cast<RealType>(10) ), // Number of failures, k
// static_cast<double>(9.9989999999889024e-41), // Boost.Math
// static_cast<float>(1.00156406e-040)
static_cast<RealType>(9.999e-41), // exact from 100 digit calculator.
2e3 * tolerance); // Note bigger tolerance needed.
#endif
// Moshier Cephes 100 digits calculator says 9.999e-41
//0.9999*pow(1-0.9999,10)
// 9.9990000000000000000000000000000000000000000000000000000000000000000000E-41
// 9.998999999988988e-041
// > formatC(dgeom(10, 0.9999), digits=17) [1] "9.9989999999889024e-41"
// p * pow(q, k) 9.9989999999889880e-041
// exp(p * k * log1p(-p)) 9.9989999999889024e-041
// 0.9999999999 * pow(1-0.9999999999,10)= 9.9999999990E-101
// > formatC(dgeom(10,0.9999999999), digits=17) [1] "1.0000008273040127e-100"
BOOST_CHECK_CLOSE_FRACTION( //
pdf(geometric_distribution<RealType>(static_cast<RealType>(0.9999999999L)),
static_cast<RealType>(10) ), //
static_cast<RealType>(9.9999999990E-101L), // 1.0000008273040179e-100
1e9 * tolerance); // Note big tolerance needed.
// 1.0000008273040179e-100 Boost.Math
// 1.0000008273040127e-100 R
// 0.9999999990000004e-100 100 digit calculator 'exact'
BOOST_CHECK_CLOSE_FRACTION( //
pdf(geometric_distribution<RealType>(static_cast<RealType>(0.00000000001L)),
static_cast<RealType>(10) ), //
static_cast<RealType>(9.999999999e-12L), // get 9.9999999989999994e-012
1 * tolerance); // Note small tolerance needed.
BOOST_CHECK_CLOSE_FRACTION( //
pdf(geometric_distribution<RealType>(static_cast<RealType>(0.00000000001L)),
static_cast<RealType>(1000) ), //
static_cast<RealType>(9.9999999e-12L), // get 9.9999998999999913e-012
tolerance); // Note small tolerance needed.
///////////////////////////////////////////////////
BOOST_CHECK_CLOSE_FRACTION( //
// > formatC(dgeom(0.0001,0.5, FALSE), digits=17) [1] " 0.5"
// R treats geom as a discrete distribution.
// But Boost.Math is continuous, so if you want R behaviour,
// make number of failures, k into an integer with the floor function.
pdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)),
static_cast<RealType>(floor(0.0001L)) ), // Number of failures, k is very small but MADE integral,
static_cast<RealType>(0.5), // nearly success probability.
tolerance);
// R switches over at about 1e7 from k = 0, returning 0.5, to k = 1, returning 0.25.
// Boost.Math does not do this, even for 0.9999999999999999
// > formatC(pgeom(0.999999,0.5, FALSE), digits=17) [1] " 0.5"
// > formatC(pgeom(0.9999999,0.5, FALSE), digits=17) [1] " 0.25"
BOOST_CHECK_CLOSE_FRACTION( // > formatC(pgeom(0.0001,0.5, TRUE), digits=17)[1] " 0.5"
// > formatC(pgeom(0.0001,0.5, FALSE), digits=17) [1] " 0.5"
// R treats geom as a discrete distribution.
// But Boost.Math is continuous, so if you want R behaviour,
// make number of failures, k into an integer with the floor function.
pdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)),
static_cast<RealType>(floor(0.9999999999999999L)) ), // Number of failures, k is very small but MADE integral,
static_cast<RealType>(0.5), // nearly success probability.
tolerance);
BOOST_CHECK_CLOSE_FRACTION( // > formatC(pgeom(0.0001,0.5, TRUE), digits=17)[1] " 0.5"
// > formatC(pgeom(0.0001,0.5, FALSE), digits=17) [1] " 0.5"
// R treats geom as a discrete distribution.
// But Boost.Math is continuous, so if you want R behaviour,
// make number of failures, k into an integer with the floor function.
pdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)),
static_cast<RealType>(floor(1. - tolerance)) ),
// Number of failures, k is very small but MADE integral,
// Need to use tolerance here,
// as epsilon is ill-defined for Real concept:
// numeric_limits<RealType>::epsilon() 0
static_cast<RealType>(0.5), // nearly success probability.
tolerance * 10);
BOOST_CHECK_CLOSE_FRACTION(
pdf(geometric_distribution<RealType>(static_cast<RealType>(0.0001L)),
static_cast<RealType>(2)), // k = 2.
static_cast<RealType>(9.99800010e-5L), // 'exact '
tolerance);
//> formatC(dgeom(2, 0.9999), digits=17) [1] "9.9989999999977806e-09"
BOOST_CHECK_CLOSE_FRACTION(
pdf(geometric_distribution<RealType>(static_cast<RealType>(0.9999L)),
static_cast<RealType>(2)), // k = 0
static_cast<RealType>(9.999e-9L), // 'exact'
1000*tolerance);
BOOST_CHECK_CLOSE_FRACTION(
pdf(geometric_distribution<RealType>(static_cast<RealType>(0.9999L)),
static_cast<RealType>(3)), // k = 3
static_cast<RealType>(9.999e-13L), // get
1000*tolerance);
BOOST_CHECK_CLOSE_FRACTION(
pdf(geometric_distribution<RealType>(static_cast<RealType>(0.9999L)),
static_cast<RealType>(5)), // k = 5
static_cast<RealType>(9.999e-21L), // 9.9989999999944947e-021
1000*tolerance);
BOOST_CHECK_CLOSE_FRACTION(
pdf(geometric_distribution<RealType>( static_cast<RealType>(0.0001L)),
static_cast<RealType>(3)), // k = 0.
static_cast<RealType>(9.99700029999e-5L), //
tolerance);
// Tests on cdf:
// MathCAD pgeom k, r, p) == failures, successes, probability.
BOOST_CHECK_CLOSE_FRACTION(cdf(
geometric_distribution<RealType>(static_cast<RealType>(0.5)), // prob 0.5
static_cast<RealType>(0) ), // k = 0
static_cast<RealType>(0.5), // probability =p
tolerance);
BOOST_CHECK_CLOSE_FRACTION(cdf(complement(
geometric_distribution<RealType>(static_cast<RealType>(0.5)), //
static_cast<RealType>(0) )), // k = 0
static_cast<RealType>(0.5), // probability =
tolerance);
BOOST_CHECK_CLOSE_FRACTION(cdf(
geometric_distribution<RealType>(static_cast<RealType>(0.25)), // prob 0.5
static_cast<RealType>(1) ), // k = 0
static_cast<RealType>(0.4375L), // probability =p
tolerance);
BOOST_CHECK_CLOSE_FRACTION(cdf(complement(
geometric_distribution<RealType>(static_cast<RealType>(0.25)), //
static_cast<RealType>(1) )), // k = 0
static_cast<RealType>(1-0.4375L), // probability =
tolerance);
BOOST_CHECK_CLOSE_FRACTION(cdf(complement(
geometric_distribution<RealType>(static_cast<RealType>(0.5)), //
static_cast<RealType>(1) )), // k = 0
static_cast<RealType>(0.25), // probability = exact 0.25
tolerance);
BOOST_CHECK_CLOSE_FRACTION( //
cdf(geometric_distribution<RealType>(static_cast<RealType>(0.5)),
static_cast<RealType>(4)), // k =4.
static_cast<RealType>(0.96875L), // exact
tolerance);
// Tests of other functions, mean and other moments ...
geometric_distribution<RealType> dist(static_cast<RealType>(0.25));
// mean:
BOOST_CHECK_CLOSE_FRACTION(
mean(dist), static_cast<RealType>((1 - 0.25) /0.25), tol5eps);
BOOST_CHECK_CLOSE_FRACTION(
mode(dist), static_cast<RealType>(0), tol1eps);
// variance:
BOOST_CHECK_CLOSE_FRACTION(
variance(dist), static_cast<RealType>((1 - 0.25) / (0.25 * 0.25)), tol5eps);
// std deviation:
// sqrt(0.75/0.125)
BOOST_CHECK_CLOSE_FRACTION(
standard_deviation(dist), //
static_cast<RealType>(sqrt((1.0L - 0.25L) / (0.25L * 0.25L))), // using 100 digit calc
tol5eps);
BOOST_CHECK_CLOSE_FRACTION(
skewness(dist), //
static_cast<RealType>((2-0.25L) /sqrt(0.75L)),
// using calculator
tol5eps);
BOOST_CHECK_CLOSE_FRACTION(
kurtosis_excess(dist), //
static_cast<RealType>(6 + 0.0625L/0.75L), //
tol5eps);
// 6.083333333333333 6.166666666666667
BOOST_CHECK_CLOSE_FRACTION(
kurtosis(dist), // true
static_cast<RealType>(9 + 0.0625L/0.75L), //
tol5eps);
// hazard:
RealType x = static_cast<RealType>(0.125);
BOOST_CHECK_CLOSE_FRACTION(
hazard(dist, x)
, pdf(dist, x) / cdf(complement(dist, x)), tol5eps);
// cumulative hazard:
BOOST_CHECK_CLOSE_FRACTION(
chf(dist, x), -log(cdf(complement(dist, x))), tol5eps);
// coefficient_of_variation:
BOOST_CHECK_CLOSE_FRACTION(
coefficient_of_variation(dist)
, standard_deviation(dist) / mean(dist), tol5eps);
// Special cases for PDF:
BOOST_CHECK_EQUAL(
pdf(
geometric_distribution<RealType>(static_cast<RealType>(0)), //
static_cast<RealType>(0)),
static_cast<RealType>(0) );
BOOST_CHECK_EQUAL(
pdf(
geometric_distribution<RealType>(static_cast<RealType>(0)),
static_cast<RealType>(0.0001)),
static_cast<RealType>(0) );
BOOST_CHECK_EQUAL(
pdf(
geometric_distribution<RealType>(static_cast<RealType>(1)),
static_cast<RealType>(0.001)),
static_cast<RealType>(0) );
BOOST_CHECK_EQUAL(
pdf(
geometric_distribution<RealType>(static_cast<RealType>(1)),
static_cast<RealType>(8)),
static_cast<RealType>(0) );
BOOST_CHECK_SMALL(
pdf(
geometric_distribution<RealType>(static_cast<RealType>(0.25)),
static_cast<RealType>(0))-
static_cast<RealType>(0.25),
2 * boost::math::tools::epsilon<RealType>() ); // Expect exact, but not quite.
// numeric_limits<RealType>::epsilon()); // Not suitable for real concept!
// Quantile boundary cases checks:
BOOST_CHECK_EQUAL(
quantile( // zero P < cdf(0) so should be exactly zero.
geometric_distribution<RealType>(static_cast<RealType>(0.25)),
static_cast<RealType>(0)),
static_cast<RealType>(0));
BOOST_CHECK_EQUAL(
quantile( // min P < cdf(0) so should be exactly zero.
geometric_distribution<RealType>(static_cast<RealType>(0.25)),
static_cast<RealType>(boost::math::tools::min_value<RealType>())),
static_cast<RealType>(0));
BOOST_CHECK_CLOSE_FRACTION(
quantile( // Small P < cdf(0) so should be near zero.
geometric_distribution<RealType>(static_cast<RealType>(0.25)),
static_cast<RealType>(boost::math::tools::epsilon<RealType>())), //
static_cast<RealType>(0),
tol5eps);
BOOST_CHECK_CLOSE_FRACTION(
quantile( // Small P < cdf(0) so should be exactly zero.
geometric_distribution<RealType>(static_cast<RealType>(0.25)),
static_cast<RealType>(0.0001)),
static_cast<RealType>(0),
tolerance);
//BOOST_CHECK( // Fails with overflow for real_concept
//quantile( // Small P near 1 so k failures should be big.
//geometric_distribution<RealType>(static_cast<RealType>(8), static_cast<RealType>(0.25)),
//static_cast<RealType>(1 - boost::math::tools::epsilon<RealType>())) <=
//static_cast<RealType>(189.56999032670058) // 106.462769 for float
//);
if(std::numeric_limits<RealType>::has_infinity)
{ // BOOST_CHECK tests for infinity using std::numeric_limits<>::infinity()
// Note that infinity is not implemented for real_concept, so these tests
// are only done for types, like built-in float, double.. that have infinity.
// Note that these assume that BOOST_MATH_OVERFLOW_ERROR_POLICY is NOT throw_on_error.
// #define BOOST_MATH_THROW_ON_OVERFLOW_POLICY == throw_on_error would throw here.
// #define BOOST_MAT_DOMAIN_ERROR_POLICY IS defined throw_on_error,
// so the throw path of error handling is tested below with BOOST_MATH_CHECK_THROW tests.
BOOST_CHECK(
quantile( // At P == 1 so k failures should be infinite.
geometric_distribution<RealType>(static_cast<RealType>(0.25)),
static_cast<RealType>(1)) ==
//static_cast<RealType>(boost::math::tools::infinity<RealType>())
static_cast<RealType>(std::numeric_limits<RealType>::infinity()) );
BOOST_CHECK_EQUAL(
quantile( // At 1 == P so should be infinite.
geometric_distribution<RealType>( static_cast<RealType>(0.25)),
static_cast<RealType>(1)), //
std::numeric_limits<RealType>::infinity() );
BOOST_CHECK_EQUAL(
quantile(complement( // Q zero 1 so P == 1 < cdf(0) so should be exactly infinity.
geometric_distribution<RealType>(static_cast<RealType>(0.25)),
static_cast<RealType>(0))),
std::numeric_limits<RealType>::infinity() );
} // test for infinity using std::numeric_limits<>::infinity()
else
{ // real_concept case, so check it throws rather than returning infinity.
BOOST_CHECK_EQUAL(
quantile( // At P == 1 so k failures should be infinite.
geometric_distribution<RealType>(static_cast<RealType>(0.25)),
static_cast<RealType>(1)),
boost::math::tools::max_value<RealType>() );
BOOST_CHECK_EQUAL(
quantile(complement( // Q zero 1 so P == 1 < cdf(0) so should be exactly infinity.
geometric_distribution<RealType>(static_cast<RealType>(0.25)),
static_cast<RealType>(0))),
boost::math::tools::max_value<RealType>());
} // has infinity
BOOST_CHECK( // Should work for built-in and real_concept.
quantile(complement( // Q near to 1 so P nearly 1, so should be large > 300.
geometric_distribution<RealType>(static_cast<RealType>(0.25)),
static_cast<RealType>(boost::math::tools::min_value<RealType>())))
>= static_cast<RealType>(300) );
BOOST_CHECK_EQUAL(
quantile( // P == 0 < cdf(0) so should be zero.
geometric_distribution<RealType>(static_cast<RealType>(0.25)),
static_cast<RealType>(0)),
static_cast<RealType>(0));
// Quantile Complement boundary cases:
BOOST_CHECK_EQUAL(
quantile(complement( // Q = 1 so P = 0 < cdf(0) so should be exactly zero.
geometric_distribution<RealType>( static_cast<RealType>(0.25)),
static_cast<RealType>(1))),
static_cast<RealType>(0)
);
BOOST_CHECK_EQUAL(
quantile(complement( // Q very near 1 so P == epsilon < cdf(0) so should be exactly zero.
geometric_distribution<RealType>(static_cast<RealType>(0.25)),
static_cast<RealType>(1 - boost::math::tools::epsilon<RealType>()))),
static_cast<RealType>(0)
);
// Check that duff arguments throw domain_error:
BOOST_MATH_CHECK_THROW(
pdf( // Negative success_fraction!
geometric_distribution<RealType>(static_cast<RealType>(-0.25)),
static_cast<RealType>(0)), std::domain_error);
BOOST_MATH_CHECK_THROW(
pdf( // Success_fraction > 1!
geometric_distribution<RealType>(static_cast<RealType>(1.25)),
static_cast<RealType>(0)),
std::domain_error);
BOOST_MATH_CHECK_THROW(
pdf( // Negative k argument !
geometric_distribution<RealType>(static_cast<RealType>(0.25)),
static_cast<RealType>(-1)),
std::domain_error);
//BOOST_MATH_CHECK_THROW(
//pdf( // check limit on k (failures)
//geometric_distribution<RealType>(static_cast<RealType>(0.25)),
//std::numeric_limits<RealType>infinity()),
//std::domain_error);
BOOST_MATH_CHECK_THROW(
cdf( // Negative k argument !
geometric_distribution<RealType>(static_cast<RealType>(0.25)),
static_cast<RealType>(-1)),
std::domain_error);
BOOST_MATH_CHECK_THROW(
cdf( // Negative success_fraction!
geometric_distribution<RealType>(static_cast<RealType>(-0.25)),
static_cast<RealType>(0)), std::domain_error);
BOOST_MATH_CHECK_THROW(
cdf( // Success_fraction > 1!
geometric_distribution<RealType>(static_cast<RealType>(1.25)),
static_cast<RealType>(0)), std::domain_error);
BOOST_MATH_CHECK_THROW(
quantile( // Negative success_fraction!
geometric_distribution<RealType>(static_cast<RealType>(-0.25)),
static_cast<RealType>(0)), std::domain_error);
BOOST_MATH_CHECK_THROW(
quantile( // Success_fraction > 1!
geometric_distribution<RealType>(static_cast<RealType>(1.25)),
static_cast<RealType>(0)), std::domain_error);
check_out_of_range<geometric_distribution<RealType> >(0.5);
// End of check throwing 'duff' out-of-domain values.
{ // Compare geometric and negative binomial functions.
using boost::math::negative_binomial_distribution;
using boost::math::geometric_distribution;
RealType k = static_cast<RealType>(2.L);
RealType alpha = static_cast<RealType>(0.05L);
RealType p = static_cast<RealType>(0.5L);
BOOST_CHECK_CLOSE_FRACTION( // Successes parameter in negative binomial is 1 for geometric.
geometric_distribution<RealType>::find_lower_bound_on_p(k, alpha),
negative_binomial_distribution<RealType>::find_lower_bound_on_p(k, static_cast<RealType>(1), alpha),
tolerance);
BOOST_CHECK_CLOSE_FRACTION( // Successes parameter in negative binomial is 1 for geometric.
geometric_distribution<RealType>::find_upper_bound_on_p(k, alpha),
negative_binomial_distribution<RealType>::find_upper_bound_on_p(k, static_cast<RealType>(1), alpha),
tolerance);
BOOST_CHECK_CLOSE_FRACTION( // Should be identical - successes parameter is not used.
geometric_distribution<RealType>::find_maximum_number_of_trials(k, p, alpha),
negative_binomial_distribution<RealType>::find_maximum_number_of_trials(k, p, alpha),
tolerance);
}
//geometric::find_upper_bound_on_p(k, alpha);
return;
} // template <class RealType> void test_spots(RealType) // Any floating-point type RealType.
BOOST_AUTO_TEST_CASE( test_main )
{
// Check that can generate geometric distribution using the two convenience methods:
using namespace boost::math;
geometric g05d(0.5); // Using typedef - default type is double.
geometric_distribution<> g05dd(0.5); // Using default RealType double.
// Basic sanity-check spot values.
// Test some simple double only examples.
geometric_distribution<double> mydist(0.25);
// success fraction == 0.25 == 25% or 1 in 4 successes.
// Note: double values (matching the distribution definition) avoid the need for any casting.
// Check accessor functions return exact values for double at least.
BOOST_CHECK_EQUAL(mydist.success_fraction(), static_cast<double>(1./4.));
//cout << numeric_limits<RealType>::epsilon() << endl;
// (Parameter value, arbitrarily zero, only communicates the floating point type).
#ifdef TEST_FLOAT
test_spots(0.0F); // Test float.
#endif
#ifdef TEST_DOUBLE
test_spots(0.0); // Test double.
#endif
#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
#ifdef TEST_LDOUBLE
test_spots(0.0L); // Test long double.
#endif
#if !BOOST_WORKAROUND(BOOST_BORLANDC, BOOST_TESTED_AT(0x582)) && !defined(BOOST_MATH_NO_REAL_CONCEPT_TESTS)
#if defined(TEST_REAL_CONCEPT) && !defined(BOOST_MATH_NO_REAL_CONCEPT_TESTS)
test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
#endif
#endif
#else
std::cout << "<note>The long double tests have been disabled on this platform "
"either because the long double overloads of the usual math functions are "
"not available at all, or because they are too inaccurate for these tests "
"to pass.</note>" << std::endl;
#endif
} // BOOST_AUTO_TEST_CASE( test_main )
/*
*/
Reference in New Issue View Git Blame Copy Permalink