how to compute binomial distribution without overflow?

[gracia]

I need to compute the value of binomial(n, k)=n!/k!(n-k)! * (1-p)^n *
p^k.

When n and k is very big (e.g. n > 160), the step to compute n!/k!(n-
k)! is always overflow even when I used unsigned long long.

Can anyone help me?

Thank you


Gracia

 

[tommy.hinks]

I need to compute the value of binomial(n, k)=n!/k!(n-k)! * (1-p)^n *
p^k.

When n and k is very big (e.g. n > 160), the step to compute n!/k!(n-
k)! is always overflow even when I used unsigned long long.

Can anyone help me?

Thank you

Gracia

Did you try using double or long double?

I am no expert in this field, but it seems you might be able to cancel
some of the terms:

e.g.

n = 10
k = 5

1*2*3*4*5*6*7*8*9*10
--------------------- = 6*7*8*9*10
1*2*3*4*5

T

 

[Christopher]

Did you try using double or long double?

I am no expert in this field, but it seems you might be able to cancel
some of the terms:

e.g.

n = 10
k = 5

1*2*3*4*5*6*7*8*9*10
--------------------- = 6*7*8*9*10
1*2*3*4*5

T

If you do not have enough bits you can always roll your own data type.
Lots of people have done a big int class, you could do a big float or
big double or search for one. Your doing math anyway, it could be
fun Either way you will have to decide before hand what your domain
and range are, make sure you can store them, and check for overflow.

 

[Michal Spalinski]

I need to compute the value of binomial(n, k)=n!/k!(n-k)! * (1-p)^n *
p^k.

When n and k is very big (e.g. n > 160), the step to compute n!/k!(n-
k)! is always overflow even when I used unsigned long long

You can use a recurrence relation which is equivalent to the definition you
cite:

long bino(int n, int k)
{
if (k==0 || n == k) return 1;
return bino(n-1,k-1) + bino(n-1,k);
}

You can speed this up considerably by caching.

 

[osmium]

I need to compute the value of binomial(n, k)=n!/k!(n-k)! * (1-p)^n *
p^k.

When n and k is very big (e.g. n > 160), the step to compute n!/k!(n-
k)! is always overflow even when I used unsigned long long.

The easiest, but probably not the best in terms of accuracy and speed, way
would be to use logarithms.

 

[Kai-Uwe Bux]

I need to compute the value of binomial(n, k)=n!/k!(n-k)! * (1-p)^n *
p^k.

I think the (1-p) exponent should be n-k.

When n and k is very big (e.g. n > 160), the step to compute n!/k!(n-
k)! is always overflow even when I used unsigned long long.

Can anyone help me?

What about:

double binomial ( unsigned int n,
unsigned int k,
double p ) {
assert( k <= n );
double result = 1;
double q = 1 - p;
if ( k > n/2 ) {
k = n - k;
swap( p, q );
}
unsigned int l;
for ( l = 1; l <= k; ++l ) {
result *= p*q*(k+l)/l;
}
for ( l = k+1; l <= n-k; ++l ) {
result *= q*(l+k)/l;
}
return ( result );
}

(not tested)


Best

Kai-Uwe Bux

 

[James Kanze]

I need to compute the value of binomial(n, k)=n!/k!(n-k)! *
(1-p)^n * p^k.

When n and k is very big (e.g. n > 160), the step to compute
n!/k!(n- k)! is always overflow even when I used unsigned long
long.

Can anyone help me?

If you have access to the C99 math functions (which will be part
of the next version of C++ as well), you can use something like:

double
fact(
double i )
{
return round( exp( gamma( (double)( i + 1 ) ) ) ) ;
}

long
binom(
long i,
long j )
{
return round( fact( i ) / ( fact( j ) * fact( i - j ) ) ) ;
}

On the other hand, I'm not sure that it will be really correct
for very large values. Typically, a double will not be able to
represent all of the needed values exactly.

 

[Lionel B]

The easiest, but probably not the best in terms of accuracy and speed,
way would be to use logarithms.

Agreed, if n is reasonably large. I use the following code, which uses the
`double lgamma(double x)' which computes the logarithm of the gamma function.
I don't think it's standard C++, but there's probably a good chance you have
it, or something equivalent (on my linux system, the docs tell me it conforms
to C99, SVr4 and 4.3BSD).

inline double lfactorial(const unsigned n)
{
return lgamma(static_cast<double>(n+1));
}

inline double lbicof(const unsigned n, const unsigned k)
{
assert(k <= n);
return lfactorial(n)-lfactorial(k)-lfactorial(n-k);
}

inline double bicof(const unsigned n, const unsigned k)
{
return std::exp(lbicof(n,k));
}

inline double lbinom(const unsigned n, const unsigned k, const double p)
{
assert(p >= 0.0 && p <= 1.0);
return lbicof(n,k) + static_cast<double>(k)*std::log(p) + static_cast<double>(n-k)*std::log(1.0-p);
}

inline double binom(const unsigned n, const unsigned k, const double p)
{
return std::exp(lbinom(n,k,p));
}

 

[iandjmsmith]

I need to compute the value ofbinomial(n, k)=n!/k!(n-k)! * (1-p)^n *
p^k.

When n and k is very big (e.g. n > 160), the step to compute n!/k!(n-
k)! is always overflow even when I used unsigned long long.

Can anyone help me?

Thank you

Gracia

See
http://groups.google.co.uk/group/comp.lang.javascript/msg/ebaefcd0d790545b?hl=en&lr=&rnum=8
for details of how to perform accurate calculations for the binomial
distribution. VBA code for these calculations can be found in
http://members.aol.com/iandjmsmith/Examples.xls

A Javascript binomial distribution calculator can be found at
http://members.aol.com/iandjmsmith/EXAMPLES.HTM


Ian Smith