プログラミングの歴史--1

1969 年頃の話，一番最初に覚えたのは FORTRAN

入力装置 5 番から（古！！）データの個数とデータを読み，合計値と平均値を計算して出力装置 6 番に出力するプログラムは，以下のようなものだった。これでも，おったまげちょんなものだったのだ。

テストデータの入ったファイル（桁ズレしていたら目も当てられないことになった）

10
    1
    2
    3
    4
    5
    6
    7
    8
    9
   10

全てが大文字なのも，当時を忍ばせる（IBM のカードパンチャーは大文字しか打てなかった）

      INTEGER N
      REAL*8 X(1000)
      REAL*8 S, M
      READ(5, 100) N
 100  FORMAT(I2)
      DO 10 I = 1, N
        READ(5, 200) X(I)
 10   CONTINUE
 200  FORMAT(F5.0)
      S = 0
      DO 20 I = 1, N
        S = S + X(I)
 20   CONTINUE
      WRITE(6, 300) S
 300  FORMAT(" TOTAL =", G15.6)
      M = S / N
      WRITE(6, 400) M
 400  FORMAT("  MEAN =", G15.6)
      STOP
      END

今だと，これをコンパイルして実行するなら

$ gfortran prog.f
$ ./a.out < test.dat
 TOTAL =    55.0000    
  MEAN =    5.50000    

 

もしもし亀よ亀さんよ.... Python の pd.corr() が，どうしてそんなに遅いのか？

「みんな大好き，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

IPA フォント(Python)

IPA フォントをインストールしたのにちゃんと使えないということだったが，

~/.matplotlib の中にある，2 つのファイルを削除する（いきなり削除するのが怖い人は，名前を変える）

-rw-r--r--  1 foo  503  104286 10 31  2017 fontList.cache0  これと
-rw-r--r--  1 foo  503  165099  9 12 15:58 fontList.json0　　これ
-rw-r--r--  1 foo  503  136304  1 11 09:28 fontlist-v300.json（これは例え削除しても再度作られる）
-rw-r--r--@ 1 foo  503   33307  1 11 09:25 matplotlibrc
drwxr-xr-x  2 foo  503      64  1 23  2014 tex.cache/

みんな大好き，R と Python の速度比較

統計計算で，R と Python のどちらが速いかというのは多くの人の関心事のようで，ネット上でもいくつかある。

PythonとR言語のプログラム処理速度を比較！
https://ai-trend.jp/programming/python/speed-comp-python-r/
というのもそれ類のページで，初見 　2018/02/01 で 2020/04/14 に更新されたようである。

それによると，どうも R の分が悪いようなので，掲載例以外の統計処理についても再試してみる。

システム構成
 mxcOS Catalina バージョン 10.15.7
 MacBook Pro (Retina, Mid 2012)
 プロセッサ　2.7 GHz クアッドコアIntel Core i7
 メモリ　16 GB 1600 MHz DDR3
R version 4.0.3
Python 3.8.1

まず，データフレームを用意する。100000行，1500列 のデータフレーム。各列は平均値50，標準偏差10の正規乱数を小数点以下2桁までで丸めた値を使う。変数名は X1, X2, ..., X1500。

CSV ファイル "test.csv" へ書き出す。

ちなみに，書き出しに要する時間は system.time() を用いて以下のように測定した。
R: system.time({write.csv(x, file="test.csv", row.names=FALSE)})  135.968 sec.
135.968 秒かかった。

データフレームへの読み込み

R:       system.time({df = read.csv("test.csv")})  78.787 sec.

列ごとの変数型を指定すると速くなる
R:      system.time({df = read.csv("test.csv", colClasses="numeric")})  33.607 sec.

data.table パッケージの fread() はもっともっと速い。Python の 5.7 倍速
R:      library(data.table)
        system.time(df = fread("test.csv"))  5.023 sec.

Python: start = time(); df = pd.read_csv('test.csv'); print(time() - start)  28.793 sec.

各変数ごとに単変量解析　　R の圧勝。Python の 5〜20 倍速

和　　R の圧勝

R:      s = colSums(df)  0.739 sec.

Python: s = df.sum()     5.416 sec.

平均値　　R の圧勝
R:      m = colMeans(df) 0.662 sec.

Python: m = df.mean()    5.352 sec.

不偏分散　　R の圧勝
R:      v = sapply(df, var) 0.552 sec.

Python: v = df.var(ddof=1)  8.466 sec.

標準偏差　　R の圧勝
R:  sd = sapply(df, sd)       0.561 sec.

Python: sd = df.std(ddof=1)  11.574 sec.

相関係数（単相関係数; 二変数相関）
R:      cor(df[,1], df[,2])                   0.011 sec.

Python: import numpy as np
        np.corrcoef(df['X1'], df['X2'])[0,1]  0.011 sec.

Python: from scipy.stats import pearsonr
        pearsonr(df['X1'], df['X2'])[0]       0.004 sec.

相関係数行列  R の圧勝　ただし，ケンドールの順位相関係数は話にならないほど遅い

以下を追記した

R ユーザ，歓喜の涙 -- ケンドールの順位相関係数の計算が Python の倍速！
https://blog.goo.ne.jp/r-de-r/e/2716d7ca3c9fb3a1eb643d83d10dfecd

以下で，予測式中の p は変数の個数（df の列数）

ピアソンの積率相関係数　R が 8.7 倍速い

R:      r = cor(df)    100.509 sec.　予測式 sec = 0.00001828 * p ** 2

Pythom: r = df.corr()  360.236 sec.　予測式 sec =  0.0001586 * p ** 2

スピアマンの順位相関係数　R が 20 倍速い

R:      rs = cor(df, method="spearman")  140.892sec.　予測式 sec = 0.00006395 * p ** 2

Python: rs = df.corr(method='spearman') 4238.496 sec.　予測式 sec = 0.001282 * p ** 2

ケンドールの順位相関係数　R は測定不能なほど遅い

R:      rk = cor(df, method="kendall")  推定値 106216279 sec. 約 3.36 年 !

Python: rk = df.corr(method='kendall')  28072.694 sec.　予測式 sec = 0.01353 * p ** 2

分散・共分散行列（R, Python とも不偏）  Python の圧勝 20倍速い

R:      v = var(df)  106.829 sec.

Python: v = df.cov()   4.987 sec.

重回帰分析　　Python の圧勝

数個の独立変数を使うとき

R:      ans = lm(X1 ~ X2+X3+X4+X5+X6+X7+X8+X9+X10, data=df)
        summary(ans)
        ans$coeff
        str(summary(ans))
        summary(ans)$r.squared  # 0.163 sec.

Python: import scipy as sp
        from sklearn.linear_model import LinearRegression
        clf = LinearRegression()
        y = df['X1']
        x = df[['X2', 'X3', 'X4', 'X5', 'X6', 'X7', 'X8', 'X9', 'X10']]
        ans = clf.fit(x, y)
        ans.coef_
        ans.intercept_
        ans.score(x, y)
        print(time() - start)  # 0.031 sec.

たくさんの独立変数を使うとき

R:      x = as.matrix(df[-1])
        ans = lm(df[, 1] ~ x)
        summary.ans = summary(ans)
        coeff = ans$coeff
        R2 = summary(ans)$r.squared  # 412.829 sec.

Python: y = df['X1']
        x = df.drop('X1', axis=1)
        ans = clf.fit(x, y)
        coef = ans.coef_
        intercept = ans.intercept_
        R2 = ans.score(x, y)        # 136.559 sec.