[Octave-bug-tracker] [bug #60528] betaincinv throws error for some valid inputs

[Michael Leitner]

Follow-up Comment #9, bug #60528 (project octave):

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.

    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?60528>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/