# Iman-Conover implementation for correlated randoms in python with numpy/scipy

I've been trying to implement the iman-conover method in python so I could generate correlated random numbers from distributions other than normal (I use a normal & uniform in my example below).

I started by following this blog and realized that I preferred the explanation/walk through at Howard Rudd's site even though it requires some translation from VBA.

There are a few functions (modified from the first blog):

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

def ic_m(n, d):
a = np.arange(1, (n+1))
p = stats.norm.ppf(a/(n+1))
p = normalize(p)
score = np.zeros((n, d))
for j in range(0, score.shape):
score[:, j] = np.random.permutation(p)
return score

def normalize(v):
norm=np.linalg.norm(v, ord=1)
return v/norm

def rank(N):
rank = np.zeros((N.shape, N.shape))
for j in range(0, N.shape):
rank[:, j] = stats.rankdata(N[:, j], method='ordinal')
return rank.astype(int) - 1

def reorder(rank, samples):
rank_samples = np.zeros((samples.shape, samples.shape))
for j in range(0, samples.shape):
s = np.sort(samples[:, j])
rank_samples[:, j] = s[rank[:,j]]
return rank_samples


And then the actual test of the method:

n, d = 1000, 2
corrTar = .2
S = np.array(([1., corrTar],
[corrTar, 1.]))
C = np.linalg.cholesky(S)
M = ic_m(n,d)
D = (1./n) * np.dot(M.T, M)
E = np.linalg.cholesky(D)
N = np.dot(np.dot(M, np.linalg.inv(E)), C)
R = rank(N)

dists = np.array((
stats.norm.ppf(np.random.uniform(0.0, 1.0, n), loc=0, scale=1),
stats.uniform.ppf(np.random.uniform(0.0, 1.0, n), 0, 1)
))
dists = reorder(R, dists.T)
np.corrcoef(dists.T)

x =dists.T
y =dists.T


This generates roughly what I'd expect when creating a scatter plot of the data (I'd post images but am limited to 2 links). If I iterate through that process 1,000 times and record the spearman correlation after each test I get a distribution that is centered around the desired correlation.

The problem occurs when I increase the desired correlation to anything above 0.6 (0.85 as an example). The resulting correlation is centered around 0.60 and never even approaches the 0.85. I do not believe this is normal behavior and haven't been able to see where I misrepresented the method in my code.

Can anyone see what I cannot?

• I'm getting nothing but nans for M... can you tell me what version of Python you're running? – JD Long Nov 9 '17 at 21:06
• @JDLong I believe I was on 3.5 at this point. I’ll dig up this project tonight to check. – vizie Nov 9 '17 at 21:17
• 3.5.4 |Anaconda, Inc.| (default, Oct 13 2017, 11:22:58) [GCC 7.2.0] numpy 1.13.3 scipy 0.19.1 – vizie Nov 9 '17 at 23:34

The issue seems to be with np.linalg.cholesky (or how its being used here and in the original post). If I run your code as posted, but substitute scipy.linalg.cholesky for np.linalg.cholesky, the results are more reasonable when corrTar >= 0.6.

Note that the numpy cholesky returns a lower triangular matrix and the scipy cholesky returns an upper triangular matrix. Transposing the numpy cholesky matrices similarly resolves the issue.

looks like you changed up ic_m() from the author's original which is as follows:

def ic_m(n, d):
p = [norm.ppf(1./(x+1)) for x in range(1, n+1)]
p = (p - np.mean(p))*(1. / np.std(p))
score = np.zeros((n, d))
for j in range(0, score.shape):
score[:, j] = np.random.permutation(p)
return score


what was driving the change? I can get the code above to seemingly work, but not yours.

• I dug into this more at the time and am having trouble remembering, but here's what I noticed with the 2 different functions when I ran them just now: 2 scatter plots. Both exhibit the same problem above where any correlation coefficient target greater than 0.6 never comes close (both those were with a target of .9). – vizie Nov 9 '17 at 23:45