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

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