[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: binocdf accuracy
From: |
CdeMills |
Subject: |
Re: binocdf accuracy |
Date: |
Sun, 7 Jul 2013 04:38:53 -0700 (PDT) |
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.