From fdbb3ed7ce85b8b199bb8428969d77cb6a46fa7b Mon Sep 17 00:00:00 2001 From: pabristow Date: Mon, 14 Jun 2021 17:46:18 +0100 Subject: [PATCH 1/2] Corrected formula for kurtosis to match Mathematica. New Test passes. --- .../math/distributions/hypergeometric.hpp | 37 +++++++++++++++---- test/__temporary_test.cpp | 22 ----------- test/test_hypergeometric_dist.cpp | 2 +- 3 files changed, 31 insertions(+), 30 deletions(-) delete mode 100644 test/__temporary_test.cpp diff --git a/include/boost/math/distributions/hypergeometric.hpp b/include/boost/math/distributions/hypergeometric.hpp index b0915c5c4..1ce3291ec 100644 --- a/include/boost/math/distributions/hypergeometric.hpp +++ b/include/boost/math/distributions/hypergeometric.hpp @@ -1,5 +1,6 @@ // Copyright 2008 Gautam Sewani // Copyright 2008 John Maddock +// Copyright 2021 Paul A. Bristow // // Use, modification and distribution are subject to the // Boost Software License, Version 1.0. @@ -231,6 +232,8 @@ namespace boost { namespace math { return static_cast(detail::hypergeometric_quantile(RealType(1 - c.param), c.param, c.dist.defective(), c.dist.sample_count(), c.dist.total(), Policy())); } // quantile + // https://www.wolframalpha.com/input/?i=kurtosis+hypergeometric+distribution + template inline RealType mean(const hypergeometric_distribution& dist) { @@ -269,13 +272,33 @@ namespace boost { namespace math { template inline RealType kurtosis_excess(const hypergeometric_distribution& dist) { - RealType r = static_cast(dist.defective()); - RealType n = static_cast(dist.sample_count()); - RealType N = static_cast(dist.total()); - RealType t1 = N * N * (N - 1) / (r * (N - 2) * (N - 3) * (N - r)); - RealType t2 = (N * (N + 1) - 6 * N * (N - r)) / (n * (N - n)) - + 3 * r * (N - r) * (N + 6) / (N * N) - 6; - return t1 * t2; + //RealType r = static_cast(dist.defective()); + //RealType n = static_cast(dist.sample_count()); + //RealType N = static_cast(dist.total()); + //RealType t1 = N * N * (N - 1) / (r * (N - 2) * (N - 3) * (N - r)); + //RealType t2 = (N * (N + 1) - 6 * N * (N - r)) / (n * (N - n)) + // + 3 * r * (N - r) * (N + 6) / (N * N) - 6; + //return t1 * t2; + + // https://www.wolframalpha.com/input/?i=kurtosis+hypergeometric+distribution shown as plain text: + // mean | (m n)/N + // standard deviation | sqrt((m n(N - m) (N - n))/(N - 1))/N + // variance | (m n(1 - m/N) (N - n))/((N - 1) N) + // skewness | (sqrt(N - 1) (N - 2 m) (N - 2 n))/((N - 2) sqrt(m n(N - m) (N - n))) + // kurtosis | ((N - 1) N^2 ((3 m(N - m) (n^2 (-N) + (n - 2) N^2 + 6 n(N - n)))/N^2 - 6 n(N - n) + N(N + 1)))/(m n(N - 3) (N - 2) (N - m) (N - n)) + + // Kurtosis[HypergeometricDistribution[n, m, N]] + RealType m = static_cast(dist.defective()); // Failures or success events. (Also symbols K or M are used). + RealType n = static_cast(dist.sample_count()); // draws or trials. + RealType n2 = n * n; // n^2 + RealType N = static_cast(dist.total()); // Total population from which n draws or trials are made. + RealType N2 = N * N; // N^2 + + // result = ((N - 1) N^2 ((3 m(N - m) (n^2 (-N) + (n - 2) N^2 + 6 n(N - n)))/N^2 - 6 n(N - n) + N(N + 1)))/(m n(N - 3) (N - 2) (N - m) (N - n)); + RealType result = ((N-1)*N2*((3*m*(N-m)*(n2*(-N)+(n-2)*N2+6*n*(N-n)))/N2-6*n*(N-n)+N*(N+1)))/(m*n*(N-3)*(N-2)*(N-m)*(N-n)); + // Agrees with kurtosis hypergeometric distribution(50,200,500) kurtosis = 2.96917 + // N[ kurtosis[hypergeometricdistribution(50,200,500)], 55] 2.969174035736058474901169623721804275002985337280263464 + return result; } // RealType kurtosis_excess(const hypergeometric_distribution& dist) template diff --git a/test/__temporary_test.cpp b/test/__temporary_test.cpp deleted file mode 100644 index f9e827328..000000000 --- a/test/__temporary_test.cpp +++ /dev/null @@ -1,22 +0,0 @@ -// Copyright John Maddock 2018 -// Distributed under 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) - -// This test file exists to output diagnostic info for tests failing in the online matrix -// for perplexing reasons, it's contents are subject to constant change!! -// -#define BOOST_MATH_INSTRUMENT - -#include -#include - -int main() -{ - std::cout << "EllintD(1, -1) = " << std::setprecision(20) << boost::math::ellint_d(1.0, -1.0) << std::endl; - std::cout << "EllintD(6.4851474761962890625e-01L , -7.6733188629150390625e+00L) = " << std::setprecision(20) << boost::math::ellint_d(6.4851474761962890625e-01, -7.6733188629150390625e+00) << std::endl; - - std::cout << "EllintD(1, -1) = " << std::setprecision(20) << boost::math::ellint_d(1.0L, -1.0L) << std::endl; - std::cout << "EllintD(6.4851474761962890625e-01L , -7.6733188629150390625e+00L) = " << std::setprecision(20) << boost::math::ellint_d(6.4851474761962890625e-01L , -7.6733188629150390625e+00L) << std::endl; -} - diff --git a/test/test_hypergeometric_dist.cpp b/test/test_hypergeometric_dist.cpp index 12e9871ee..0c0d01559 100644 --- a/test/test_hypergeometric_dist.cpp +++ b/test/test_hypergeometric_dist.cpp @@ -468,7 +468,7 @@ void test_spots(RealType /*T*/, const char* type_name) BOOST_CHECK_CLOSE(mode(d), static_cast(20), tolerance); BOOST_CHECK_CLOSE(variance(d), static_cast(10.821643286573146292585170340681L), tolerance); BOOST_CHECK_CLOSE(skewness(d), static_cast(0.048833071022952084732902910189366L), tolerance); - BOOST_CHECK_CLOSE(kurtosis_excess(d), static_cast(2.5155486690782804816404001878293L), tolerance); + BOOST_CHECK_CLOSE(kurtosis_excess(d), static_cast(2.969174035736058474901169623721804275002985337280263464L), tolerance); BOOST_CHECK_CLOSE(kurtosis(d), kurtosis_excess(d) + 3, tolerance); BOOST_CHECK_EQUAL(quantile(d, 0.5f), median(d)); From 4382d2a7c63a816887978ee943a0cb45358821d5 Mon Sep 17 00:00:00 2001 From: pabristow Date: Wed, 16 Jun 2021 11:42:37 +0100 Subject: [PATCH 2/2] More changes to hypergeometric distribution documentation. --- doc/distributions/hypergeometric.qbk | 26 +++++++++- doc/html/index.html | 2 +- doc/html/indexes/s01.html | 6 ++- doc/html/indexes/s02.html | 2 +- doc/html/indexes/s03.html | 2 +- doc/html/indexes/s04.html | 17 ++++++- doc/html/indexes/s05.html | 30 ++++++++++- doc/html/math_toolkit/barycentric.html | 12 ++--- doc/html/math_toolkit/conventions.html | 2 +- .../dist_ref/dists/hypergeometric_dist.html | 37 ++++++++++++-- doc/html/math_toolkit/navigation.html | 2 +- doc/sf/hypergeometric.qbk | 6 ++- .../math/distributions/hypergeometric.hpp | 23 +++------ test/test_hypergeometric_dist.cpp | 50 +++++++++++-------- 14 files changed, 158 insertions(+), 59 deletions(-) diff --git a/doc/distributions/hypergeometric.qbk b/doc/distributions/hypergeometric.qbk index f89418992..89b9ce46f 100644 --- a/doc/distributions/hypergeometric.qbk +++ b/doc/distributions/hypergeometric.qbk @@ -15,7 +15,7 @@ typedef RealType value_type; typedef Policy policy_type; // Construct: - hypergeometric_distribution(unsigned r, unsigned n, unsigned N); + hypergeometric_distribution(unsigned r, unsigned n, unsigned N); // r=defective/failures/success, n=trials/draws, N=total population. // Accessors: unsigned total()const; unsigned defective()const; @@ -73,6 +73,15 @@ Returns the number of objects /r/ in population /N/ which are defective. Returns the number of objects /n/ which are sampled from the population /N/. +[warning Both naming/symbol and order of parameters is confusing with no two implementations the same! +See Wolfram Mathematica +[@https://mathworld.wolfram.com/HypergeometricDistribution.html Hypergeometric Distribution] +and Wikipedia +[@https://en.wikipedia.org/wiki/Hypergeometric_distribution Hypergeometric Distribution] +and Python +[@https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.hypergeom.html#scipy.stats.hypergeom scipy.stats.hypergeom]. +] + [h4 Non-member Accessors] All the [link math_toolkit.dist_ref.nmp usual non-member accessor functions] @@ -137,6 +146,9 @@ our implementation against some spot values computed using the online calculator here [@http://stattrek.com/Tables/Hypergeometric.aspx http://stattrek.com/Tables/Hypergeometric.aspx]. Finally we test accuracy against some high precision test data using this implementation and NTL::RR. +Spot test values for moments (mean to kurtosis) are from Mathematica [@https://mathworld.wolfram.com/HypergeometricDistribution.html Hypergeometric Distribution] +and agree with an implementation of Wikipedia [@https://en.wikipedia.org/wiki/Hypergeometric_distribution Hypergeometric Distribution] +and Python [@https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.hypergeom.html#scipy.stats.hypergeom scipy.stats.hypergeom]. [h4 Implementation] @@ -198,7 +210,7 @@ lgamma. However, in this area where N > 104729, the user should expect to lose around log[sub 10]N decimal digits during the calculation in the worst case. -The CDF and its complement is calculated by directly summing the PDF's. +The CDF and its complement is calculated by directly summing the PDFs. We start by deciding whether the CDF, or its complement, is likely to be the smaller of the two and then calculate the PDF at /k/ (or /k+1/ if we're calculating the complement) and calculate successive PDF values via the @@ -219,6 +231,16 @@ calculated via: [equation hypergeometric6] +[note The kurtosis formula above is not quite correct and kurtosis excess is now calculated +from +[@https://www.wolframalpha.com/input/?i=kurtosis+hypergeometric+distribution Wolfram Alpha hypergeometric distribution kurtosis]. +(The hypergeometric distribution kurtosis excess is mentioned in +[@https://mathworld.wolfram.com/HypergeometricDistribution.html Wolfram Hypergeometric distribution] +but coyly only described as ['\"the kurtosis excess is given by a complicated expression.\"]). +This has been found numerically equivalent to the +[@https://en.wikipedia.org/wiki/Hypergeometric_distribution Wikipedia hypergeometric kurtosis excess formula]. +] + [endsect] [/ hypergeometric.qbk diff --git a/doc/html/index.html b/doc/html/index.html index f930626cd..5f71a8555 100644 --- a/doc/html/index.html +++ b/doc/html/index.html @@ -134,7 +134,7 @@ This manual is also available in -

Last revised: March 30, 2021 at 17:44:57 GMT

+

Last revised: June 16, 2021 at 10:35:24 GMT

diff --git a/doc/html/indexes/s01.html b/doc/html/indexes/s01.html index 674fd20ff..ed5fcc865 100644 --- a/doc/html/indexes/s01.html +++ b/doc/html/indexes/s01.html @@ -24,7 +24,7 @@

-Function Index

+Function Index

1 2 4 A B C D E F G H I J K L M N O P Q R S T U V W X Y Z

@@ -2502,6 +2502,10 @@
  • +

    means

    +
    +
    • +

  • means_and_covariance

    • diff --git a/doc/html/indexes/s02.html b/doc/html/indexes/s02.html index de7348e8d..117811624 100644 --- a/doc/html/indexes/s02.html +++ b/doc/html/indexes/s02.html @@ -24,7 +24,7 @@

    -Class Index

    +Class Index

    A B C D E F G H I K L M N O P Q R S T U V W

    diff --git a/doc/html/indexes/s03.html b/doc/html/indexes/s03.html index 6c570fcb3..b5a850a11 100644 --- a/doc/html/indexes/s03.html +++ b/doc/html/indexes/s03.html @@ -24,7 +24,7 @@

    -Typedef Index

    +Typedef Index

    A B C D E F G H I K L N P R S T U V W

    diff --git a/doc/html/indexes/s04.html b/doc/html/indexes/s04.html index 3d40d41e5..e93374726 100644 --- a/doc/html/indexes/s04.html +++ b/doc/html/indexes/s04.html @@ -24,7 +24,7 @@

    -Macro Index

    +Macro Index

    B F

    @@ -305,6 +305,13 @@
  • +

    BOOST_MATH_STANDALONE

    +
    +
    • +

  • BOOST_MATH_STD_USING

    F diff --git a/doc/html/indexes/s05.html b/doc/html/indexes/s05.html index e06c70a53..6e8b4bbce 100644 --- a/doc/html/indexes/s05.html +++ b/doc/html/indexes/s05.html @@ -23,7 +23,7 @@

    • -Index

    +Index

    1 2 4 5 7 A B C D E F G H I J K L M N O P Q R S T U V W X Y Z

    @@ -1218,6 +1218,13 @@
  • +

    BOOST_MATH_STANDALONE

    +
    +
    • +

  • BOOST_MATH_STD_USING

    • @@ -3485,6 +3501,7 @@

  • Error Handling

    @@ -3494,6 +3511,7 @@
    • @@ -4214,6 +4232,10 @@
  • +

    F Distribution Examples

    +
    +
    • +

  • Facets for Floating-Point Infinities and NaNs

    • nonfinite_num_get

      • @@ -4934,6 +4956,7 @@
    • +

      means

      +
      +
      • +

    • means_and_covariance

      • @@ -9083,6 +9110,7 @@

    • BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS

    • BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS

    • BOOST_MATH_NO_REAL_CONCEPT_TESTS

      • +

    • BOOST_MATH_STANDALONE

    • bug

    • constants

    • cpp_dec_float

      • diff --git a/doc/html/math_toolkit/barycentric.html b/doc/html/math_toolkit/barycentric.html index 042e1b6af..8ac6e284f 100644 --- a/doc/html/math_toolkit/barycentric.html +++ b/doc/html/math_toolkit/barycentric.html @@ -32,7 +32,7 @@ 
      #include <boost/math/interpolators/barycentric_rational.hpp>
 
-namespace boost{ namespace math{
+namespace boost{ namespace math{ namespace interpolators{
     template<class Real>
     class barycentric_rational
     {
@@ -53,7 +53,7 @@
         std::vector<Real>&& return_y();
     };
 
-}}
+}}}

      @@ -78,7 +78,7 @@ 
      std::vector<double> x(500);
 std::vector<double> y(500);
 // populate x, y, then:
-boost::math::barycentric_rational<double> interpolant(std::move(x), std::move(y));
+boost::math::interpolators::barycentric_rational<double> interpolant(std::move(x), std::move(y));

      This implicitly calls the constructor with approximation order 3, and hence @@ -87,7 +87,7 @@ A call to the constructor with an explicit approximation order is demonstrated below

      -
      boost::math::barycentric_rational<double> interpolant(std::move(x), std::move(y), 5);
+
      boost::math::interpolators::barycentric_rational<double> interpolant(std::move(x), std::move(y), 5);
 
      
 
      
       To evaluate the interpolant, simply use
@@ -168,7 +168,7 @@
     // We'll skip the code for filling the above vectors with data for now...
 
     // Now we want to interpolate this potential at any r:
-    boost::math::barycentric_rational<double> b(r.data(), mrV.data(), r.size());
+    boost::math::interpolators::barycentric_rational<double> b(r.data(), mrV.data(), r.size());
 
     for (size_t i = 1; i < 8; ++i)
     {
@@ -242,7 +242,7 @@
     // start by creating 2 ranges for the x and y values:
     auto x_range = boost::adaptors::keys(r);
     auto y_range = boost::adaptors::values(r);
-    boost::math::barycentric_rational<double> b(x_range.begin(), x_range.end(), y_range.begin());
+    boost::math::interpolators::barycentric_rational<double> b(x_range.begin(), x_range.end(), y_range.begin());
     //
     // We'll use a lambda expression to provide the functor to our root finder, since we want
     // the abscissa value that yields 3, not zero.  We pass the functor b by value to the
diff --git a/doc/html/math_toolkit/conventions.html b/doc/html/math_toolkit/conventions.html
index 17177b15c..5b815d0c0 100644
--- a/doc/html/math_toolkit/conventions.html
+++ b/doc/html/math_toolkit/conventions.html
@@ -27,7 +27,7 @@
 Document Conventions

    - +

    This documentation aims to use of the following naming and formatting conventions. diff --git a/doc/html/math_toolkit/dist_ref/dists/hypergeometric_dist.html b/doc/html/math_toolkit/dist_ref/dists/hypergeometric_dist.html index 1426a6ea6..7b11aa984 100644 --- a/doc/html/math_toolkit/dist_ref/dists/hypergeometric_dist.html +++ b/doc/html/math_toolkit/dist_ref/dists/hypergeometric_dist.html @@ -41,7 +41,7 @@ typedef RealType value_type; typedef Policy policy_type; // Construct: - hypergeometric_distribution(unsigned r, unsigned n, unsigned N); + hypergeometric_distribution(unsigned r, unsigned n, unsigned N); // r=defective/failures/success, n=trials/draws, N=total population. // Accessors: unsigned total()const; unsigned defective()const; @@ -128,6 +128,18 @@ Returns the number of objects n which are sampled from the population N.

    +
    + + + + + +
    [Warning]Warning

    + Both naming/symbol and order of parameters is confusing with no two implementations + the same! See Wolfram Mathematica Hypergeometric + Distribution and Wikipedia Hypergeometric + Distribution and Python scipy.stats.hypergeom. +

    Non-member @@ -218,7 +230,10 @@ We also sanity check our implementation against some spot values computed using the online calculator here http://stattrek.com/Tables/Hypergeometric.aspx. Finally we test accuracy against some high precision test data using this - implementation and NTL::RR. + implementation and NTL::RR. Spot test values for moments (mean to kurtosis) + are from Mathematica Hypergeometric + Distribution and agree with an implementation of Wikipedia Hypergeometric + Distribution and Python scipy.stats.hypergeom.

    @@ -299,7 +314,7 @@ log10N decimal digits during the calculation in the worst case.

    - The CDF and its complement is calculated by directly summing the PDF's. + The CDF and its complement is calculated by directly summing the PDFs. We start by deciding whether the CDF, or its complement, is likely to be the smaller of the two and then calculate the PDF at k (or k+1 if we're calculating the complement) and calculate @@ -327,6 +342,22 @@

    +
    + + + + + +
    [Note]Note

    + The kurtosis formula above is not quite correct and kurtosis excess is + now calculated from Wolfram + Alpha hypergeometric distribution kurtosis. (The hypergeometric + distribution kurtosis excess is mentioned in Wolfram + Hypergeometric distribution but coyly only described as "the + kurtosis excess is given by a complicated expression."). + This has been found numerically equivalent to the Wikipedia + hypergeometric kurtosis excess formula. +

    diff --git a/doc/html/math_toolkit/navigation.html b/doc/html/math_toolkit/navigation.html index dc305c118..7748e49df 100644 --- a/doc/html/math_toolkit/navigation.html +++ b/doc/html/math_toolkit/navigation.html @@ -27,7 +27,7 @@ Navigation

    - +

    Boost.Math documentation is provided in both HTML and PDF formats. diff --git a/doc/sf/hypergeometric.qbk b/doc/sf/hypergeometric.qbk index 7882fd55b..9978cf14f 100644 --- a/doc/sf/hypergeometric.qbk +++ b/doc/sf/hypergeometric.qbk @@ -275,7 +275,9 @@ readers will have to refer to the code for the transitions between methods, as w In almost all cases where /z < 0/ we use [@https://en.wikipedia.org/wiki/Confluent_hypergeometric_function Kummer's relation] to make /z/ positive: -[:[/\large $$_1F_1(a, b; -z) = e^{-z} {_1F_1}(a, b; z)$$] +[:[/\large $$_1F_1(a, b; -z) = e^{-z} {_1F_1}(b-a, b; z)$$] +[/ https://github.com/boostorg/math/issues/638] +[/See https://dlmf.nist.gov/13.2#E39] [$../equations/hypergeometric_1f1_12.svg]] The main series representation @@ -532,4 +534,4 @@ error inherent in calculating the N'th term via logarithms. Distributed under 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). -] \ No newline at end of file +] diff --git a/include/boost/math/distributions/hypergeometric.hpp b/include/boost/math/distributions/hypergeometric.hpp index 1ce3291ec..8ec62ccc4 100644 --- a/include/boost/math/distributions/hypergeometric.hpp +++ b/include/boost/math/distributions/hypergeometric.hpp @@ -17,7 +17,6 @@ #include #include - namespace boost { namespace math { template > @@ -27,7 +26,7 @@ namespace boost { namespace math { typedef RealType value_type; typedef Policy policy_type; - hypergeometric_distribution(unsigned r, unsigned n, unsigned N) // Constructor. + hypergeometric_distribution(unsigned r, unsigned n, unsigned N) // Constructor. r=defective/failures/success, n=trials/draws, N=total population. : m_n(n), m_N(N), m_r(r) { static const char* function = "boost::math::hypergeometric_distribution<%1%>::hypergeometric_distribution"; @@ -40,7 +39,7 @@ namespace boost { namespace math { return m_N; } - unsigned defective()const + unsigned defective()const // successes/failures/events { return m_r; } @@ -85,9 +84,9 @@ namespace boost { namespace math { private: // Data members: - unsigned m_n; // number of items picked - unsigned m_N; // number of "total" items - unsigned m_r; // number of "defective" items + unsigned m_n; // number of items picked or drawn. + unsigned m_N; // number of "total" items. + unsigned m_r; // number of "defective/successes/failures/events items. }; // class hypergeometric_distribution @@ -272,32 +271,22 @@ namespace boost { namespace math { template inline RealType kurtosis_excess(const hypergeometric_distribution& dist) { - //RealType r = static_cast(dist.defective()); - //RealType n = static_cast(dist.sample_count()); - //RealType N = static_cast(dist.total()); - //RealType t1 = N * N * (N - 1) / (r * (N - 2) * (N - 3) * (N - r)); - //RealType t2 = (N * (N + 1) - 6 * N * (N - r)) / (n * (N - n)) - // + 3 * r * (N - r) * (N + 6) / (N * N) - 6; - //return t1 * t2; - // https://www.wolframalpha.com/input/?i=kurtosis+hypergeometric+distribution shown as plain text: // mean | (m n)/N // standard deviation | sqrt((m n(N - m) (N - n))/(N - 1))/N // variance | (m n(1 - m/N) (N - n))/((N - 1) N) // skewness | (sqrt(N - 1) (N - 2 m) (N - 2 n))/((N - 2) sqrt(m n(N - m) (N - n))) // kurtosis | ((N - 1) N^2 ((3 m(N - m) (n^2 (-N) + (n - 2) N^2 + 6 n(N - n)))/N^2 - 6 n(N - n) + N(N + 1)))/(m n(N - 3) (N - 2) (N - m) (N - n)) - // Kurtosis[HypergeometricDistribution[n, m, N]] RealType m = static_cast(dist.defective()); // Failures or success events. (Also symbols K or M are used). RealType n = static_cast(dist.sample_count()); // draws or trials. RealType n2 = n * n; // n^2 RealType N = static_cast(dist.total()); // Total population from which n draws or trials are made. RealType N2 = N * N; // N^2 - // result = ((N - 1) N^2 ((3 m(N - m) (n^2 (-N) + (n - 2) N^2 + 6 n(N - n)))/N^2 - 6 n(N - n) + N(N + 1)))/(m n(N - 3) (N - 2) (N - m) (N - n)); RealType result = ((N-1)*N2*((3*m*(N-m)*(n2*(-N)+(n-2)*N2+6*n*(N-n)))/N2-6*n*(N-n)+N*(N+1)))/(m*n*(N-3)*(N-2)*(N-m)*(N-n)); // Agrees with kurtosis hypergeometric distribution(50,200,500) kurtosis = 2.96917 - // N[ kurtosis[hypergeometricdistribution(50,200,500)], 55] 2.969174035736058474901169623721804275002985337280263464 + // N[kurtosis[hypergeometricdistribution(50,200,500)], 55] 2.969174035736058474901169623721804275002985337280263464 return result; } // RealType kurtosis_excess(const hypergeometric_distribution& dist) diff --git a/test/test_hypergeometric_dist.cpp b/test/test_hypergeometric_dist.cpp index 0c0d01559..f1a4b8b40 100644 --- a/test/test_hypergeometric_dist.cpp +++ b/test/test_hypergeometric_dist.cpp @@ -1,6 +1,6 @@ // Copyright John Maddock 2008 -// Copyright Paul A. Bristow -// Copyright Gautam Sewani +// Copyright Paul A. Bristow 2021 +// Copyright Gautam Sewani 2008 // Use, modification and distribution are subject to the // Boost Software License, Version 1.0. @@ -165,8 +165,8 @@ void do_test_hypergeometric(const T& data, const char* type_name, const char* te // test hypergeometric against data: // result = boost::math::tools::test_hetero( - data, - bind_func(funcp, 0, 1, 2, 3), + data, + bind_func(funcp, 0, 1, 2, 3), extract_result(4)); handle_test_result(result, data[result.worst()], result.worst(), type_name, "hypergeometric PDF", test_name); @@ -180,8 +180,8 @@ void do_test_hypergeometric(const T& data, const char* type_name, const char* te // test hypergeometric against data: // result = boost::math::tools::test_hetero( - data, - bind_func(funcp, 0, 1, 2, 3), + data, + bind_func(funcp, 0, 1, 2, 3), extract_result(5)); handle_test_result(result, data[result.worst()], result.worst(), type_name, "hypergeometric CDF", test_name); @@ -195,8 +195,8 @@ void do_test_hypergeometric(const T& data, const char* type_name, const char* te // test hypergeometric against data: // result = boost::math::tools::test_hetero( - data, - bind_func(funcp, 0, 1, 2, 3), + data, + bind_func(funcp, 0, 1, 2, 3), extract_result(6)); handle_test_result(result, data[result.worst()], result.worst(), type_name, "hypergeometric CDF complement", test_name); std::cout << std::endl; @@ -234,7 +234,7 @@ void do_test_hypergeometric_quantile(const T& data, const char* type_name, const (void)ccp; #if !defined(TEST_QUANT) || (TEST_QUANT == 1) - boost::math::hypergeometric_distribution > > du(r, n, N); if((cp < 0.9) && (cp > boost::math::tools::min_value())) @@ -247,7 +247,7 @@ void do_test_hypergeometric_quantile(const T& data, const char* type_name, const } #endif #if !defined(TEST_QUANT) || (TEST_QUANT == 2) - boost::math::hypergeometric_distribution > > dl(r, n, N); if((cp < 0.9) && (cp > boost::math::tools::min_value())) { @@ -259,7 +259,7 @@ void do_test_hypergeometric_quantile(const T& data, const char* type_name, const } #endif #if !defined(TEST_QUANT) || (TEST_QUANT == 3) - boost::math::hypergeometric_distribution > > dn(r, n, N); if((cp < 0.9) && (cp > boost::math::tools::min_value())) @@ -272,7 +272,7 @@ void do_test_hypergeometric_quantile(const T& data, const char* type_name, const } #endif #if !defined(TEST_QUANT) || (TEST_QUANT == 4) - boost::math::hypergeometric_distribution > > dou(r, n, N); if((cp < 0.9) && (cp > boost::math::tools::min_value())) @@ -299,7 +299,7 @@ void do_test_hypergeometric_quantile(const T& data, const char* type_name, const } #endif #if !defined(TEST_QUANT) || (TEST_QUANT == 5) - boost::math::hypergeometric_distribution > > di(r, n, N); if((cp < 0.9) && (cp > boost::math::tools::min_value())) @@ -331,7 +331,7 @@ void do_test_hypergeometric_quantile(const T& data, const char* type_name, const template -void test_spot(unsigned x, unsigned n, unsigned r, unsigned N, +void test_spot(unsigned x, unsigned n, unsigned r, unsigned N, RealType p, RealType cp, RealType ccp, RealType tol) { // @@ -366,18 +366,18 @@ void test_spot(unsigned x, unsigned n, unsigned r, unsigned N, { using namespace boost::math::policies; - boost::math::hypergeometric_distribution > > du(r, n, N); BOOST_CHECK_EX(quantile(du, cp) >= x); BOOST_CHECK_EX(quantile(complement(du, ccp)) >= x); - boost::math::hypergeometric_distribution > > dl(r, n, N); BOOST_CHECK_EX(quantile(dl, cp) <= x); BOOST_CHECK_EX(quantile(complement(dl, ccp)) <= x); - boost::math::hypergeometric_distribution > > dn(r, n, N); BOOST_CHECK_EX(quantile(dn, cp) == x); @@ -443,7 +443,7 @@ void test_spots(RealType /*T*/, const char* type_name) static_cast(2e-16L), // limit of test data boost::math::tools::epsilon()); cout<<"Absolute tolerance:"<(20), tolerance); BOOST_CHECK_CLOSE(mode(d), static_cast(20), tolerance); - BOOST_CHECK_CLOSE(variance(d), static_cast(10.821643286573146292585170340681L), tolerance); - BOOST_CHECK_CLOSE(skewness(d), static_cast(0.048833071022952084732902910189366L), tolerance); + // Test values and code revised to correct kurtosis using Mathematica algorithm and test values. + // Also checked using boost::multiprecision::cpp_bin_float_quad at 128-bit precision. + // And compared to an implementation of the Wikipedia equation at https://en.wikipedia.org/wiki/Hypergeometric_distribution. + // https://github.com/boostorg/math/issues/639 hypergeometric kurtosis differs from Wolfram and scipy and Wikipedia. + // N[variance[hypergeometricdistribution(200,50,500)], 55] + BOOST_CHECK_CLOSE(variance(d), static_cast(10.82164328657314629258517034068136272545090180360721443L), tolerance); + // N[skewness[hypergeometricdistribution(200,50,500)], 55] + BOOST_CHECK_CLOSE(skewness(d), static_cast(0.04883307102295208473290291018936563084537649911337978351L), tolerance); + // https://www.wolframalpha.com/input/?i=N%5Bkurtosis%5Bhypergeometricdistribution%28200%2C50%2C500%29%5D%2C+55%5D+ + // N[kurtosis[hypergeometricdistribution(200,50,500)], 55] = 2.969174035736058474901169623721804275002985337280263464 BOOST_CHECK_CLOSE(kurtosis_excess(d), static_cast(2.969174035736058474901169623721804275002985337280263464L), tolerance); BOOST_CHECK_CLOSE(kurtosis(d), kurtosis_excess(d) + 3, tolerance); BOOST_CHECK_EQUAL(quantile(d, 0.5f), median(d)); @@ -496,6 +504,6 @@ BOOST_AUTO_TEST_CASE( test_main ) "to pass." << std::endl; #endif - + } // BOOST_AUTO_TEST_CASE( test_main )