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.

2 comments:

  1. Erwin, 2 questions:
    this post seems to be similar to the Random Number Generator code you built back to 2005, but you added up a statement of the failure of the condition
    abort$(abs(m0 - b0*c0) > 0.001) "calculation of c0 was not correct";

    But since b0=m0/c0, m0-b0*c0 should always be 0, isn't it? how is it possible abs(m0 - b0*c0) > 0.001? I think i've missed some point... could you please help me with that?

    2. it comes up to my mind another question that what if the parameter needed to generate the random number, say, mean or lambda, is a variable rather than a known parameter? can we still generat a random number through this? (in the sense that the random number generated is a function of the unknown parameter, rather than a specific number).Does the code above work?

    ReplyDelete
  2. The test was just to make sure I did not make a mistake in the expression to calculate the scale parameter from a known mean. Nothing to get excited about.

    Point 2 does not make sense to me.

    ReplyDelete