Pythonによる単回帰分析入門（数値計算/数式処理系の6ライブラリの比較）

# 概要

Pythonで**単回帰分析をしたいとき**、つまり、**回帰直線を求めたいとき**、それを簡単に叶えてくれるライブラリが複数用意されています。ここでは「numpy」「scipy」「scikit-learn」「statsmodels」「sympy」「tensorflow」の6つライブラリを使って回帰直線（$y=ax+b$ のパラメータ $a$ と $b$）を求めるための方法／コードについて比較・整理をしました。

また、おまけとしてデータ可視化ライブラリ「seaborn」で自動描画される回帰直線から強引に回帰式パラメータを求める方法についても記載しました。

コードの動作確認は、GoogleColab.（Python 3.6.9）により行なっています。

# 準備

各種ライブラリをインポートをします。そして、一次関数 $8x+10$ に対して、正規分布 $N(0,6^2)$ に従う乱数 $\varepsilon$ を加えた $y=8x+10+\varepsilon$ からサンプルデータ（$8$個のセット）を生成します。

```py:データの準備
import numpy as np
import pandas as pd
import scipy
import sklearn.linear_model
import statsmodels
import sympy
import tensorflow as tf
import matplotlib.pyplot as plt
import seaborn as sns

np.random.seed(10) # ランダムシートの固定
n = 8
x_range = (10,60)
x = np.round(np.linspace(*x_range,n),2)
e = np.random.normal(0,6,n)
y = np.round(0.8 * x + 10 + e,2)
df = pd.DataFrame({'x':x,'y':y})
display(df)
```

実行結果は、次のようになります。<br>![2020-12-29_21h15_28.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/152805/dd6b364f-ade0-a17c-13ce-c05ceccdca4d.png)

このデータを、matplotlib でプロットしていきます。

```py:データのプロット
fig,ax = plt.subplots(1,1)
ax.scatter(x,y)
ax.set_ylim(0,80)       # Y軸の描画範囲
ax.grid()               # グリッドをオン
ax.set_axisbelow(True)  # グリッドを背面に移動
plt.show()              # 出力セルに描画
```

実行結果は、次のようになります。<br>![2020-12-29_21h18_29.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/152805/1910107c-7995-03ac-d456-ee5cae28bb1a.png)

上記にプロットしたデータ（入力 $x$、出力 $y$）に対して、回帰式 $y=ax+b$  のパラメータ $a$ を $b$ を求めていきます。

# 回帰分析：statsmodels 編

`statsmodels.formula.api.ols(...)` を利用して回帰式を求めます。

- 公式リファレンス：https://www.statsmodels.org/stable/generated/statsmodels.formula.api.ols.html

```py:statsmodels
print(f'> statsmodels {statsmodels.__version__}')
c = statsmodels.formula.api.ols(formula='y ~ x',data=df).fit()
print(f'a = {c.params["x"]:.2f}')
print(f'b = {c.params["Intercept"]:.2f}')
```

```:実行結果
> statsmodels 0.10.2
a = 0.72
b = 13.27
```

また、上記につづけて、`print(c.summary())` とすれば、**決定係数**や**パラメータの信頼区間**といった情報を確認することができます。

```:c.summary()
                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.856
Model:                            OLS   Adj. R-squared:                  0.832
Method:                 Least Squares   F-statistic:                     35.69
Date:                Tue, 29 Dec 2020   Prob (F-statistic):           0.000987
Time:                        12:32:03   Log-Likelihood:                -23.986
No. Observations:                   8   AIC:                             51.97
Df Residuals:                       6   BIC:                             52.13
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     13.2717      4.676      2.838      0.030       1.830      24.713
x              0.7230      0.121      5.974      0.001       0.427       1.019
==============================================================================
Omnibus:                        5.195   Durbin-Watson:                   2.078
Prob(Omnibus):                  0.074   Jarque-Bera (JB):                1.759
Skew:                          -1.142   Prob(JB):                        0.415
Kurtosis:                       3.242   Cond. No.                         91.3
==============================================================================
```

# 回帰分析：scikit-learn 編

`sklearn.linear_model.LinearRegression()` を利用して回帰式を求めます。単回帰の場合であっても、入力は行列（ndim=2）にして与える必要があります。今回の場合であれば、`X.shape` が `(8,)` ではダメで、`(8,1)` に変形しないといけません。

- 公式リファレンス：https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html

 
```py:scikit-learn
print(f'> sklearn {sklearn.__version__}')
X = df['x'].values.reshape(len(df),1) # X.shape -> (8,1)
y = df['y'].values                    # y.shape -> (8,)
lr = sklearn.linear_model.LinearRegression()
lr.fit(X,y)
print(f'a = {lr.coef_[0]:.2f}')
print(f'b = {lr.intercept_:.2f}')
```

```:実行結果
> sklearn 0.22.2.post1
a = 0.72
b = 13.27
```

# 回帰分析：numpy 編

`numpy.polyfit(...)` を利用して回帰式を求めます。

- 公式リファレンス：https://numpy.org/doc/stable/reference/generated/numpy.polyfit.html

```py:numpy
print(f'> numpy {np.__version__}')
x = df['x'].values
y = df['y'].values
c = np.polyfit(x,y,1)
print(f'a = {c[0]:.2f}')
print(f'b = {c[1]:.2f}')
```

```:実行結果
> numpy 1.19.4
a = 0.72
b = 13.27
```

# 回帰分析：scipy 編

`scipy.optimize.leastsq(...)` あるいは `scipy.optimize.curve_fit(...)` を利用して回帰式を求めます。

- 公式リファレンス：https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.html

- 公式リファレンス：https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html

```py:scipy.optimize.leastsq(...) 
print(f'> scipy {scipy.__version__}')
x = df['x'].values
y = df['y'].values
c = scipy.optimize.leastsq(lambda c,x,y:y-(c[0]*x+c[1]),(0,0),args=(x,y))
print(f'a = {c[0][0]:.2f}')
print(f'b = {c[0][1]:.2f}')
```

```py:scipy.optimize.curve_fit(...) 
print(f'> scipy {scipy.__version__}')
x = df['x'].values
y = df['y'].values
c = scipy.optimize.curve_fit(lambda x,a,b:a*x+b,x,y)
print(f'a = {c[0][0]:.2f}')
print(f'b = {c[0][1]:.2f}')
```

実行結果は、いずれも同じで次のようになります。

```:実行結果
> scipy 1.4.1
a = 0.72
b = 13.27
```

# 回帰分析：sympy 編

sympyは数値計算ではなく、数式処理のライブラリになります。よって、最小二乗法に基づく解析解を求めることができます。

コードのなかの `J` は、次のように残差の平方和（最小二乗法における最小化対象）になります。

$$ J(a,b) = \sum_{i=1}^{n}\left( y - ( ax+b ) \right)^2$$

また `dJa` と `dJb` は、それぞれ $J(a,b)$ の偏微分の解析解です。また、`eq1` と `eq2` は、次のように、それらをイコールゼロとした方程式になります。

$$ \frac{\partial J(a,b)}{\partial a}=0,\ \ \  \frac{\partial J(a,b)}{\partial b}=0 $$

そして、`sympy.solve((eq1,eq2),(a,b))` で、これらの連立方程式を解いています。

```py:sympy
print(f'> sympy {sympy.__version__}')
rx = df['x'].values
ry = df['y'].values

a, b = sympy.symbols('a,b')
n = sympy.symbols('n', integer=True, positive=True)
i = sympy.Idx('i')
x = sympy.IndexedBase('x')
y = sympy.IndexedBase('y')

xi = [ (x[i+1],rx[i]) for i in range(len(rx)) ]
yi = [ (y[i+1],ry[i]) for i in range(len(ry)) ]

J = sympy.Sum((y[i]-(a*x[i]+b))**2,(i,1,n)) # 残差平方和

dJa = sympy.diff(J,a) # aで偏微分
dJb = sympy.diff(J,b) # bで偏微分

dJa = dJa.subs(n,len(df)).doit().subs(xi).subs(yi) # xi,yiの代入
dJb = dJb.subs(n,len(df)).doit().subs(xi).subs(yi) # xi,yiの代入

eq1, eq2 = sympy.Eq(dJa,0), sympy.Eq(dJb,0) # 連立方程式を解く

r = sympy.solve((eq1,eq2),(a,b))
print(f'a = {float(r[a].subs(xi,yi)):.2f}')
print(f'b = {float(r[b].subs(xi,yi)):.2f}')
```

```:実行結果
> sympy 1.1.1
a = 0.72
b = 13.27
```

# 回帰分析：tensorflow 編

TensorFlow（2.X）の低レイヤの関数 `tf.GradientTape(...)` を使った勾配降下法により、回帰分析を行ないます。いい感じの結果を得るためには、学習率や反復回数などの調整が必要です。

```py:tensorflow
print(f'> tensorflow {tf.__version__}')

x = df['x'].values
y = df['y'].values

a = tf.Variable(0.)
b = tf.Variable(0.)

for i in range(50000):
  with tf.GradientTape() as tape:
    y_pred = a*x+b
    loss = tf.math.reduce_mean(tf.square(y_pred-y))
  gradients = tape.gradient(loss,[a,b])
  a.assign_sub(0.0005 * gradients[0]) # 学習率 0.0005
  b.assign_sub(0.0005 * gradients[1])

  if (i + 1) % 10000 == 0:
    print(f'Step:{i+1} a = {a.numpy():.2f} b = {b.numpy():.2f}')

print(f'a = {a.numpy():.2f}')
print(f'b = {b.numpy():.2f}')
```

```:実行結果
> tensorflow 2.4.0
Step:10000 a = 0.77 b = 11.07
Step:20000 a = 0.73 b = 12.90
Step:30000 a = 0.72 b = 13.21
Step:40000 a = 0.72 b = 13.26
Step:50000 a = 0.72 b = 13.27
a = 0.72
b = 13.27
```

# 回帰分析：seaborn 編（おまけ）

可視化ライブラリ `seaborn.regplot(...)` では、次のように回帰直線と$95\%$信頼区間をデフォルトで描画してくれます。

![2020-12-29_22h50_31.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/152805/b47aae06-ec70-868e-c1d8-2ff159319102.png)

しかし残念ながら、それら回帰直線や信頼区間のパラメータを参照する手段は用意されていないようなのです（[参考](https://stackoverflow.com/questions/33490833/display-regression-equation-in-seaborn-regplot)）。そこで、回帰直線のプロットデータを取得し、そこから式のパラメータを計算することができます。

```py:seaborn
print(f'> seaborn {sns.__version__}')
fig,ax = plt.subplots(1,1)
d = sns.regplot(x='x', y='y', data=df, ax=ax) # プロットデータの取得
ax.set_ylim(0,80)
ax.grid() 
ax.set_axisbelow(True)
plt.show()

xr = d.get_lines()[0].get_xdata() # 図内の回帰直線のXデータ
yr = d.get_lines()[0].get_ydata() # 図内の回帰直線のYデータ
a = (yr[-1]-yr[0])/(xr[-1]-xr[0])
b = yr[0]-a*xr[0]
print(f'a = {a:.2f}')
print(f'b = {b:.2f}')
```

実行結果は、次のようになります。

```:実行結果
> seaborn 0.11.0
a = 0.72
b = 13.27
```

# 参考

- 「Pythonデータ分析／機械学習のための基本コーティング pandasライブラリ活用入門」@インプレス
- 「TensorFlow2 プログラミング実装ハンドブック」@秀和システム
- 「[最小二乗法でカーブフィッティング。関数３つを使い比べ](http://ailaby.com/least_square/)」@[コード7区](http://ailaby.com/)