Re: binocdf inaccuracy in Octave

[Daniel J Sebald]

On 07/08/2013 03:36 PM, Pascal Dupuis wrote:
> 2013/7/8 Daniel J Sebald<address@hidden>:
> > Rik, that string of zeros at the front of binocdf () has me wondering if the
> > issue here is that the FORTRAN routine computing betainc () is single
> > precision float and not double precision float.  Note that the binopdf.m
> > script uses the following computations:
> >
> >      pdf(k) = exp (gammaln (n+1) - gammaln (x(k)+1) - gammaln (n-x(k)+1)
> >                    + x(k)*log (p) + (n-x(k))*log (1-p));
> >
> > Are these all double precision operations?
>
>
> In file liboctave/cruft/slatec-fn/dbetai.f, all variables are double
> precision.

Yes, I saw that too.  However, exactly what is going on in source code, 
I'm not sure.  I searched for the code in SLATEC but couldn't find it.


>  The standard way to avoid accuracy loss when summing a lot
> of terms spanning many order of magnitude, is to sum smallest
> magnitude terms first and terminates by the greatest terms. So maybe
> this routine is optimised for small p and has issues with p close to 1
> as the magnitudes of the terms are reversed ?  We could try one of the
> equivalence formula to transform p close to one into p close to 0 ?

Good point.  Perhaps that is what the issue is.  This:

cdf(k) = 1 - betainc (p, tmp + 1, n - tmp);

is subtracting a numerically tiny value from 1.0

> > 1 - 1e-50
ans =  1
> > 1 - 1e-45
ans =  1
> > 1 - 1e-40
ans =  1
> > 1 - 1e-35
ans =  1
> > 1 - 1e-30
ans =  1
> > 1 - 1e-25
ans =  1
> > 1 - 1e-20
ans =  1
> > 1 - 1e-15
ans =  0.999999999999999
> > 1 - 1e-12
ans =  0.999999999999000
> > 1 - 1e-9
ans =  0.999999999000000
> > 1 - 1e-6
ans =  0.999999000000000
> > 1 - 1e-3
ans =  0.999000000000000
> > 1 - 1e-0
ans = 0

So, although betainc() can converge to anything within the range [0,1] 
(let's assume for the moment the core library routine is always doing 
something proper as far as accuracy), it's numerical value can be so 
small that arithmetic on the results is dodgy.

[Note, I'm going to write a follow-up thread about "eps", numerical 
accuracy.]

So, we should never really be doing any such kind of

ret = 1 - numlibfunc (...)

However, do a search

  [sebald@ scripts]$ grep " 1 -" * -r

and the scripts, at least, are littered with all sorts of operations.

Rik, you really discovered something here...looks like you have your 
work cut out for you!  :-)

Dan