*!version 1.1 1999 Joseph Hilbe
* version 1.0.0  1993 Joseph Hilbe, Walter Linde-Zwirble        (sg44: STB-28)
* Inverse Gaussian random number generator with mu variable specified
* Example rndivgx mu, s(0.5); mean=mu; sigma=0.5;
*      mu=1/sqrt(eta); variance is sigma^2*mu^3

program define rndivgx
   version 3.1
   set type double
   local varlist "req ex"
   local options "`options' Sigma(real 1.0)"
   parse "`*'"
   parse "`varlist'", parse(" ")
   qui  {
      local mu `1'
      mac shift
      local s=`sigma'*`sigma'*`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 "Constructed inverse Gaussian random variable"
      set type float
  }
end