Qiita なんかで,数週間おきに,直線回帰についての記事が上がったりする。
いずれも sklearn だっけ?を使うか,あるいはそれこそ各変数の平均値を求める段階から自分でプログラムを書くかという両極端。後者は,そこまでやらなくてもという,前者はたかが直線回帰でそんなにたくさん記述が必要かという。いずれも,残念感満載の記事。それが何回も何回も出てくる。もう,うんざりだ。
Python に用意されているほどほどの機能を使って,端的に直線回帰を行う関数を書いてみる。
scipy.stats の linregress や pearsonr は使わないことにして(だって,それ使ったらむっちゃ簡単ですやん。Qiita 成り立たへん)
import numpy as np
def lm(x, y):
vx, vxy, _, vy = np.ravel(np.cov(x, y))
b = vxy / vx
return np.mean(y) - b * np.mean(x), b, vxy / np.sqrt(vx * vy)
独立変数のベクトル x と,従属変数のベクトル y を与えれば,回帰直線の切片と傾きと,ついでに変数間の相関係数を返す。必要最小限の機能だ。
引数は,numpy ndarray でも,リストでも,pandas データフレームから取り出した pandas.Series でもよい。
>>> def lm(x, y):
... vx, vxy, _, vy = np.ravel(np.cov(x, y))
... b = vxy / vx
... return np.mean(y) - b * np.mean(x), b, vxy / np.sqrt(vx * vy)
>>> lm(range(5), [3, 2, 6, 10, 15])
(0.7999999999999998, 3.2, 0.9444501377615284)
>>> import pandas as pd
>>> df = pd.DataFrame({"x": range(5), "y": [3, 2, 6, 10, 15]})
>>> lm(df.x, df.y)
(0.7999999999999998, 3.2, 0.9444501377615284)
一番,オーソドックスには
>>> x = np.array(range(5))
>>> y = np.array([3, 2, 6, 10, 15])
>>> intercept, slope, r = lm(x, y)
(0.7999999999999998, 3.2, 0.9444501377615284)
重相関係数が欲しいなら
>>> r**2
0.89198606271777
予測値が欲しいなら
>>> intercept + slope * x
array([ 0.8, 4. , 7.2, 10.4, 13.6])
散布図と回帰直線が欲しいなら,
import matplotlib.pyplot as plt
plt.scatter(x, y)
x2 = np.array([min(x), max(x)])
y2 = intercept + slope * x2
plt.plot(x2, y2)
plt.show()
他に,何が欲しい?