binocdf function apparently not functioning with very small numbers (near 1e-15)
6 views (last 30 days)
Show older comments
Having difficulty debugging what's going on with the following code. I'm wondering whether I need to use a special variable or package, given the very small magnitude of the variable values?
Basically, the code is to compute the k-tail of the (1/(s!)) binomial distribution, namely:
\sum_{i=k}^n (1/(s!))^i (1-(1/(s!)))^(n-1)
After taking in a fixed n, m, and statistical significance level alpha (and it's correct that I'm desiring VERY small values here), I'm looking for the first level k for which the value usk drops below alpha.
When I get to a value of s = 9, for example, the code should exit around k = 5 or so, but does not (it reaches 612, as the value of usk seems stuck around 1e-14 > 1e-15).
Can someone point me in the right direction? It does work for smaller values of s.
I'm using MATLAB R2012A, but do have access to R2013A on Unix (though both are giving the same strange result).
Thanks!
function StatSig
n = 612;
m = 18;
alph = 1e-15;
%Starting at level s = 3, go through each possible level s, up to m
for s=3:m
%Re-initialize variable k
k = 2;
%Keep incrementing the number of rows k until we drop below the level of statistical significance alpha
keep_going = true;
while (keep_going)
usk = 1 - binocdf(k,n,(1/factorial(s)));
%if condition is met...
if (usk < alph)
keep_going = false;
else
k = k + 1;
end
end
fprintf('%d\n', k);
end
end
0 Comments
Accepted Answer
Roger Stafford
on 21 Aug 2014
Your problem is one of accuracy at the upper end. You have written:
usk = 1 - binocdf(k,n,(1/factorial(s)));
so that a single bit error in binocdf down from an exact one would be 1-2^(-53), and when you subtract this from 1, you get 2^(-53) = 1.110e-16. You are expecting more of matlab's double precision than it can deliver.
You should be taking advantage of the 'upper' option of 'binocdf' where you won't face this issue. With that you don't subtract the result from 1 and consequently get much great relative accuracy. Read the documentation.
0 Comments
More Answers (2)
See Also
Categories
Find more on Matrix Indexing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!