Re: binocdf accuracy

[CdeMills]

Dr. Alexander Klein-2 wrote
> One of my experiments included binocdf (0:50,50,1-1e-3 )', and while I
> admit that the example is somewhat extreme, I suspect that almost all
> values returned by Octave are inaccurate, while binopdf (0:50,50,1-1e-3 )'
> produces results nearly identical to my implementation in Erlang, see
> attachments.
> 
> I think this would be easy to patch if only binocdf was afflicted, but
> binocdf uses the incomplete beta function, and so maybe other functions
> might have problems, too.

Hello again,

there are two sources of numerical inaccuracies in the betainc function:
1) the use of bincoeff
2) the use of power

I added a betainc implementation to the multi-precision package where those
two operations are handled with extended numerical resolution. In your case,
bincoeff(50, something) maximum value is factorial(50)/factorial(25)^2.

I enclose a file with the 'exact' values, i.e. computed with 600 bits of
precision; this was required to get the smallest values:

x=[0:49]; n= 50;p=1-mp(1, 600)/1000;

the important point is the probability value: you use 1-1e-3; which is not a
number with finite binary representation. So I take the integer 1, convert
it to a mp number on 600 bits; divide by the integer 1000 and substract this
from one.

 cdfmp=1-betainc(p, 1+floor(x), n-floor(x));

in my implementation, all the computations are performed at the precision of
the input with the highest one, i.e. 600

To save to a file:
fid=fopen('binocdf-mp.txt', 'wt');
fprintf(fid, '%.15e\n', double(cdfmp.'));
fclose(fid);

Regards

Pascal

binocdf-mp.txt
<http://octave.1599824.n4.nabble.com/file/n4655283/binocdf-mp.txt>  





--
View this message in context: 
http://octave.1599824.n4.nabble.com/binocdf-accuracy-tp4655269p4655283.html
Sent from the Octave - Maintainers mailing list archive at Nabble.com.