Re: binocdf inaccuracy in Octave

[Rik]

On 07/08/2013 02:56 PM, Daniel J Sebald wrote:
> Yes, there is a mixing of large and small values such that the small
>> value is swamped by the limited precision of the double type.
>>
>> The example binocdf (1,50, 0.999) translates to
>>
>> 1 - betainc (.999, 2, 49) => 1 - 1 => 0
>>
>> If we don't use the mathematical identity to swap A and B in the inputs
>> to betainc then the expression is
>>
>> betainc (.001, 49, 2)  => 4.99510000000012e-146
>>
>> which is correct.
>>
>> It turns out that the Fortran code in xdbetai.f is making use of the
>> same identity.  Check out this code at the start of the function which
>> swaps A and B and reverses x to 1-x.
>>
>> Y = X
>> P = PIN
>> Q = QIN
>> IF (Q.LE.P .AND. X.LT.0.8D0) GO TO 20
>> IF (X.LT.0.2D0) GO TO 20
>> Y = 1.0D0 - Y
>> P = QIN
>> Q = PIN
>>
>> At the end of the function it sorts things back out again
>>
>> 70   IF (Y.NE.X .OR. P.NE.PIN) DBETAI = 1.0D0 - DBETAI
>>
>> So, when betainc is called as betainc (.999, 2, 49) it reverses the
>> arguments and eventually reaches line 70 where it takes 1.0 - 4.995e-146
>> which equals 1.
>
> I couldn't find that.  All I see in xdbetai.f is:
Sorry, the file is dbetai.f in liboctave/cruft/slatec-fn.

--Rik
>
>       subroutine xdbetai (x, a, b, result)
>       external dbetai
>       double precision x, a, b, result, dbetai
>       result = dbetai (x, a, b)
>       return
>       end
>
>