Generating random variables having a GAMMA distribution

In GAMS models I often use the following simple algorithm to generate random variable:

1. Generate u = Uniform(0,1)
2. Solve F(x) = u

E.g. for the GAMMA distribution:

$ontext    generate Gamma variates    erwin@amsterdamoptimization.com $offtext set i /i1*i25/; * draw from gamma distribution with mean = 100 scalars   c0  'aka k'      /1.5/   b0  'aka theta'   m0  'mean'       /100/ ; * m=bc b0 = m0/c0; display b0,c0,m0; abort$(abs(m0 - b0*c0) > 0.001) "calculation of c0 was not correct"; parameter u(i); u(i) = uniform(0.001,0.999); display u; variable x(i); equation f(i); f(i).. gammareg(x(i)/b0,c0) =e= u(i); model m /f/; x.lo(i) = 0; x.l(i) = m0; option cns=conopt; solve m using cns; scalar avg; avg = sum(i,x.l(i))/card(i); display avg;

Using a numerical procedure is always dangerous, as we see here:

C O N O P T 3   version 3.14T    Copyright (C)   ARKI Consulting and Development A/S                    Bagsvaerdvej 246 A                    DK-2880 Bagsvaerd, Denmark Using default options. Reading Data   Iter Phase Ninf   Infeasibility   RGmax    NSB   Step InItr MX OK      0   0        6.8475654921E+00 (Input point) ** An equation in the pre-triangular part of the model cannot    be solved because the pivot is too small.    Adding a bound or initial value may help.

The CDF of the Gamma distribution looks quite innocent:

Two simple ways of working around this. First PATH seems to handle this a little bit better (and CONOPT agrees with the solution). This can be done with:

option cns=path; solve m using cns; option cns=conopt; solve m using cns;

Indeed CONOPT new says:

C O N O P T 3   version 3.14T    Copyright (C)   ARKI Consulting and Development A/S                    Bagsvaerdvej 246 A                    DK-2880 Bagsvaerd, Denmark Using default options. Reading Data   Iter Phase Ninf   Infeasibility   RGmax    NSB   Step InItr MX OK      0   0        7.9597856362E-09 (Input point)                                Pre-triangular equations:       25                                Post-triangular equations:       0      1   0        0.0000000000E+00 (After pre-processing)      2   0        0.0000000000E+00 (After scaling)      2            0.0000000000E+00 ** Feasible solution to a square system.

Another solution would be to use a different starting point. We used the mean which is an obvious starting point, but with

x.l(i) = 1;

instead of

x.l(i) = m0;

CONOPT can also solve the equations directly.

This illustrates that, although we can work around things, it would be better to have a larger number of statistical functions in GAMS, including some inverse CDF’s.