Pythonでの線形回帰（statmodels, scikit-learn, PyMC3）

# はじめに

Pythonではいくつか線形回帰をするために使えるライブラリがあります。個人的に線形回帰をする必要にせまられ、そのための方法を調べたのでメモを兼ねてシェアしたいと思います。使ったライブラリは以下：
- statmodels 
- scikit-learn
- PyMC3


## データ準備

まず適当なデータを用意します。今回は以下の式にノイズ $\epsilon$ を加えた人工データを使います。

$$y = \beta_0 + \beta_1 x_1 + \beta_2 x2 + \epsilon$$

ここで推定するのは、$\beta_0, \beta_1, \beta_2$の値。

```python
import numpy as np
import pandas as pd
import random
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Generate random data
beta = [1.2, 0.5]
# prep data
x1 = np.random.random(size=1000)*5 
x2 = np.random.random(size=1000)*10 
x = np.transpose([x1, x2])
y = np.dot(x, beta)  + np.random.normal(scale=0.5, size=1000)
#data = dict(x=x, y=y)
data = dict(x1=x1, x2=x2, y=y)
df = pd.DataFrame(data)
# visualisation
plt.scatter(x1, y, color='b')
plt.scatter(x2, y, color='orange')
# 3D
fig = plt.figure()
ax = Axes3D(fig)
ax.scatter(x1, x2, y)
plt.show()

```

![output_2_0.png](https://qiita-image-store.s3.amazonaws.com/0/154129/a94c6ac7-7ec0-59c4-f111-9b4a124a5cde.png)
![output_2_1.png](https://qiita-image-store.s3.amazonaws.com/0/154129/2deaa2f7-a70a-48e7-ac14-f3596fba87ef.png)



## Statmodelsを使う場合

```python
import statsmodels.api as sm

x = sm.add_constant(x)
results = sm.OLS(endog=y, exog=x).fit()
results.summary()
```


<table class="simpletable">
<caption>OLS Regression Results</caption>
<tr>
  <th>Dep. Variable:</th>            <td>y</td>        <th>  R-squared:         </th> <td>   0.951</td>
</tr>
<tr>
  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.951</td>
</tr>
<tr>
  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   9745.</td>
</tr>
<tr>
  <th>Date:</th>             <td>Fri, 10 Mar 2017</td> <th>  Prob (F-statistic):</th>  <td>  0.00</td> 
</tr>
<tr>
  <th>Time:</th>                 <td>09:58:59</td>     <th>  Log-Likelihood:    </th> <td> -724.14</td>
</tr>
<tr>
  <th>No. Observations:</th>      <td>  1000</td>      <th>  AIC:               </th> <td>   1454.</td>
</tr>
<tr>
  <th>Df Residuals:</th>          <td>   997</td>      <th>  BIC:               </th> <td>   1469.</td>
</tr>
<tr>
  <th>Df Model:</th>              <td>     2</td>      <th>                     </th>     <td> </td>   
</tr>
<tr>
  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>   
</tr>
</table>
<table class="simpletable">
<tr>
    <td></td>       <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th> <th>[95.0% Conf. Int.]</th> 
</tr>
<tr>
  <th>const</th> <td>    0.0499</td> <td>    0.042</td> <td>    1.181</td> <td> 0.238</td> <td>   -0.033     0.133</td>
</tr>
<tr>
  <th>x1</th>    <td>    1.1823</td> <td>    0.011</td> <td>  107.081</td> <td> 0.000</td> <td>    1.161     1.204</td>
</tr>
<tr>
  <th>x2</th>    <td>    0.4983</td> <td>    0.005</td> <td>   91.004</td> <td> 0.000</td> <td>    0.488     0.509</td>
</tr>
</table>
<table class="simpletable">
<tr>
  <th>Omnibus:</th>       <td> 0.654</td> <th>  Durbin-Watson:     </th> <td>   2.079</td>
</tr>
<tr>
  <th>Prob(Omnibus):</th> <td> 0.721</td> <th>  Jarque-Bera (JB):  </th> <td>   0.599</td>
</tr>
<tr>
  <th>Skew:</th>          <td>-0.059</td> <th>  Prob(JB):          </th> <td>   0.741</td>
</tr>
<tr>
  <th>Kurtosis:</th>      <td> 3.023</td> <th>  Cond. No.          </th> <td>    17.2</td>
</tr>
</table>



## scikit-learnを使う場合


```python
from sklearn import linear_model
# compare different regressions
lreg = linear_model.LinearRegression()
lreg.fit(x, y)
print("Linear regression: \t", lreg.coef_)
breg = linear_model.BayesianRidge()
breg.fit(x, y)
print("Bayesian regression: \t", breg.coef_)
ereg = linear_model.ElasticNetCV()
ereg.fit(x, y)
print("Elastic net:  \t\t", ereg.coef_)
print("true parameter values: \t", beta)
```

    Linear regression: 	 [ 1.18232244  0.49832431]
    Bayesian regression: 	 [ 1.18214701  0.49830501]
    Elastic net:  		 [ 1.17795855  0.49756084]
    true parameter values: 	 [1.2, 0.5]


## PyMC3を使う場合

上２つの場合では最尤推定を使ってパラメータを推定していましたが、PyMC3ではベイズ定理に基づいてマルコフ連鎖モンテカルロ法（MCMC)でパラメータの確率分布を求めます。


```python
import pymc3 as pm

with pm.Model() as model_robust:
    family = pm.glm.families.Normal()
    #pm.glm.glm('y ~ x', data, family=family)
    pm.glm.glm('y ~ x1+x2', df, family=family)
    start = pm.find_MAP()
    step = pm.NUTS(scaling=start)
    #step = pm.Metropolis()
    trace_robust = pm.sample(25000, step)
```

    Optimization terminated successfully.
             Current function value: 764.008811
             Iterations: 19
             Function evaluations: 30
             Gradient evaluations: 30


    100%|██████████| 25000/25000 [00:32<00:00, 761.52it/s]


```python
# show results
pm.traceplot(trace_robust[2000:])
#pm.summary(trace_robust[2000:])
plt.show()
pm.df_summary(trace_robust[2000:])
```
![output_9_0.png](https://qiita-image-store.s3.amazonaws.com/0/154129/abe56f1c-aa48-a75b-1d72-d6024454b1be.png)
　

<div>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>mean</th>
      <th>sd</th>
      <th>mc_error</th>
      <th>hpd_2.5</th>
      <th>hpd_97.5</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Intercept</th>
      <td>-0.022296</td>
      <td>0.040946</td>
      <td>0.000564</td>
      <td>-0.100767</td>
      <td>0.057460</td>
    </tr>
    <tr>
      <th>x1</th>
      <td>1.199371</td>
      <td>0.011235</td>
      <td>0.000126</td>
      <td>1.177191</td>
      <td>1.221067</td>
    </tr>
    <tr>
      <th>x2</th>
      <td>0.500502</td>
      <td>0.005346</td>
      <td>0.000057</td>
      <td>0.490258</td>
      <td>0.511003</td>
    </tr>
    <tr>
      <th>sd</th>
      <td>0.490271</td>
      <td>0.010949</td>
      <td>0.000072</td>
      <td>0.469599</td>
      <td>0.512102</td>
    </tr>
  </tbody>
</table>
</div>


## もうちょっと複雑な場合でやってみる

PyMC3がどこまでできるのか見るために、もうちょっと複雑な場合でやってみます。元になる式は以下：

$y = \beta_0 + \beta_1 x_1 + \beta_2 x_2+ \beta_3 x_3 + \beta_4 x_4 + \beta_5 x_5 + \epsilon$



```python
# Generate random data
beta = [1.2, 0.5, 9.8, 0.2, 1.08]
# prep data
x1 = np.random.random(size=1000)*5 
x2 = np.random.random(size=1000)*10 
x3 = np.random.random(size=1000)
x4 = np.random.normal(size=1000)*2 
x5 = np.random.random(size=1000)*60 
x = np.transpose([x1, x2, x3, x4, x5])
y = np.dot(x, beta)  + np.random.normal(scale=0.5, size=1000)
#data = dict(x=x, y=y)
data = dict(x1=x1, x2=x2, x3=x3, x4=x4, x5=x5, y=y)
df = pd.DataFrame(data)

with pm.Model() as model_robust:
    family = pm.glm.families.Normal()
    pm.glm.glm('y ~ x1+x2+x3+x4+x5', df, family=family)
    start = pm.find_MAP()
    step = pm.NUTS(scaling=start)
    #step = pm.Metropolis()
    trace_robust = pm.sample(25000, step)
```


```python
# show results
pm.traceplot(trace_robust[2000:])
plt.show()
print("true parameter values are:", beta)
pm.df_summary(trace_robust[2000:])
```
![output_13_0.png](https://qiita-image-store.s3.amazonaws.com/0/154129/13c85d7b-5046-25a6-4a17-c80a437df00e.png)


    true parameter values are: [1.2, 0.5, 9.8, 0.2, 1.08]


<div>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>mean</th>
      <th>sd</th>
      <th>mc_error</th>
      <th>hpd_2.5</th>
      <th>hpd_97.5</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Intercept</th>
      <td>0.041924</td>
      <td>0.059770</td>
      <td>0.000737</td>
      <td>-0.080421</td>
      <td>0.154130</td>
    </tr>
    <tr>
      <th>x1</th>
      <td>1.199973</td>
      <td>0.011466</td>
      <td>0.000106</td>
      <td>1.177061</td>
      <td>1.222395</td>
    </tr>
    <tr>
      <th>x2</th>
      <td>0.494488</td>
      <td>0.005656</td>
      <td>0.000053</td>
      <td>0.483624</td>
      <td>0.505661</td>
    </tr>
    <tr>
      <th>x3</th>
      <td>9.699889</td>
      <td>0.056527</td>
      <td>0.000484</td>
      <td>9.587273</td>
      <td>9.809352</td>
    </tr>
    <tr>
      <th>x4</th>
      <td>0.197271</td>
      <td>0.008424</td>
      <td>0.000052</td>
      <td>0.181196</td>
      <td>0.214320</td>
    </tr>
    <tr>
      <th>x5</th>
      <td>1.081120</td>
      <td>0.000922</td>
      <td>0.000008</td>
      <td>1.079339</td>
      <td>1.082917</td>
    </tr>
    <tr>
      <th>sd</th>
      <td>0.514316</td>
      <td>0.011438</td>
      <td>0.000067</td>
      <td>0.492296</td>
      <td>0.536947</td>
    </tr>
  </tbody>
</table>
</div>


# おわりに

短いコードで書けるので特に解説はしませんでしたが、どれもうまくパラメータを推定できていることがわかります。
もし何か質問があれば遠慮なくコメントください。

**（追記）** 出力画像が間違っていたので修正。