code: #!/usr/bin/env python
  
  """
  Random variate generator for the generalized inverse Gaussian distribution.
  Reference: L Devroye. Random variate generation for the generalized inverse Gaussian distribution.
             Statistics and Computing, 24(2):239-246, 2014.
  
  """
  
  
  import math
  from scipy import random
  
  
  def psi(x, alpha, lam):
      f = -alpha*(math.cosh(x)-1)-lam*(math.exp(x)-x-1)
      return f
  
  
  def dpsi(x, alpha, lam):
      f = -alpha*math.sinh(x)-lam*(math.exp(x)-1)
      return f
  
  
  def g(x, sd, td, f1, f2):
      if (x >= -sd) and (x <= td):
          f = 1
      elif x > td:
          f = f1
      elif x < -sd:
          f = f2
  
      return f
  
  
  def gigrnd(p, a, b):
      # setup -- sample from the two-parameter version gig(lam,omega)
      p = float(p); a = float(a); b = float(b)
      lam = p
      omega = math.sqrt(a*b)
  
      if lam < 0:
          lam = -lam
          swap = True
      else:
          swap = False
  
      alpha = math.sqrt(math.pow(omega,2)+math.pow(lam,2))-lam
  
      # find t
      x = -psi(1, alpha, lam)
      if (x >= 1/2) and (x <= 2):
          t = 1
      elif x > 2:
          t = math.sqrt(2/(alpha+lam))
      elif x < 1/2:
          t = math.log(4/(alpha+2*lam))
  
      # find s
      x = -psi(-1, alpha, lam)
      if (x >= 1/2) and (x <= 2):
          s = 1
      elif x > 2:
          s = math.sqrt(4/(alpha*math.cosh(1)+lam))
      elif x < 1/2:
          if alpha == 0:
              s = 1/lam
          else:
              if lam==0:
                  s =  math.log(1+1/alpha+math.sqrt(1/math.pow(alpha,2)+2/alpha))
              else:
                  s = min(1/lam, math.log(1+1/alpha+math.sqrt(1/math.pow(alpha,2)+2/alpha)))
  
      # find auxiliary parameters
      eta = -psi(t, alpha, lam)
      zeta = -dpsi(t, alpha, lam)
      theta = -psi(-s, alpha, lam)
      xi = dpsi(-s, alpha, lam)
  
      p = 1/xi
      r = 1/zeta
  
      td = t-r*eta
      sd = s-p*theta
      q = td+sd
  
      # random variate generation
      while True:
          U = random.random()
          V = random.random()
          W = random.random()
          if U < q/(p+q+r):
              rnd = -sd+q*V
          elif U < (q+r)/(p+q+r):
              rnd = td-r*math.log(V)
          else:
              rnd = -sd+p*math.log(V)
  
          f1 = math.exp(-eta-zeta*(rnd-t))
          f2 = math.exp(-theta+xi*(rnd+s))
          if W*g(rnd, sd, td, f1, f2) <= math.exp(psi(rnd, alpha, lam)):
              break
  
      # transform back to the three-parameter version gig(p,a,b)
      rnd = math.exp(rnd)*(lam/omega+math.sqrt(1+math.pow(lam,2)/math.pow(omega,2)))
      if swap:
          rnd = 1/rnd
  
      rnd = rnd/math.sqrt(a/b)
      return rnd
  
  
  
apis: scipy.random.random