[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: binocdf inaccuracy in Octave
From: |
Pascal Dupuis |
Subject: |
Re: binocdf inaccuracy in Octave |
Date: |
Mon, 8 Jul 2013 22:36:48 +0200 |
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