Inverse Gaussian Generator

 

The fact that lamba*( -mu)2/(mu2 * x) is distributed as a chi-squared we have:

 

real *8 function wald(iseed,u,rlambda)

implicit real*8 (a-h,o-z)

chi2 = gauss(iseed)**2.

wald = u*(2.*rlambda+u*chi2-

        sqrt(4.*rlambda*u*chi2+(u*chi2)**2.))/(2.*rlambda)

xx = rand(int(iseed*10./9.))

if(xx.gt.u/(u+wald))wald=(u**2.)/wald

return

end

 

Where rand is a uniform generator, and gauss is a normal generator.