Inverse Gaussian Generator

The fact that lamba*( -mu)²/(mu² * 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.