Re: binocdf inaccuracy in Octave

[Pascal Dupuis]

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. 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 ?

Regards

Pascal