'Simulate an inverse Gaussian distributed random variable
Public Function IG(ByVal mu As Double, ByVal lambda As Double) As Double
'Reference 1 : Devroye, Non-Uniform Random Variate Generation
'Reference 2 : Michael Schucany Haas

	Dim N, x, y, u, temp As Double

	N = Gauss()
	y = N * N
	x = mu + mu * mu * y / (2 * lambda) - mu * Sqr(4 * mu * lambda * y + mu * mu * y * y) / (2 * lambda)
	u = alea()
	If (u > mu / (mu + x)) Then
    		IG = mu * mu / x
	Else
    		IG = x
	End If
End Function