*!version 1.1 1999 Joseph Hilbe
* version 1.0.0 04/14/95 Joseph Hilbe                       (sg44: STB-28)
* Binomial distribution random number generator
* EX: rndbinx mu den : mu= variable with p values; den=case denominator
program define rndbinx
  version 4.0
  set type double
  local varlist "ex req"
  parse "`*'"
  parse "`varlist'",parse(" ")
  local pp "`1'"
  mac shift
  local n "`1'"
  mac shift
  qui {
  tempvar m am ran1 ran2 ds ts sum1 e y en em bnl mu pc p plog pclog oldg 
  tempvar s t j g
  gen `p'=`pp'
  replace `p'=1.0-`pp' if `pp' > 0.5
  gen `am'=`n'*`p'
  gen `m'=3
  replace `m'=1 if `n' < 25
  replace `m'=2 if ((`am'<1.0) & (`n'>=25))
  gen `mu'=`n'*`p'
  gen `pc'=1.0-`p'
  gen `plog'=log(`p')
  gen `pclog'=log(`pc')
  gen `bnl'=0
  gen `en'=`n'
  gen `oldg'=lngamma(`en' +1.0)
  gen `em'=-1
  gen `e'=-1
  gen `ran1'=uniform()
  gen `ran2'=uniform()
  gen `j'=0
  gen `t'=1.0
  gen `g'=exp(-`am')
  gen `ds'=1
  gen `ts'=1
  gen `y'=-1
  gen `s'=sqrt(2.0*`am'*`pc')
  egen `sum1'=sum(`ds')
  noi di in gr "( Generating " _c 
  while `sum1' > 0 {
* method 1 for small denominators
  replace `bnl'=`bnl'+1 if ((`ran1' < `p') & (`j'<`n')& (`m'==1) )
  replace `j'=`j'+1     if ((`j'<`n') & (`m'==1))
  replace `ds'=0        if ((`j'==`n') & (`m'==1))
* method 2 for small expected numbers
  replace `em'=`em'+1    if ((`ds'==1) & (`m'==2))
  replace `t'=`t'*`ran1' if ((`ds'==1) & (`m'==2))
  replace `ds'=0         if ((`g'>`t') & (`m'==2))
* method 3 for general case
  replace `y'=sin(_pi*`ran1')/cos(_pi*`ran1') if ((`ds'==1) & (`m'==3))
  replace `em'=`s'*`y' +`am'                  if ((`ds'==1) & (`m'==3))
  replace `ts'=0 if (((0 >`em') | (`em' >=(`en'+1.0))) & (`ds'==1) & (`m'==3))
  #delimit ;
  replace `e'=1.2*`s'*(1.0+(`y'*`y'))*exp(`oldg'-lngamma(`em'+1.0)
    -lngamma(`en'-`em'+1.0)+(`em'*`plog')+(`en'-`em')*`pclog') if
    ((`ds'==1)&(`ts'==1) & (`m'==3));
  #delimit cr
  replace `ds'=0 if ((`ran2'<`e') & (`ds'==1) & (`ts'==1) & (`m'==3))
* get ready for the next time through the loop
  noi di in gr "." _c 
  replace `ran1'=uniform()
  replace `ran2'=uniform()
  replace `e'=-1
  replace `ts'=1
  drop `sum1'
  egen `sum1'=sum(`ds')
  }
* clean up for exit
  noi di in gr " )" 
  replace `bnl'=`em' if (`m'==2)
  replace `bnl'=`n' if ((`m'==2) & (`em'>`n'))
  replace `bnl'=int(`em'+0.5) if (`m'==3)
  replace `bnl'=`n'-`bnl' if `p' ~= `pp'
  gen bnlx=`bnl'
  noi di in gr "Variable " in ye "bnlx " in gr "created."
  lab var bnlx "Constructed binomial random variable"
  set type float
}
end