computing binomial distribution tails

[Johan Kullstam]

I wish to compute the tail of a binomial distribution but I keep
running into numerical difficulties.

For example, given probability p (small) what is the chance of no
success for m bernouli trials?

Q = 1-(1-p)^m

has trouble when p is near eps in magnitude and fails utterly when p <
eps.

using

Q = 1-binomial(0,m,p)

doesn't help.

In C, using glibc, I note that there are log1p and expm1 functions
which allows

Q = -expm1(m*log1p(-p))

and this works just fine, but, alas, there are no log1p or expm1
functions in octave.  (I did gin up a couple of C++ octfiles though.)

Another related problem is computing the binomial_cdf tail when p is
small.

binomial_cdf(k,m,p) is very near to 1 and hence what I really want is
a kind of complementary cdf like erfc is to erf.  Does such a thing
exist for the binomial distribution in octave?

I have been resorting to

s = 0;
for ii = m:-1:k+1,
  ## accumulate the small fry first
  s += binomial_pdf(ii,m,p);
endfor

instead of
s = 1-binomial_cdf(k,m,p);

but it is rather slow when m-k is large (greater than say 500).

-- 
Johan KULLSTAM


-- 
Johan KULLSTAM



-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------