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 |
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 |
No comments:
Post a Comment