Modify

Ticket #5095 (closed Bugs: fixed)

Opened 3 years ago

Last modified 3 years ago

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

hyper.cpp Download (458 bytes) - added by David Koes <dkoes@…> 3 years ago.
example
BMIdump Download (16.1 KB) - added by David Koes <dkoes@…> 3 years ago.
program output with BOOST_MATH_INSTRUMENT defined
hypergeometric_pdf.hpp Download (15.5 KB) - added by anonymous 3 years ago.
BMIdumpmore.gz Download (8.8 KB) - added by David Koes <dkoes@…> 3 years ago.
even more instrumentation (bad value only, gzipped)
hypergeometric_pdf.2.hpp Download (15.5 KB) - added by anonymous 3 years ago.

Change History

Changed 3 years ago by David Koes <dkoes@…>

example

comment:1 Changed 3 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 3 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 3 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 3 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 3 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.

Changed 3 years ago by David Koes <dkoes@…>

program output with BOOST_MATH_INSTRUMENT defined

comment:6 Changed 3 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 3 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 3 years ago by anonymous

comment:8 Changed 3 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.

Changed 3 years ago by David Koes <dkoes@…>

even more instrumentation (bad value only, gzipped)

comment:9 in reply to: ↑ 7 Changed 3 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 3 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 3 years ago by anonymous

comment:11 Changed 3 years ago by David Koes <dkoes@…>

Bingo. That fixed it. Thanks!

comment:12 Changed 3 years ago by johnmaddock

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

comment:13 Changed 3 years ago by anonymous

  • Status changed from new to closed
  • Resolution set to fixed

This was fixed in revision [68346].

comment:14 Changed 3 years ago by johnmaddock

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

View

Add a comment

Modify Ticket

Change Properties
<Author field>
Action
as closed
The resolution will be deleted. Next status will be 'reopened'
Author


E-mail address and user name can be saved in the Preferences.

 
Note: See TracTickets for help on using tickets.