*!version 1.1 1999 Joseph Hilbe
* version 1.0.0  1993 Joseph Hilbe, Walter Linde-Zwirble        (sg44: STB-28)
* gamma distribution random number generator with shape & scale variables
* Example: rndgamx mu, s(1) 

program define rndgamx
  version 4.0
  set type double
  qui  {
    local varlist "req ex"
    local options "`options' Shape(real 1.0)"
    parse "`*'"
    parse "`varlist'",parse(" ")
    tempvar ia theta mu
    gen `mu'= `1'
    gen `ia' = `shape'
    gen `theta'=`mu'*`ia'  /* scale parameter */
    noi di in gr "( Generating " _c
    if `ia'<1 { error }
    if (`ia'<6 & `ia'==int(`ia')) {   /* if shape 1-5 */
       tempvar x
       gen `x' = 1.0
       local i = 1
       while `i'<=`ia'  {
          replace `x'=`x'*uniform()
          local i = `i'+1
       noi di in gr "." _c
       }
       replace `x' = -log(`x')
    }
    else {   
       tempvar xr e ran1 ran2 ds ts sum1 y x
       local am = `ia'-1
       local s = sqrt(2.0*`am'+1.0)
       gen `x' = 1.0
       gen `xr' = -1.0
       gen `e'  = -1.0
       gen `ran1'= uniform()
       gen `ran2'= uniform()
       gen `ds'= 1
       gen `ts'= 1
       gen `y' = -1
       egen `sum1' = sum(`ds')
       while `sum1'>0   {
          replace `y'=sin(_pi*`ran1')/cos(-_pi*`ran1')
          replace `xr'=`s'*`y' + `am' if (`ds'==1)
          replace `ts'=0 if (`xr' <= 0 & `ds'==1)
          #delimit ;
          replace `e' = (1.0+`y'*`y')*
           (exp(`am'*log(`xr'/`am')-(`s'*`y'))) if (`ds'==1 & `ts'==1);
          #delimit cr
          replace `ds'=0 if (`e' >= `ran2' & `ds'==1 & `ts'==1)
          replace `ran1'= uniform()
          replace `ran2'= uniform()
          replace `e' = -1
          replace `ts' = 1
          drop `sum1'
          egen `sum1' = sum(`ds')
          noi di in gr "." _c
       }
       replace `x' = `xr'
  }
  noi di in gr " )"
  gen xg = `x' * `theta'
  noi di in bl "Variable " in ye "xg " in bl "created."
  lab var xg "Constructed gamma random variable"
  set type float
}
end