pearson_vectors.py

# Author: denis.engemann@gmail.com
# License: simplified BSD (3 clause)
# Note: code is based on scipy.stats.pearsonr

def ss(a, axis):
    return np.sum(a * a, axis=axis)

def compute_corr(x, y):
    x = np.asarray(x)
    y = np.asarray(y)
    mx = x.mean(axis=-1)
    my = y.mean(axis=-1)
    xm, ym = x - mx[..., None], y - my[..., None]
    r_num = np.add.reduce(xm * ym, axis=-1)
    r_den = np.sqrt(ss(xm, axis=-1) * ss(ym, axis=-1))
    r = r_num / r_den
    return r