docs: def rand_imancon(X, rho):
      """Iman-Conover Method to generate random ordinal variables
          (Implementation adopted from Ekstrom, 2005)
  
      x : ndarray
          <obs x cols> matrix with "cols" ordinal variables
          that are uncorrelated.
  
      rho : ndarray
          Spearman Rank Correlation Matrix
  
      Links
      * Iman, R.L., Conover, W.J., 1982. A distribution-free approach to
          inducing rank correlation among input variables. Communications
          in Statistics - Simulation and Computation 11, 311–334.
          https://doi.org/10.1080/03610918208812265
      * Ekstrom, P.-A., n.d. A Simulation Toolbox for Sensitivity Analysis 57.
          http://ecolego.facilia.se/ecolego/files/Eikos_thesis.pdf
      """
      import numpy as np
      import scipy.stats as sstat
      import scipy.special as sspec
  
      # data prep
      n, d = X.shape
  
      # Ordering
      T = np.corrcoef(X, rowvar=0)  # T
      Q = np.linalg.cholesky(T)  # Q
      invQ = np.linalg.inv(Q)  # inv(Q)
      P = np.linalg.cholesky(rho)  # P
      S = np.dot(invQ, P)  # S=P*inv(Q)
  
      # get ranks of 'X'
      rnks = np.nan * np.empty((n, d))
      for k in range(0, d):
          rnks[:, k] = sstat.rankdata(X[:, k], method='average')
      # create Rank Scores
      rnkscore = -np.sqrt(2.0) * sspec.erfcinv(2.0 * rnks / (n + 1))
  
      # the 'Y' variables have the same correlation matrix
      Y = np.dot(rnkscore, S.T)
  
      # get ranks of 'Y'
      rnks = np.nan * np.empty((n, d))
      for k in range(0, d):
          rnks[:, k] = sstat.rankdata(Y[:, k], method='average')
      rnks = rnks.astype(int)
  
      # Sort X what will decorrelated X
      X = np.sort(X, axis=0)
      # Rerank X
      Z = np.nan * np.empty((n, d))
      for k in range(0, d):
          Z[:, k] = X[rnks[:, k] - 1, k]
  
      # done
      return Z