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


bug #60682: betainc is inaccurate

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

Add a New Comment (Rich Markup)

You are not logged in

Please log in, so followups can be emailed to you.


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 m_ginesi (Updated the item)
  • -email is unavailable- added by mleitner (Submitted the item)

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.




    Follows 1 latest change.

    Date Changed by Updated Field Previous Value => Replaced by
    2021-05-28 m_ginesi Attached File- => Added betainc_trivial.diff, #51492

    Back to the top

    Powered by Savane 3.6