「みんな大好き，R と Python の速度比較 」で書いたんだけど，
「Python での積率相関係数行列（ピアソンの積率相関係数）を求める df.corr() が R に比べて 9 倍ほども遅い！」ということであったが，
一方で，「R の分散共分散行列は R より20 倍も速い」ということ。
100000x1500 のデータフレームを使って計算してみる。
df.corr() を使うと，366.456 sec. かかってしまう。
import stat2_lib
import numpy as np
from time import time
import pandas as pd
import os
start = time(); df = pd.read_csv("test.csv"); print('***** read', time() - start)
nr, nc = df.shape
start = time()
corr = df.corr()
print('\n***** df.corr()', time() - start)
print(np.array(corr.iloc[:3, :3]))
[[ 1.00000000e+00 -1.77002905e-04 -2.76211732e-03]
[-1.77002905e-04 1.00000000e+00 1.54937687e-03]
[-2.76211732e-03 1.54937687e-03 1.00000000e+00]]
df.corr() を使うと，366.456 sec. かかってしまう。
なんでこうなるのか。
分散共分散行列を S2xy，変数 x, y の標準偏差を SDx, SDy とすれば，相関係数行列 r は Sxy / (SDx * SDy) なのだから，
start = time()
cov = df.cov()
std = np.sqrt(np.diag(cov))
r = np.array(cov)
nr, nc = r.shape
for i in range(nc):
for j in range(nc):
r[i, j] = r[i, j] / std[i] / std[j]
print('\n***** df.cov() to r', time() - start)
print(cov.iloc[:3, :3])
print(std[:3])
print(r[:3, :3])
これによれば，結果は 6.840 sec. で得られる。
X1 X2 X3 # 分散・共分散行列
X1 99.530159 -0.017654 -0.275388
X2 -0.017654 99.949528 0.154801
X3 -0.275388 0.154801 99.873711
[9.97648027 9.99747606 9.99368355] # 標準偏差
[[ 1.00000000e+00 -1.77002905e-04 -2.76211732e-03] # 相関係数行列
[-1.77002905e-04 1.00000000e+00 1.54937687e-03]
[-2.76211732e-03 1.54937687e-03 1.00000000e+00]]
なお，numpy には corrcoef() というものがあって，これでも同じく，あっという間に 4.021 sec. 解が得られる。
corrcoef_r = np.corrcoef(df.T)
print(corrcoef_r[:3, :3])
[[ 1.00000000e+00 -1.77002905e-04 -2.76211732e-03]
[-1.77002905e-04 1.00000000e+00 1.54937687e-03]
[-2.76211732e-03 1.54937687e-03 1.00000000e+00]]
なおなお，自前の fortran プログラムでさえも 88.453 sec. で答えを出せるよ。
[[ 1.00000000e+00 -1.77002905e-04 -2.76211732e-03]
[-1.77002905e-04 1.00000000e+00 1.54937687e-03]
[-2.76211732e-03 1.54937687e-03 1.00000000e+00]]
subroutine cor(nr, nc , x, r, method)
implicit none
integer, intent(in) :: nr, nc
real(8), dimension (:, :), intent(inout) :: x
real(8), dimension (nc, nc), intent(out) :: r
character(*), optional, intent(in) :: method
real(8), dimension (nc) :: mean, sd, w
real(8) fn, fn1, factor, vx
real(8), dimension (nr) :: data
integer i, j, k
if (method == 'spearman') then
do i = 1, nc
data = x(:, i)
call rank(nr, data)
x(:, i) = data
end do
end if
do i = 1, nc
mean(i) = 0
do j = 1, nc
r(i, j) = 0
end do
end do
do i = 1, nr
do j = 1, nc
w(j) = x(i, j)
fn = 1 / dble(i)
factor = (i - 1) * fn
w(j) = w(j) - mean(j)
vx = w(j)
mean(j) = mean(j) + vx * fn
vx = vx * factor
do k = 1, j
r(k, j) = r(k, j) + vx * w(k)
end do
end do
end do
fn = nr
fn1 = nr - 1
do i = 1, nc
do j = 1, i
r(j, i) = r(j, i) / fn1
end do
sd(i) = sqrt(r(i, i))
end do
do i = 1, nc
do j = 1, i
if (i == j) then
r(i, i) = 1
else
r(j, i) = r(j, i) / (sd(j) * sd(i))
r(i, j) = r(j, i)
end if
end do
end do
end subroutine cor
