bugGNU Octave - Bugs: bug #60528, betaincinv throws error for some...

 
 

bug #60528: betaincinv throws error for some valid inputs

Submitter:  Paul Netsaver <netsaver>
Submitted:  Wed 05 May 2021 10:45:01 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Unexpected Error or Warning
Status:  Fixed Assigned to:  None
Originator Name:  Netsaver Paul 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

Sun 22 Aug 2021 05:58:26 PM UTC, comment #10: 

I used the style


f(finv (x)) == x


for all of the BIST tests except for the ones specifically for this bug report.  From betaincinv.m


%!assert <*60528> (betaincinv (1e-6, 1, 3), 3.3333344444450617e-7, 2*eps)
%!assert <*60528> (betaincinv (1-1e-6, 3, 1), 0.9999996666665555, 2*eps)
%!assert (betainc (betaincinv (0.9, 1e-3, 1), 1e-3, 1), 0.9, 2*eps)


I think it is important that the first BIST test, which is a regression test for the example code in comment #0, be left as a comparison to an exact value.  The second test is not actually a regression test for this bug report and so it can be modified to use the other styling and to drop the bug 60528 tag.

Rik <rik5>
Group administrator
Sun 22 Aug 2021 04:37:13 PM UTC, comment #9: 

Yes, I am also of the opinion that in betaincinv (apart from input validation and the classes of the output) we should test only whether it inverts betainc, thus there should be no magical values anywhere here. And I think that it should even be possible to test for equality.

There is just one point to note that is relevant for the numerical inversion of any generic function f: let's say we start with x. In exact arithmetics, setting y=f(x) we have finv(y)==x. Now the point is that either f or finv (exactly one of the two) will contract a given interval with respect to the local floating-point discretization.

What I mean is that (now in finite-precision arithmetics) if you start with the small interval [x1,x2] that maps to [y1,y2], it can be that there are ten distinct floating-point values between x1 and x2, but only five between y1 and y2. Thus, if you start with one of the ten x, you get one of the five y, and if you then go back through the inverse function, there is no way that you can guarantee to arrive at the correct starting point -- in this example, you have to expect a deviation of one eps, but for more extreme functions, this can easily be much more.

However, what you can do is to choose initially one of the five y. Then you will get an x by finv, and putting it through f again could (and should) give you exactly the initial y (if the inversion works okay, which is what we want to test here).

So I would propose to add here only tests of the kind assert(betaincinv(betainc(x)),x,eps(x)) or the other way around. You can use trial and error to find out which of the two. The tests with hard-coded numerical values should go exclusively into betainc. And yes, one should use test there with respect to the deviations from the closest possible true value (i.e., Wolfram Alpha's 20-digit result, which will be rounded anyway by the interpreter). Sadly, it seems that the continued fraction expansion for general a and b is intrinsically quite inaccurate in this case, so it can indeed be possible that quite large tolerances are required.

Michael Leitner <mleitner>
Sun 22 Aug 2021 10:14:26 AM UTC, comment #8: 

If the goal is to test if Octave's `betainc` and `betaincinv` behave like inverse functions to each other even for these corner case input values, maybe we could add BISTs that do that (additionally)?

If we compare the result of the function to fix values, I'd argue that these values should be what one would expect from a "pure maths" point of view (not including "arbitrary" round off errors). Otherwise, the BISTs might suggest an error when an eventual change at some point in the future actually increases the precision.

Markus Mützel <mmuetzel>
Group administrator
Sat 21 Aug 2021 09:44:19 PM UTC, comment #7: 

Ideally, betaincinv() should produce whatever the pure math answer is.  But, since Octave has limited precision the actual goal is that betainc and betaincinv are inverses of each other such that


x == betainc (becaincinv (x));


For the test case, this holds


octave:21> betainc (betaincinv (1-1e-6, 3, 1, "upper"), 3, 1, "upper")
ans = 0.999999000000000
octave:22> ans - (1-1e-6)
ans = 0


So the issue is that betainc() for both .01 and (.01 + 3.977e-14) both yield 1-1e-6 exactly.  The case of a 'b' input of 1 is a special case and the value for 'upper' tail is


y(b_one) = 1 - x(b_one) .^ a(b_one)
which for this case reduces to
y = 1 - x .^ 3;


Indeed, calculating this result shows that the two values yield the same result.

Rik <rik5>
Group administrator
Sat 21 Aug 2021 07:58:04 PM UTC, comment #6: 

Fwiw, in Matlab R2021a:

>> betaincinv (1-1e-6, 3, 1)

ans =

   0.999999666666556


And according to Wolfram Alpha it should be 0.9999996666665555554938271193415336076582837792003093193653324836...
https://www.wolframalpha.com/input/?i2d=true&i=InverseBetaRegularized%281-1e-6%5C%28%2C%29+3%5C%28%2C%29+1%29

So that looks good.

The second example should be exactly 1/100 according to Wolfram Alpha:
https://www.wolframalpha.com/input/?i2d=true&i=InverseBetaRegularized%281e-6%5C%28%2C%29+3%5C%28%2C%29+1%29

In Matlab R2021a:

>> betaincinv (1-1e-6, 3, 1, "upper")

ans =

   0.010000000000096


So, we seem to be closer than Matlab both on Windows and even more so on Linux.

I updated the reference values from Wolfram Alpha and increased the tolerance for one of the BISTs here:
https://hg.savannah.gnu.org/hgweb/octave/rev/40157299fdeb

Closing report again.

Markus Mützel <mmuetzel>
Group administrator
Sat 21 Aug 2021 06:44:36 PM UTC, comment #5: 

Just increase the tolerance to have the tests pass.  I don't have access to a Win32 machine so I don't know what those values are, but it a few eps are generally fine.  What we've seen repeatedly is that the math libraries on Windows are less accurate than those on Linux.

Rik <rik5>
Group administrator
Sat 21 Aug 2021 05:15:40 PM UTC, comment #4: 

I'm seeing two regressions on Windows with the newly added tests:

>>>>> processing C:\Program Files\GNU Octave\Octave-7.0.0\mingw64\share\octave\7.0.0\m\specfun\betaincinv.m
***** assert <*60528> (betaincinv (1-1e-6, 3, 1), 0.999999666666556, 2*eps)
!!!!! regression: https://octave.org/testfailure/?60528
ASSERT errors for:  assert (betaincinv (1 - 1e-6, 3, 1),0.999999666666556,2 * eps)

  Location  |  Observed  |  Expected  |  Reason
     ()           1            1         Abs err 5.5511e-16 exceeds tol 4.4409e-16 by 1e-16
***** assert <*60528> (betaincinv (1-1e-6, 3, 1,"upper"), .01000000000003977, 2*eps)
!!!!! regression: https://octave.org/testfailure/?60528
ASSERT errors for:  assert (betaincinv (1 - 1e-6, 3, 1, "upper"),.01000000000003977,2 * eps)

  Location  |  Observed  |  Expected  |  Reason
     ()          0.01         0.01       Abs err 1.4074e-14 exceeds tol 4.4409e-16 by 1e-14


The first one is quite close to the tolerance. So, could we just relax it?

The second one is further off from the expected value. (Is this a literature value?):

>> betaincinv (1 - 1e-6, 3, 1, "upper")
ans = 1.000000000005384e-02


Markus Mützel <mmuetzel>
Group administrator
Fri 20 Aug 2021 04:00:12 PM UTC, comment #3: 

The betaincinv function has been completely overhauled.  See the discussion over on bug #60539.  The algorithm is much more accurate than before.  I added a few BIST tests for extreme values in this changeset http://hg.savannah.gnu.org/hgweb/octave/rev/374e68721324.

Marking bug as Fixed and closing report.

Rik <rik5>
Group administrator
Thu 06 May 2021 04:45:01 PM UTC, comment #2: 

On the development branch I further cleaned up the code here: http://hg.savannah.gnu.org/hgweb/octave/rev/6707db623d9d.  Also, see bug #60539 that I filed about the performance of this function.

Rik <rik5>
Group administrator
Wed 05 May 2021 09:31:49 PM UTC, comment #1: 

Confirmed.  The code in betaincinv uses a Newton's Method search.  As part of that process intermediate guesses can fall below zero.  I fixed that in this changeset http://hg.savannah.gnu.org/hgweb/octave/rev/26ba91f0eea7.

The code I used was


  x(x(todo) < 0) = eps;  # Avoid negative x in betainc() call, bug #60528


where I moved any intermediate guess that was negative to be just greater than 0.

Of course, that made me wonder if it was possible for intermediate values to exceed 1.0 and that was also true.  So I had to add another line


  x(x(todo) > 1) = 1-eps;


See additional changeset http://hg.savannah.gnu.org/hgweb/octave/rev/69b6b783a8ab.

Marking as Ready for Test.  You can download the m-file from the Mercurial repository and copy it over your existing betaincinv.m to see the changes.

Rik <rik5>
Group administrator
Wed 05 May 2021 10:45:01 AM UTC, original submission:  

OCTAVE 6.2
betaincinv(1e-6,1,3)
% test printouts from the called 'betainc' follow:
x = 0.0000000, a = 1, b = 3
x = 1.0000000, a = 1, b = 3
x = 0.5000000, a = 1, b = 3
x = 0.2500000, a = 1, b = 3
x = 0.1250000, a = 1, b = 3
x = 0.0625000, a = 1, b = 3
x = 0.0312500, a = 1, b = 3
x = 0.0156250, a = 1, b = 3
x = 0.0078125, a = 1, b = 3
x = 0.0039062, a = 1, b = 3
x = 0.0019531, a = 1, b = 3
x = 0.0009766, a = 1, b = 3
x = 0.0009766, a = 1, b = 3
x = -0.0000006, a = 1, b = 3
error: betainc: X must be in the range [0, 1]
error: called from
    betainc at line 99 column 5
    betaincinv>@<anonymous> at line 172 column 25
    betaincinv>newton_method at line 236 column 15
    betaincinv at line 181 column 14

MATLAB
betaincinv(1e-6,1,3)
ans = 3.333333e-07


out = betaincinv(p,m,n)
It is seen that the internal iterative procedure calls the direct beta function with negative value of p and this cause the error.
It happens with p >= 1e-5 and m=1 and n>2 (not so uncommon values by the way, e.g. p=1e-6, m=1, n=3 may happen in radar applications).
It seems that Matlab intercepts this case and get the results directly dividing p by n: out = p/n.

Regards,
Paul Netsaver





Paul Netsaver <netsaver>

 

(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 mleitner (Posted a comment)
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by rik5
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by netsaver (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 11 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2021-08-21 mmuetzel StatusIn Progress Fixed
        Open/ClosedOpen Closed
    2021-08-21 mmuetzel StatusFixed In Progress
        Open/ClosedClosed Open
    2021-08-20 rik5 StatusReady For Test Fixed
        Open/ClosedOpen Closed
        Release6.2.0 dev
        Carbon-Copy- Added mleitner
    2021-05-05 rik5 Summaryspecial function: betaincinv betaincinv throws error for some valid inputs
    2021-05-05 rik5 StatusNone Ready For Test
        Operating SystemMicrosoft Windows Any

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code