1
$\begingroup$

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[1]):
        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[0], N.shape[1]))
    for j in range(0, N.shape[1]):
        rank[:, j] = stats.rankdata(N[:, j], method='ordinal')
    return rank.astype(int) - 1

def reorder(rank, samples):
    rank_samples = np.zeros((samples.shape[0], samples.shape[1]))
    for j in range(0, samples.shape[1]):
        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[0]
y =dists.T[1]

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?

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

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.

$\endgroup$
0
$\begingroup$

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[1]):
    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.

$\endgroup$
  • 1
    $\begingroup$ 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). $\endgroup$ – vizie Nov 9 '17 at 23:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.