*!version 1.1 1999 Joseph Hilbe
* version 1.0.0  1993 Joseph Hilbe, Walter Linde-Zwirble        (sg44: STB-28)
* Inverse Gaussian random number generator
* Example rndivg 1000 10 0.5; set obs 1000; mean=10; sigma=0.5; 
*         variance is sigma^2*mu^3

program define rndivg
   version 4.0
   set type double
   cap drop xig
   qui  {
      local cases `1'
      set obs `cases'
      mac shift
      local mu `1'
      mac shift
      local sig `1'
      local s=`sig'*`sig'*`mu'
      local peak = `mu'*((sqrt(4+9*`s'*`s')-(3*`s'))/2)
      local sq=sqrt(`s'*`mu'*`mu')
      #delimit ;
      local maxx=exp(-(`peak'-`mu')*(`peak'-`mu')/(2*`peak'*`mu'*`s')) 
           / sqrt(2*`s'*_pi*(`peak'^3)/`mu');
      #delimit cr
      tempvar ran1 ran2 ds ts sum1 t em y
      gen `ran1'=uniform()
      gen `ran2'=uniform()
      gen byte `ds'=1
      gen byte `ts'=1
      gen `em'=-1
      gen `t'=-1
      gen `y'=-1
      egen `sum1'=sum(`ds')
      noi di in gr "( Generating " _c
      while `sum1' > 0 {
         replace `y' = sin(_pi*`ran1')/cos(_pi*`ran1')
         replace `em'= `sq'*`y'+`peak' if `ds'==1
         replace `ts'=0 if ((0 > `em') & (`ds'==1))
         #delimit ;
         replace `t' = 0.66*(1+`y'*`y')*exp(-(`em'-`mu')*(`em'-`mu') /
            (2*`em'*`mu'*`s')) / sqrt(2*`s'*_pi*(`em'^3)/`mu') /`maxx' 
            if ((`ds'==1) & (`ts'==1));
         #delimit cr
         replace `ds'=0 if ((`ran2'<`t') & (`ds'==1) & (`ts'==1))
         replace `ran1'=uniform()
         replace `ran2'=uniform()
         replace `t' = -1
         replace `ts'=1
         drop `sum1'
         egen `sum1'=sum(`ds')
         noi di in gr "." _c
      }
      noi di in gr " )"
      gen xig = `em' 
      noi di in bl "Variable " in ye "xig " in bl "created."
      lab var xig "Inverse Gaussian random variable"
      set type float
   }
end