Binomial

In a nonlinear model I needed to use the following expression in a model equation:

\[y = \binom{x}{k}p^k(1-p)^{(x-k)}\]

Here \(x\) is the only decision variable. The integer version of “n choose k” is well known:

\[\binom{n}{k} = \frac{n!}{k!(n-k)!}\]

For fractional values we can generalize this to:

\[\binom{x}{y} = \frac{\Gamma(x+1)}{\Gamma(y+1)\Gamma(x-y+1)}\]

where \(\Gamma(.)\) indicates the gamma function, defined by

\[\Gamma(z) = \int_0^{\infty}x^{z-1}e^{-x} dx\]

In GAMS we have the binomial(x,y) function for this. The problem is that \(x\) can become rather large, while \(k\) is a small integer, say \(k=0,\dots,4\). For some values of \(x\) this leads to very inaccurate results and even overflows in subsequent calculations.

E.g.: 

scalars   x /112.8/   k /2/ ; parameter y; y("binomial") = binomial(x,k); y("log") = exp(loggamma(x+1)-loggamma(k+1)-loggamma(x-k+1)); display y;

yields the wrong result:

----      8 PARAMETER y  binomial 1.0000E+299,    log         6305.520

The logarithmic version is correct.

Instead of using the logarithmic version:

\[\binom{x}{k} = \exp\left(\log\Gamma(x+1)-\log\Gamma(k+1)-\log\Gamma(x-k+1)\right)\]

we can also exploit that \(k\) is a small integer. From:

\[\Gamma(z+1)=z\Gamma(z)\]

we have:

\[\begin{align}\binom{x}{k}&=\frac{x(x-1)\cdots(x-k+1)\Gamma(x-k+1)}{k!\Gamma(x-k+1)}\\ &=\frac{x(x-1)\cdots(x-k+1)}{k!}\\ &=\frac{\displaystyle\prod_{i=1}^k (x-i+1)}{k!}\end{align}\]

In GAMS:

scalars   x /112.8/   k /2/ ; parameter y; y("binomial") = binomial(x,k); y("log") = exp(loggamma(x+1)-loggamma(k+1)-loggamma(x-k+1)); set i /1*2/; y("prod") = prod(i,x-i.val+1)/fact(k); display y;

with results:

----     10 PARAMETER y  binomial 1.0000E+299,    log         6305.520,    prod        6305.520

I used this last “product” version in the model.

Note that the overflow in the binomial function does not happen for all (or even most) large values:

scalars   x /1112.8/   k /2/ ; parameter y; y("binomial") = binomial(x,k); y("log") = exp(loggamma(x+1)-loggamma(k+1)-loggamma(x-k+1)); set i /1*2/; y("prod") = prod(i,x-i.val+1)/fact(k); display y;

gives:

----     10 PARAMETER y  binomial 618605.520,    log      618605.520,    prod     618605.520