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

 
 

bug #39636: [octave forge] (statistics) nbinpdf has floating point error

Submitter:  None
Submitted:  Wed 31 Jul 2013 07:10:55 PM UTC
   
 
Category:  Octave Package Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Wont Fix Assigned to:  None
Originator Name:  ken massey 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

Thu 28 Jul 2022 02:16:15 PM UTC, comment #11: 

Andreas (statistics package maintainer) decided to close this item.

Please use the new tracker if there is anything worth to be followed up: https://github.com/gnu-octave/statistics/issues

Kai Torben Ohlhus <siko1056>
Group Member
Sun 28 Feb 2016 03:44:59 PM UTC, comment #10: 

Andrew, I do not remember the details, but my proposal fixes the original bug and all successive variations, isn't it?

Marco Caliari <caliari>
Group Member
Sat 27 Feb 2016 10:46:54 AM UTC, comment #9: 

Marco, how would the change you propose help?  There would still be a discontinuity in bincoeff.

Lachlan Andrew <lachlan>
Thu 01 Aug 2013 02:02:57 PM UTC, comment #8: 

Hi,

why not to change the relevant part in nbinpdf.m to


if (isscalar (n) && isscalar (p))
  pdf(k) = bincoeff (x(k) + n - 1, x(k)) .* (p ^ n) .* ((1 - p) .^ x(k));
else
  pdf(k) = bincoeff (x(k) + n(k) - 1, x(k)) .* (p(k) .^ n(k)) .* ((1 - p(k)) .^ x(k));
endif


according to the definition given in comment #2?
About bincoeff.m, I do not understand the formula given on line 92. Where does it come from?

Marco

Marco Caliari <caliari>
Group Member
Thu 01 Aug 2013 10:54:13 AM UTC, comment #7: 

Thanks for your help with this issue.  I believe the R statistical software might be helpful.

http://stat.ethz.ch/R-manual/R-patched/library/stats/html/NegBinomial.html

Following the source code, I see:

http://svn.r-project.org/R/trunk/src/nmath/pnbinom.c

http://svn.r-project.org/R/trunk/src/nmath/pbeta.c

http://svn.r-project.org/R/trunk/src/nmath/toms708.c
(bratio function for incomplete beta)

Maybe there is something in there that could help you patch the discontinuities.

Anonymous
Thu 01 Aug 2013 03:18:17 AM UTC, comment #6: 

Ok, I see the point you are making. There are similar discontinuities at other integer values of N +/- a few epsilon. Visually,


octave:44> n = 3;
octave:45> d = -1e3:1e3;
octave:46> plot (d, nbinpdf (0, n + d*eps, .5));


Mike Miller <mtmiller>
Group Member
Thu 01 Aug 2013 01:48:15 AM UTC, comment #5: 

We obviously can't just compute the binomial coefficient using the gamma function for all arguments, since it will fail for everything but the tiniest numbers. If you look at the actual implementation, you'll see that we're using  bincoeff to compute this, which in turn computes the binomial coefficients using the usual log-gamma approach.

The bug, if there is one, is in the bincoeff function then. And again, the problem is that we are essentially computing the value of the beta function around a removable singularity, using sums of log-gamma functions. You'll note that on line 92 here,

    http://hg.savannah.gnu.org/hgweb/octave/file/26589abbc78d/scripts/miscellaneous/bincoeff.m#l91

since we've determined that n is not an integer, then we can take the lgamma of it. But then we're actually very close to a singularity of the gamma function. The whole thing still sort of manages to cancel and not give a completely ridiculous answer, but I don't know if there is any sort of code that we can implement here to make it better.

Do you have a suggestion?

Jordi GutiƩrrez Hermoso <jordigh>
Group Member
Thu 01 Aug 2013 12:47:47 AM UTC, comment #4: 

I understand about floating point issues and whether something is an integer or not, but if I want to nbinpdf for real arguments, then I would expect continuity.

octave:2> nbinpdf(0, 3, .5)
ans =  0.12500
octave:3> nbinpdf(0, 2.9999999999999, .5)
ans =  0.12495
octave:4> nbinpdf(0, 1/(1-2/3), .5)
ans =  0.18110           # discontinuous by bad floating pt luck
octave:5> nbinpdf(0, 3.0000000000001, .5)
ans =  0.12475


By definition in terms of the gamma function,

nbinpdf(x,n,p) = gamma(x+n) / (gamma(x+1)*gamma(n)) p^n(1-p)^x

octave:10> x=0; n=1/(1-2/3); p=.5; gamma(x+n) / (gamma(x+1)*gamma(n)) p^n(1-p)^x
ans =  0.12500

octave:11> x=0; n=1/(1-2/3); p=.5; gamma(x+n) / (gamma(x+1)*gamma(n)) p^n(1-p)^x
ans =  0.12500   # gives the right answer


Anonymous
Wed 31 Jul 2013 11:53:02 PM UTC, comment #3: 

No, if that's your reasoning then there is no bug here. The problem is that you think your junk variable is still exactly equal to 3, which it is not.


octave:1> junk = 3;
octave:2> nbinpdf (0, junk, .5)
ans =  0.12500
octave:3> junk = 1 / (1 - 2/3)
junk =  3.0000
octave:4> junk == 3
ans = 0
octave:5> abs (junk - 3)
ans =    4.4409e-16
octave:6> fix (junk)
ans =  2
octave:7> nbinpdf (0, junk, .5)
ans =  0.18110
octave:8> nbinpdf (0, round (junk), .5)
ans =  0.12500


The help for nbinpdf clearly states that it will accept values for N that are either integer or non-integer. It is up to the caller to ensure that the value is actually an integer if that's what you wanted to call it with.

If I misunderstood, you can respond and clarify what you think the bug is, but I am closing this as invalid.

Mike Miller <mtmiller>
Group Member
Wed 31 Jul 2013 11:42:12 PM UTC, comment #2: 

The PMF for negative binomial should be for my example

  nbinpdf(0,3,.5)

((0+3-1) choose 0) (.5)^3 (.5)^0 = .125

see http://en.wikipedia.org/wiki/Negative_binomial

Since 1/(1-2/3) = 1/(1/3) = 3, this should also give the same answer, but it gives .18110 instead.

   nbinpdf(0,1/(1-2/3),.5)

Perhaps the implementation needs to round the second argument to an exact integer?

Anonymous
Wed 31 Jul 2013 07:40:20 PM UTC, comment #1: 

Is this a bug? Ultimately it boils down to discontinuities of the beta function at negative integers, where the gamma function has its poles.

Jordi GutiƩrrez Hermoso <jordigh>
Group Member
Wed 31 Jul 2013 07:10:55 PM UTC, original submission:  


> nbinpdf(0,3,.5)

ans =  0.12500    # this is correct

> junk=1/(1-2/3)  # equals 3 in exact arithmetic

junk =  3.0000

> nbinpdf(0,junk,.5) 

ans =  0.18110    # this is incorrect

Anonymous

 

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

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by siko1056 (Posted a comment)
  • -email is unavailable- added by lachlan (Posted a comment)
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by jordigh (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 14 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-07-28 siko1056 StatusConfirmed Wont Fix
        Open/ClosedOpen Closed
    2019-06-04 mtmiller Summary[Octave-Forge] nbinpdf has floating point error [octave forge] (statistics) nbinpdf has floating point error
    2019-06-04 mtmiller Carbon-CopyRemoved 80942 -
    2019-06-04 rik5 Carbon-CopyRemoved 72865 -
    2019-06-04 rik5 CategoryOctave Function Octave Package
        Summarynbinpdf has floating point error [Octave-Forge] nbinpdf has floating point error
    2016-02-27 lachlan Item GroupIncorrect Result Inaccurate Result
    2013-08-01 mtmiller StatusInvalid / Not an Octave Bug Confirmed
        Open/ClosedClosed Open
        Release3.6.1 dev
        Operating SystemGNU/Linux Any
    2013-07-31 mtmiller StatusNone Invalid / Not an Octave Bug
        Open/ClosedOpen Closed

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code