incorrect results from hypergeometric pdf

[David Koes <dkoes@…>]

#include <boost/math/distributions/hypergeometric.hpp> #include <boost/math/policies/policy.hpp> #include <iostream>

using namespace std; using namespace boost;

int main() {

unsigned N = 16086184; unsigned n = 256004; unsigned Q = 251138; math::hypergeometric_distribution<double> hyper(n, Q, N); cout << math::pdf<double>(hyper, 4000) << " " << math::pdf<double>(hyper, 4001) << " " << math::pdf<double>(hyper, 4002) << "\n";

return 0;

}

Output: 0.00640003 1.11519e-09 0.00638443 The value for 4001 is incorrect (according to ​http://stattrek.com/Tables/Hypergeometric.aspx). In fact, every value where k is odd appears to be incorrect.

[David Koes <dkoes@…>]

example

[David Koes <dkoes@…>]

As a workaround, if the lgamma backup version of hypergeometric_pdf_lanczos_imp is used instead, correct results are obtained.

In /usr/local/include/boost/math/distributions/detail/hypergeometric_pdf.hpp: 414c414 < result = detail::hypergeometric_pdf_lanczos_imp(value_type(), x, r, n, N, evaluation_type(), forwarding_policy()); ---

result = detail::hypergeometric_pdf_lanczos_imp(value_type(), x, r, n, N, lanczos::undefined_lanczos(), forwarding_policy());

[David Koes <dkoes@…>]

Properly formatted diff of boost/math/distributions/detail/hypergeometric_pdf.hpp

414c414
<       result = detail::hypergeometric_pdf_lanczos_imp(value_type(), x, r, n, N, evaluation_type(), forwarding_policy());
---
>       result = detail::hypergeometric_pdf_lanczos_imp(value_type(), x, r, n, N, lanczos::undefined_lanczos(), forwarding_policy());

[John Maddock]

I've been unable to reproduce here either with VC10 or with gcc-4.4 on Ubuntu Linux, what compiler/platform are you on?

The output I see from your program is: 0.00640003 0.00639305 0.00638443 which I believe agrees with the online calculator.

[anonymous]

gcc version 4.4.3 (Ubuntu 4.4.3-4ubuntu5) on Ubuntu 10.4

Also, on Ubuntu 10.10 with boost version 1.40 and gcc version 4.4.5 (Ubuntu/Linaro 4.4.4-14ubuntu5) and OS X (10.6) with boost version 1.41 and gcc version 4.2.1 (Apple Inc. build 5664).

Perhaps for some reason your installation is already falling back on the lgamma version of hypergeometric_pdf_lanczos_imp?

[anonymous]

No it's calling for example:

hypergeometric_pdf_lanczos_imp<long double,boost::math::lanczos::lanczos13m53,boost::math::policies::policy<boost::math::policies::promote_float<0>,boost::math::policies::promote_double<0>,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy,boost::math::policies::default_policy> >(double formal, unsigned int x, unsigned int r, unsigned int n, unsigned int N, double formal, double formal);

I notice that you're using an out of date Boost distribution, can you:

• Please try the last release and if that still fails, then:
• Let me have the program output when BOOST_MATH_INSTRUMENT is #defined when building?

Thanks, John.

[David Koes <dkoes@…>]

program output with BOOST_MATH_INSTRUMENT defined

[David Koes <dkoes@…>]

The Ubuntu 10.04 box has boost 1.45 on it. I also tried it on two other machines with the version of boost that they already have. I have also grabbed and built the svn trunk on the OS X machine, and tried removing /usr/local/include/boost and reinstalling on the 10.4 machine. I still get the failure.

I have attached the output with BOOST_MATH_INSTRUMENT defined.

Thanks.

[anonymous]

Thanks for trying that, I really appreciate the help you're putting in to get to the bottom of this, unfortunately I still can't see what's gone wrong (other than something has), I'm also confused because there's really nothing in the code that should rely on the evenness of the random variable!

I'm attaching an updated hypergeometric_pdf.hpp with a lot more instrumentation, can I get you to try again? Just the failing case this time - otherwise the output is going to be huge! :(

Many thanks,

Still confused yours, John.

[anonymous]

Forgot to say - if the program output is too large to attach here - zip it up and mail me direct at john at johnmaddock.co.uk

Thanks! John.

[David Koes <dkoes@…>]

even more instrumentation (bad value only, gzipped)

[David Koes <dkoes@…>]

Attached. Here are some interesting results that may partially explain the reproducibility problem:

dkoes@quasar:~/tmp$ g++  hyper.cpp; ./a.out 
1.11519e-09
dkoes@quasar:~/tmp$ g++ -m32 hyper.cpp; ./a.out 
0.00639305
dkoes@quasar:~/tmp$ g++ -mfpmath=387 hyper.cpp; ./a.out 
0.00639305
dkoes@quasar:~/tmp$ g++ -mfpmath=sse hyper.cpp; ./a.out 
1.11519e-09
dkoes@quasar:~/tmp$ g++ -mno-sse hyper.cpp; ./a.out 
In file included from /usr/local/include/boost/config/no_tr1/cmath.hpp:21,
                 from /usr/local/include/boost/math/policies/error_handling.hpp:15,
                 from /usr/local/include/boost/math/distributions/detail/common_error_handling.hpp:12,
                 from /usr/local/include/boost/math/distributions/hypergeometric.hpp:12,
                 from hyper.cpp:1:
/usr/include/c++/4.4/cmath: In function ‘double std::abs(double)’:
/usr/include/c++/4.4/cmath:94: error: SSE register return with SSE disabled
1.11519e-09
dkoes@quasar:~/tmp$ g++ -mno-sse2 hyper.cpp; ./a.out 
1.11519e-09
dkoes@quasar:~/tmp$ g++ -mno-sse3 hyper.cpp; ./a.out 
1.11519e-09
dkoes@quasar:~/tmp$ g++ -mno-sse4 hyper.cpp; ./a.out 
1.11519e-09
dkoes@quasar:~/tmp$ g++ -DBOOST_MATH_SPECIAL_FUNCTIONS_LANCZOS_SSE2 hyper.cpp; ./a.out 
1.11519e-09

[anonymous]

Thanks for that, some of the values in the "exponents" table are getting truncated, I'm attaching an updated header with some extra typecasts present to try and prevent this, can you give this a try?

Thanks, John.

[David Koes <dkoes@…>]

Bingo. That fixed it. Thanks!

[John Maddock]

(In [68347]) Added note about fixed bug. Refs #5095.

[anonymous]

This was fixed in revision [68346].

[John Maddock]

(In [68367]) Merge fix for issue 5095. Refs #5095.