Testing the difference between two coefficients of correlation

gistfile1.py

"""Testing the difference between two coefficients of correlation

Following:

Thöni, H. (1977), Testing the Difference Between two Coefficients of Correlation.
Biometrical Journal, 19: 355–359. doi: 10.1002/bimj.4710190506
http://onlinelibrary.wiley.com/doi/10.1002/bimj.4710190506/abstract
"""

from math import log, sqrt
import numpy as np
from scipy import linalg
from scipy.stats import norm

n_samples = 10
n_trials = 2000
rho = 0.5
p_val = 0.05

def z_transform(r):
    return 0.5 * log((1.0 + r) / (1.0 - r))

def u(z1, z2):
    return (z1 - z2) / sqrt(1.0 / (n_samples - 3) + 1.0 / (n_samples - 3))

def corr_stat(r1, r2):
    return u(z_transform(r1), z_transform(r2))

def rvs(n_samples, rho):
    corr_mat = [[1.0, rho], [rho, 1.0]]  # correlation matrix
    L = linalg.cholesky(corr_mat)  # cholesky factorization
    X = np.dot(L, np.random.randn(2, n_samples))
    return X

n_reject = 0
for _ in range(n_trials):
    X1 = rvs(n_samples, rho)
    X2 = rvs(n_samples, rho)
    r1 = np.corrcoef(X1)[0, 1]
    r2 = np.corrcoef(X2)[0, 1]
    q = norm.cdf(corr_stat(r1, r2))  # quantile of normal distribution
    reject = (q <= (p_val / 2)) or (q >= (1.0 - p_val / 2))
    if reject:
        n_reject += 1

pct_reject = n_reject / float(n_trials)
print pct_reject  # should be equal to p_val