Saturday, May 15, 2010

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:

gammacdf

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.