binocdf function apparently not functioning with very small numbers (near 1e-15)

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

Answer 1

[Roger Stafford]

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.

Answer 2

[Andrew]

Bump... Can someone help out with this? Does such a small level of alpha produce precision issues? Or maybe there is something else wrong?

Answer 3

[Andrew]

Hi Roger, thanks for your answer and explanation of what's going on with MATLAB's double precision. When I try running with 'upper' switch, I get:

* Error using binocdf Too many input arguments. *

Moreover, my MATLAB (2012A) documentation does not contain the 'upper' option. It looks like this was added more recently (perhaps with 2014A?).