bugGNU Octave - Bugs: bug #35827, [octave forge] (statistics)...

 
 

bug #35827: [octave forge] (statistics) hygepdf and hygecdf functions are numerically unstable for large populations

Submitter:  None
Submitted:  Tue 13 Mar 2012 01:19:38 PM UTC
   
 
Category:  Octave Package Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  Joachim Arts Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * dev
Operating System:  * Any Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Tue 26 Jul 2022 06:36:42 PM UTC, comment #9: 

no concerns with the pushed function for 8 months. closing report.

Nicholas Jankowski <nrjank>
Group Member
Thu 02 Dec 2021 01:34:38 AM UTC, comment #8: 
Nir Krakauer <nir_krakauer>
Tue 30 Nov 2021 10:39:57 PM UTC, comment #7: 

following the comment #2 and #3 recommendations, attached is a patch to hygepdf that replaces the pdf(k) lines with the direct calls to gammaln vs calling bincoeff.  as hygecdf just calls hygepdf, i think this takes care of both?

modification passes all tests after I added eps tolerances to the BISTs.  also noticed it's 3-5x faster.

I just used the comment #2 formula directly (and modified for the array case). Didn't check to see if there are any algebraic simplifications available there, but none were obvious.

I didn't see any other corner cases in the comments to test, so those with more intimate knowledge of the function may want to ensure this entirely fixes the bug.

(file #52393)

Nicholas Jankowski <nrjank>
Group Member
Wed 24 Nov 2021 05:55:21 PM UTC, comment #6: 

it doesnt seem the comment #2 fix has been implemented yet, as the issue is still present in Octave 6.4.0, Statistics 1.4.2

Nicholas Jankowski <nrjank>
Group Member
Sun 20 Nov 2016 02:51:26 PM UTC, comment #5: 

These functions should be part of the Octave-Forge statistics package which will happen at the next release.

Rik <rik5>
Group administrator
Sat 19 Nov 2016 08:33:37 PM UTC, comment #4: 

This issue is still present in Octave 4.2.0.

Hartmut <hardy>
Wed 14 Mar 2012 04:09:53 PM UTC, comment #3: 

Marco,

That's the expression I used as well, but perhaps I mistyped something.  At any rate, I think this is worth adopting as it is definitely more robust as a numerical algorithm.

Rik <rik5>
Group administrator
Wed 14 Mar 2012 09:37:59 AM UTC, comment #2: 

Rik,

the correct answer is 7.60282189600926e-03 with Octave too.


pdf(k) = exp(gammaln(m+1)-gammaln(m-x(k)+1)-gammaln(x(k)+1)+...
gammaln(t-m+1)-gammaln(t-m-n+x(k)+1)-gammaln(n-x(k)+1)-...
gammaln(t+1)+gammaln(t-n+1)+gammaln(n+1));


but is it safe to always use this expression? I think so, the difference should be O(eps).

Marco

Marco Caliari <caliari>
Group Member
Tue 13 Mar 2012 11:22:33 PM UTC, comment #1: 

Just to be clear, I get the following when running the sample code


hygecdf(10,10000,200,1000)
error: discrete_cdf: P must be a nonzero, non-negative vector


If I replace the cdf with pdf then Octave produces NaN.

What does Matlab return for


nchoosek (9800, 990)


I'm guessing it returns Inf which indicates that they are not using the straightforward implementation of two binomial coefficients divided by a third.

The algorithm for hygepdf is pretty straightforward with two binomial coefficients divided by a third (http://en.wikipedia.org/wiki/Hypergeometric_distribution).  Octave implements this directly which you can see by looking at the code in hygepdf.m.

The actual binomial coefficient in bincoeff.m is the exponential of some gammaln functions.  Thus, it would be pretty easy to pull the definition from bincoeff into hygepdf.  You could then make use of the rules for multiplying and dividing exponentials to add or subtract terms BEFORE applying the final exponential.

I just did a quick test with the example you provided and it does not produce NaN.  Instead the answer is 0.0043171.  This is a lot better but at this point I'm not sure which is more accurate: Matlab or Octave?  How did you determine that the correct answer was 0.0076?


Rik <rik5>
Group administrator
Tue 13 Mar 2012 01:19:38 PM UTC, original submission:  

The hypergeometric pdf and cdf return NaN for large populations.
For example hygecdf(10,10000,200,1000) returns NaN, when the correct result, also given in MATLAB is 0.007602821895.
Fixing by returning the binomial approximation for large populations is dangerous here because for the example above, the binomial approximation is binocdf(10,1000,200/10000)=0.010237. Thus this approximation can be quite poor already.

Anonymous

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #52393:  octave-statistics-hygepdf-bug35827.diff added by nrjank (3KiB - application/octet-stream - patch for hygepdf using comment #2 formula)
file #34233:  numHygeCdf.m added by ajf (203B - text/x-objcsrc - Wrapped the expression given in comment 2 into a ready-to-use numerical hygecdf function file)

 

Digest:
   bug dependencies.

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by nir_krakauer (Posted a comment)
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by hardy (Posted a comment)
  • -email is unavailable- added by ajf (Updated the item)
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by None (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 16 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-07-26 nrjank StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2021-12-22 nrjank StatusPatch Submitted Ready For Test
    2021-11-30 nrjank Attached File- Added octave-statistics-hygepdf-bug35827.diff, #52393
        StatusConfirmed Patch Submitted
    2019-06-04 mtmiller StatusPostponed Confirmed
        Summary[Octave Forge] The hygepdf and hygecdf functions are numerically unstable for large populations [octave forge] (statistics) hygepdf and hygecdf functions are numerically unstable for large populations
    2019-06-04 mtmiller Carbon-CopyRemoved 80942 -
    2019-06-04 rik5 Carbon-CopyRemoved 72865 -
    2019-06-04 rik5 CategoryLibraries Octave Package
        SummaryThe hygepdf and hygecdf functions are numerically unstable for large populations [Octave Forge] The hygepdf and hygecdf functions are numerically unstable for large populations
    2015-06-16 ajf Attached File- Added numHygeCdf.m, #34233
    2015-06-04 mtmiller StatusNone Postponed
        Release3.6.1 dev
        Operating SystemMicrosoft Windows Any
    2012-03-13 rik5 Dependencies- Depends on bugs #34362

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code