bugGNU Octave - Bugs: bug #60682, betainc is inaccurate


bug #60682: betainc is inaccurate

Submitter:  Michael Leitner <mleitner>
Submitted:  Fri 28 May 2021 08:51:45 AM UTC
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * GNU/Linux
* Mandatory Fields

Add a New Comment (Rich Markup)

Comment Type & Canned Response:


( Jump to the original submission )

Mon 18 Apr 2022 09:03:44 PM UTC, comment #7: 

This change has introduced a new problem reported in bug #62329.

Rik <rik5>
Project Administrator
Thu 19 Aug 2021 09:09:46 PM UTC, comment #6: 

I made some modifications for speed and then checked in Michele's code here: http://hg.savannah.gnu.org/hgweb/octave/rev/ad6a57b215e8.

Michael was in favor of adopting code for these special (although important) cases if the overhead was not too large.

I benchmarked using a 1M point vector generated with "rand (1e6, 1)".  Hence, no value of input 'a' or 'b' was actually equal to 1.  Thus, original code and new code were both using continued fraction expansion for calculation of output, and differences in timing were solely down to the extra input processing.

For original code, mean run time (N = 25) was 380 milliseconds.

For Michele's code shown below, mean run time was 396 milliseconds.

  a_b_one = (a == 1) & (b == 1);
  a_one = (a == 1) & (b != 1);
  b_one = (a != 1) & (b == 1);
  non_trivial = (a != 1) & (b != 1);

An increase of 4.2% in run time is probably okay, but I decided to re-code to trade memory for speed since most machines have copious memory these days.

New code:

  ## Trivial cases (long code here trades memory for speed)
  a_one = (a == 1);
  b_one = (b == 1);
  a_b_one = a_one & b_one;
  a_not_one = ! a_one;
  b_not_one = ! b_one;
  non_trivial = a_not_one & b_not_one;
  a_one &= b_not_one;
  b_one &= a_not_one;

Mean run time is now 382 milliseconds.  Percentage increase is just 0.5% so I checked in the change.

Marking as Fixed and closing report.

Rik <rik5>
Project Administrator
Thu 12 Aug 2021 07:35:55 AM UTC, comment #5: 

As to the comparison between matlab and octave without the fix: it seems that matlab is a bit more accurate (the jump in the result around 2^-18 seems to be only have as large, and the deviations for  x equal two to the power of an integer are on average somewhat smaller), but there is no qualitative difference, I would say. Yes, the fix works for b equal to one.

But I don't understand your last remark "I'm not seeing the accuracy issue".

Michael Leitner <mleitner>
Wed 11 Aug 2021 11:47:07 PM UTC, comment #4: 

curious how matlab compares, here are the outputs of the two commands below in Matlab 2021a:

>> x=2.^(-(0:20))';[x x.^2 betainc(x,2,1)./x.^2-1]

ans =

     1.000000000000000e+00     1.000000000000000e+00                         0
     5.000000000000000e-01     2.500000000000000e-01     2.220446049250313e-16
     2.500000000000000e-01     6.250000000000000e-02     2.220446049250313e-16
     1.250000000000000e-01     1.562500000000000e-02     2.220446049250313e-16
     6.250000000000000e-02     3.906250000000000e-03     2.220446049250313e-16
     3.125000000000000e-02     9.765625000000000e-04    -4.440892098500626e-16
     1.562500000000000e-02     2.441406250000000e-04     8.881784197001252e-16
     7.812500000000000e-03     6.103515625000000e-05    -1.110223024625157e-16
     3.906250000000000e-03     1.525878906250000e-05     1.110223024625157e-15
     1.953125000000000e-03     3.814697265625000e-06    -7.771561172376096e-16
     9.765625000000000e-04     9.536743164062500e-07    -4.440892098500626e-16
     4.882812500000000e-04     2.384185791015625e-07     4.440892098500626e-16
     2.441406250000000e-04     5.960464477539063e-08                         0
     1.220703125000000e-04     1.490116119384766e-08                         0
     6.103515625000000e-05     3.725290298461914e-09     2.220446049250313e-15
     3.051757812500000e-05     9.313225746154785e-10    -2.220446049250313e-16
     1.525878906250000e-05     2.328306436538696e-10     1.776356839400250e-15
     7.629394531250000e-06     5.820766091346741e-11     6.661338147750939e-16
     3.814697265625000e-06     1.455191522836685e-11     4.440892098500626e-16
     1.907348632812500e-06     3.637978807091713e-12     2.220446049250313e-16
     9.536743164062500e-07     9.094947017729282e-13                         0

>> x=2.^-18*(1+eps(1)*(0:9)');[x x.^2 betainc(x,2,1) betainc(x,2,1)./x.^2-1]

ans =

     3.814697265625000e-06     1.455191522836685e-11     1.455191522836686e-11     4.440892098500626e-16
     3.814697265625001e-06     1.455191522836686e-11     1.455191522836686e-11                         0
     3.814697265625002e-06     1.455191522836686e-11     1.455191522836686e-11    -4.440892098500626e-16
     3.814697265625003e-06     1.455191522836687e-11     1.455191522836686e-11    -8.881784197001252e-16
     3.814697265625003e-06     1.455191522836688e-11     1.455191522836686e-11    -1.332267629550188e-15
     3.814697265625004e-06     1.455191522836688e-11     1.455191522836691e-11     1.776356839400250e-15
     3.814697265625005e-06     1.455191522836689e-11     1.455191522836691e-11     1.332267629550188e-15
     3.814697265625006e-06     1.455191522836690e-11     1.455191522836691e-11     8.881784197001252e-16
     3.814697265625007e-06     1.455191522836690e-11     1.455191522836691e-11     4.440892098500626e-16
     3.814697265625008e-06     1.455191522836691e-11     1.455191522836691e-11                         0

using Octave 6.3.0:

>> x=2.^(-(0:20))';[x x.^2 betainc(x,2,1)./x.^2-1]
ans =

   1.000000000000000e+00   1.000000000000000e+00                       0
   5.000000000000000e-01   2.500000000000000e-01   2.220446049250313e-16
   2.500000000000000e-01   6.250000000000000e-02                       0
   1.250000000000000e-01   1.562500000000000e-02   4.440892098500626e-16
   6.250000000000000e-02   3.906250000000000e-03   2.220446049250313e-16
   3.125000000000000e-02   9.765625000000000e-04                       0
   1.562500000000000e-02   2.441406250000000e-04   6.661338147750939e-16
   7.812500000000000e-03   6.103515625000000e-05   4.440892098500626e-16
   3.906250000000000e-03   1.525878906250000e-05   4.440892098500626e-16
   1.953125000000000e-03   3.814697265625000e-06   2.220446049250313e-16
   9.765625000000000e-04   9.536743164062500e-07                       0
   4.882812500000000e-04   2.384185791015625e-07  -1.110223024625157e-16
   2.441406250000000e-04   5.960464477539062e-08   1.554312234475219e-15
   1.220703125000000e-04   1.490116119384766e-08  -2.331468351712829e-15
   6.103515625000000e-05   3.725290298461914e-09   1.110223024625157e-15
   3.051757812500000e-05   9.313225746154785e-10  -2.664535259100376e-15
   1.525878906250000e-05   2.328306436538696e-10   6.661338147750939e-16
   7.629394531250000e-06   5.820766091346741e-11   6.661338147750939e-16
   3.814697265625000e-06   1.455191522836685e-11  -3.108624468950438e-15
   1.907348632812500e-06   3.637978807091713e-12   2.220446049250313e-16
   9.536743164062500e-07   9.094947017729282e-13                       0

>> x=2.^-18*(1+eps(1)*(0:9)');[x x.^2 betainc(x,2,1) betainc(x,2,1)./x.^2-1]
ans =

   3.814697265625000e-06   1.455191522836685e-11   1.455191522836681e-11  -3.108624468950438e-15
   3.814697265625001e-06   1.455191522836686e-11   1.455191522836681e-11  -3.552713678800501e-15
   3.814697265625002e-06   1.455191522836686e-11   1.455191522836681e-11  -3.996802888650564e-15
   3.814697265625003e-06   1.455191522836687e-11   1.455191522836681e-11  -4.440892098500626e-15
   3.814697265625003e-06   1.455191522836688e-11   1.455191522836681e-11  -4.884981308350689e-15
   3.814697265625004e-06   1.455191522836688e-11   1.455191522836691e-11   1.776356839400250e-15
   3.814697265625005e-06   1.455191522836689e-11   1.455191522836691e-11   1.332267629550188e-15
   3.814697265625006e-06   1.455191522836690e-11   1.455191522836691e-11   8.881784197001252e-16
   3.814697265625007e-06   1.455191522836690e-11   1.455191522836691e-11   4.440892098500626e-16
   3.814697265625008e-06   1.455191522836691e-11   1.455191522836691e-11                       0

I thought maybe the patch was already pushed but

Octave 6.3.0:

>> betainc (0.5, 1, Inf)
ans = NaN

after applying the patch i get:

>> x=2.^(-(0:20))';[x x.^2 betainc(x,2,1)./x.^2-1]
ans =

   1.000000000000000e+00   1.000000000000000e+00                       0
   5.000000000000000e-01   2.500000000000000e-01                       0
   2.500000000000000e-01   6.250000000000000e-02                       0
   1.250000000000000e-01   1.562500000000000e-02                       0
   6.250000000000000e-02   3.906250000000000e-03                       0
   3.125000000000000e-02   9.765625000000000e-04                       0
   1.562500000000000e-02   2.441406250000000e-04                       0
   7.812500000000000e-03   6.103515625000000e-05                       0
   3.906250000000000e-03   1.525878906250000e-05                       0
   1.953125000000000e-03   3.814697265625000e-06                       0
   9.765625000000000e-04   9.536743164062500e-07                       0
   4.882812500000000e-04   2.384185791015625e-07                       0
   2.441406250000000e-04   5.960464477539062e-08                       0
   1.220703125000000e-04   1.490116119384766e-08                       0
   6.103515625000000e-05   3.725290298461914e-09                       0
   3.051757812500000e-05   9.313225746154785e-10                       0
   1.525878906250000e-05   2.328306436538696e-10                       0
   7.629394531250000e-06   5.820766091346741e-11                       0
   3.814697265625000e-06   1.455191522836685e-11                       0
   1.907348632812500e-06   3.637978807091713e-12                       0
   9.536743164062500e-07   9.094947017729282e-13                       0

>> x=2.^-18*(1+eps(1)*(0:9)');[x x.^2 betainc(x,2,1) betainc(x,2,1)./x.^2-1]
ans =

   3.814697265625000e-06   1.455191522836685e-11   1.455191522836685e-11                       0
   3.814697265625001e-06   1.455191522836686e-11   1.455191522836686e-11                       0
   3.814697265625002e-06   1.455191522836686e-11   1.455191522836686e-11                       0
   3.814697265625003e-06   1.455191522836687e-11   1.455191522836687e-11                       0
   3.814697265625003e-06   1.455191522836688e-11   1.455191522836688e-11                       0
   3.814697265625004e-06   1.455191522836688e-11   1.455191522836688e-11                       0
   3.814697265625005e-06   1.455191522836689e-11   1.455191522836689e-11                       0
   3.814697265625006e-06   1.455191522836690e-11   1.455191522836690e-11                       0
   3.814697265625007e-06   1.455191522836690e-11   1.455191522836690e-11                       0
   3.814697265625008e-06   1.455191522836691e-11   1.455191522836691e-11                       0

>> betainc (0.5, 1, Inf)
ans = 1.000000000000000e+00

the trivial case fix is a good fix.  but I'm not seeing the accuracy issue

Nicholas Jankowski <nrjank>
Project Member
Fri 28 May 2021 07:43:11 PM UTC, comment #3: 

It seems that it is just not possible to compute betainc with double precision using the classical continued fraction to more than a few eps accuracy -- I was thinking that it would become more accurate by evaluating the continued fraction directly as opposed to something like Lentz' algorithm, but it does not seem to be the case. There would however be ways so that the computation finishes after l=b rounds if the parameter b is an integer, which supposedly would give a more accurate answer specifically for small b, but it is questionable whether it is worth the effort, both with respect to the programmer's as well as the computer's.

If Michele's fix works and does not introduce measureable overhead, I would say that is about the best we can do.

Michael Leitner <mleitner>
Fri 28 May 2021 04:43:54 PM UTC, comment #2: 

That's probably a good idea, as I can imagine that quite often such trivial a and b will be used. However, the inaccuracy exists also in the general case (it's only that I then don't have a way do demonstrate it quickly in octave), so this is only a partial fix.

Michael Leitner <mleitner>
Fri 28 May 2021 04:36:54 PM UTC, comment #1: 

This can be easily solved by computing the three cases (a=1,b=1), (a=1, b!=1) and (a!=1, b=1) separately and using the continued fractions only on non-trivial cases.
I attach the diff implementing this change.

Coincidentally, test <*34405> now would fail. However, the new behavior is Matlab compatible. Before the change

betainc (0.5, 1, Inf)
>> NaN


betainc (0.5, 1, Inf)
>> 1

which is the same answer given by Matlab (2020a)

(file #51492)

Michele Ginesi <m_ginesi>
Fri 28 May 2021 08:51:45 AM UTC, original submission:  

For the second parameter b equal to one, betainc(x,a,1) should be equal to x^a. Double precision corresponds to eps(1)=2.2204e-16, so the third column of the output when doing

format long e
x=2.^(-(0:20))';[x x.^2 betainc(x,2,1)./x.^2-1]

should not be much larger than this value. Note that for x equal 2 to the power of an integer and also a an integer, x^a is exactly representable, so it could also be zero. However, in this display of not so extreme values I see deviations of up to 14 eps, namely at x=2^-18, and if you go to smaller x or larger a, you get even larger deviations.

Further, if you do

format long e
x=2.^-18*(1+eps(1)*(0:9)');[x x.^2 betainc(x,2,1) betainc(x,2,1)./x.^2-1]

you see that indeed the second column is continuously growing, but the third shows only a single jump.

All of the numbers I report here were obtained on /octave-online.net.

I did not look into the code, but betainc will be computed by some series expansion or continued fraction, and it seems that the aborting condition of this computation is a bit too tight.

Michael Leitner <mleitner>


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

Attach Files:

Attached Files
file #51492:  betainc_trivial.diff added by m_ginesi (3KiB - text/x-patch)


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by m_ginesi (Updated the item)
  • -email is unavailable- added by mleitner (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 project members can vote.





    Follow 4 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2021-08-19 rik5 StatusNone Fixed
        Open/ClosedOpen Closed
        Release6.1.0 dev
    2021-05-28 m_ginesi Attached File- Added betainc_trivial.diff, #51492

    Back to the top

    Powered by Savane 3.9