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.