#5095 closed Bugs (fixed)
incorrect results from hypergeometric pdf
Reported by: | David Koes <dkoes@…> | Owned by: | johnmaddock |
---|---|---|---|
Milestone: | To Be Determined | Component: | math |
Version: | Boost 1.45.0 | Severity: | Problem |
Keywords: | Cc: |
Description
#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.
Attachments (5)
Change History (19)
Changed 6 years ago by David Koes <dkoes@…>
comment:1 Changed 6 years ago by 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());
comment:2 Changed 6 years ago by 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());
comment:3 Changed 6 years ago by johnmaddock
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.
comment:4 Changed 6 years ago by 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?
comment:5 Changed 6 years ago by 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.
comment:6 Changed 6 years ago by 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.
comment:7 follow-up: ↓ 9 Changed 6 years ago by 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.
Changed 6 years ago by anonymous
comment:8 Changed 6 years ago by 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.
comment:9 in reply to: ↑ 7 Changed 6 years ago by 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
comment:10 Changed 6 years ago by 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.
Changed 6 years ago by anonymous
comment:11 Changed 6 years ago by David Koes <dkoes@…>
Bingo. That fixed it. Thanks!
comment:12 Changed 6 years ago by johnmaddock
comment:13 Changed 6 years ago by anonymous
- Resolution set to fixed
- Status changed from new to closed
This was fixed in revision [68346].
example