Scipyには、ピアソン相関係数を計算するための関数、scipy.stats.pearson
というものがあるのですが、残念ながらスパースな行列(scipy.sparse
)には対応していません。
実際、実装を見てみると(stats.py)、
mx = x.mean() my = y.mean() xm, ym = x - mx, y - my
という実装があり、単純に実装すると疎行列としては効率が悪くなります。
仮にxが疎行列だったとして、x - mx
は要素がほぼ -mx
な密行列になってしまうからです。
これは、定義の に当たりますが、式変形をすると、陽に平均を引く必要がなくなります。(Wikipediaより引用)
ということで、これをコードに直して、
import numpy as np import scipy.sparse def pearsons(a, b): if not (scipy.sparse.issparse(a) and scipy.sparse.issparse(b)): raise ValueError("only sparse arrays are supported") if a.shape != b.shape: raise ValueError("shape of sparse arrays must be same") if a.shape[0] != 1 or b.shape[0] != 1: raise ValueError("size of dimention 1 must be one") n = a.shape[1] a_sum = a.sum() b_sum = b.sum() nmr = n * a.multiply(b).sum() - a_sum * b_sum return (nmr / np.sqrt(n * a.multiply(a).sum() - (a_sum) ** 2) / np.sqrt(n * b.multiply(b).sum() - (b_sum) ** 2))
実際に検証してみると、以下のように同じ値になります。
import scipy.sparse from scipy.stats import pearsonr a = scipy.sparse.csr_matrix([1, 2, 1, 0, 1]) b = scipy.sparse.csr_matrix([1, 0, 1, 4, 2]) print pearsons(a, b) # -0.93250480824031368 print scipy.stats.pearsonr(a.toarray()[0], b.toarray()[0])[0] # -0.93250480824031379