code: def inverse_cdf(su, W):
      
      j = 0
      s = W[0]
      M = su.shape[0]
      A = np.empty(M, 'int')
      for n in range(M):
          while su[n] > s:
              j += 1
              s += W[j]
          A[n] = j
      return A
docs: Inverse CDF algorithm for a finite distribution. 
  
          Parameters
          ----------
          su: (M,) ndarray
              M sorted uniform variates (i.e. M ordered points in [0,1]). 
          W: (N,) ndarray 
              a vector of N normalized weights (>=0 and sum to one)
  
          Returns
          ------- 
          A: (M,) ndarray
              a vector of M indices in range 0, ..., N-1