線形単回帰を Python で理解する


# 準備

## 必要なライブラリのインポート


```python
# 数値計算やデータフレーム操作に関するライブラリをインポートする
import numpy as np
import pandas as pd
```


```python
# URL によるリソースへのアクセスを提供するライブラリをインポートする。
import urllib.request 
```


```python
# 図やグラフを図示するためのライブラリをインポートする。
import matplotlib.pyplot as plt
%matplotlib inline
```


```python
# 線形回帰を行なうライブラリ
from sklearn import linear_model
```

## データのダウンロードと整形

### 赤ワインのデータ

データは <a href="http://archive.ics.uci.edu/ml/index.php">UC Irvine Machine Learning Repository</a> から取得したものを少し改変しました。

赤ワイン https://raw.githubusercontent.com/chemo-wakate/tutorial-6th/master/beginner/data/winequality-red.txt

白ワイン https://raw.githubusercontent.com/chemo-wakate/tutorial-6th/master/beginner/data/winequality-white.txt

* fixed acidity : 不揮発酸濃度（ほぼ酒石酸濃度）
* volatile acidity : 揮発酸濃度（ほぼ酢酸濃度）
* citric acid : クエン酸濃度
* residual sugar : 残存糖濃度
* chlorides : 塩化物濃度
* free sulfur dioxide : 遊離亜硫酸濃度
* total sulfur dioxide : 亜硫酸濃度
* density : 密度
* pH : pH
* sulphates : 硫酸塩濃度
* alcohol : アルコール度数
* quality (score between 0 and 10) : 0-10 の値で示される品質のスコア


```python
# ウェブ上のリソースを指定する
url = 'https://raw.githubusercontent.com/chemo-wakate/tutorial-6th/master/beginner/data/winequality-red.txt'
# 指定したURLからリソースをダウンロードし、名前をつける。
urllib.request.urlretrieve(url, 'winequality-red.txt') 
```




    ('winequality-red.txt', <http.client.HTTPMessage at 0x11de96198>)




```python
# データの読み込み
df = pd.read_csv('winequality-red.txt', sep='\t', index_col=0)
```


```python
df.head()
```




<div>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>fixed acidity</th>
      <th>volatile acidity</th>
      <th>citric acid</th>
      <th>residual sugar</th>
      <th>chlorides</th>
      <th>free sulfur dioxide</th>
      <th>total sulfur dioxide</th>
      <th>density</th>
      <th>pH</th>
      <th>sulphates</th>
      <th>alcohol</th>
      <th>quality</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>7.4</td>
      <td>0.70</td>
      <td>0.00</td>
      <td>1.9</td>
      <td>0.076</td>
      <td>11.0</td>
      <td>34.0</td>
      <td>0.9978</td>
      <td>3.51</td>
      <td>0.56</td>
      <td>9.4</td>
      <td>5</td>
    </tr>
    <tr>
      <th>1</th>
      <td>7.8</td>
      <td>0.88</td>
      <td>0.00</td>
      <td>2.6</td>
      <td>0.098</td>
      <td>25.0</td>
      <td>67.0</td>
      <td>0.9968</td>
      <td>3.20</td>
      <td>0.68</td>
      <td>9.8</td>
      <td>5</td>
    </tr>
    <tr>
      <th>2</th>
      <td>7.8</td>
      <td>0.76</td>
      <td>0.04</td>
      <td>2.3</td>
      <td>0.092</td>
      <td>15.0</td>
      <td>54.0</td>
      <td>0.9970</td>
      <td>3.26</td>
      <td>0.65</td>
      <td>9.8</td>
      <td>5</td>
    </tr>
    <tr>
      <th>3</th>
      <td>11.2</td>
      <td>0.28</td>
      <td>0.56</td>
      <td>1.9</td>
      <td>0.075</td>
      <td>17.0</td>
      <td>60.0</td>
      <td>0.9980</td>
      <td>3.16</td>
      <td>0.58</td>
      <td>9.8</td>
      <td>6</td>
    </tr>
    <tr>
      <th>4</th>
      <td>7.4</td>
      <td>0.70</td>
      <td>0.00</td>
      <td>1.9</td>
      <td>0.076</td>
      <td>11.0</td>
      <td>34.0</td>
      <td>0.9978</td>
      <td>3.51</td>
      <td>0.56</td>
      <td>9.4</td>
      <td>5</td>
    </tr>
  </tbody>
</table>
</div>



### 品質の良いワインだけ考える


```python
# 品質の良いワインだけ選択し df1 に入れる
df1 = df[df['quality'] > 6]
```


```python
df1.head()
```




<div>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>fixed acidity</th>
      <th>volatile acidity</th>
      <th>citric acid</th>
      <th>residual sugar</th>
      <th>chlorides</th>
      <th>free sulfur dioxide</th>
      <th>total sulfur dioxide</th>
      <th>density</th>
      <th>pH</th>
      <th>sulphates</th>
      <th>alcohol</th>
      <th>quality</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>7</th>
      <td>7.3</td>
      <td>0.65</td>
      <td>0.00</td>
      <td>1.2</td>
      <td>0.065</td>
      <td>15.0</td>
      <td>21.0</td>
      <td>0.9946</td>
      <td>3.39</td>
      <td>0.47</td>
      <td>10.0</td>
      <td>7</td>
    </tr>
    <tr>
      <th>8</th>
      <td>7.8</td>
      <td>0.58</td>
      <td>0.02</td>
      <td>2.0</td>
      <td>0.073</td>
      <td>9.0</td>
      <td>18.0</td>
      <td>0.9968</td>
      <td>3.36</td>
      <td>0.57</td>
      <td>9.5</td>
      <td>7</td>
    </tr>
    <tr>
      <th>16</th>
      <td>8.5</td>
      <td>0.28</td>
      <td>0.56</td>
      <td>1.8</td>
      <td>0.092</td>
      <td>35.0</td>
      <td>103.0</td>
      <td>0.9969</td>
      <td>3.30</td>
      <td>0.75</td>
      <td>10.5</td>
      <td>7</td>
    </tr>
    <tr>
      <th>37</th>
      <td>8.1</td>
      <td>0.38</td>
      <td>0.28</td>
      <td>2.1</td>
      <td>0.066</td>
      <td>13.0</td>
      <td>30.0</td>
      <td>0.9968</td>
      <td>3.23</td>
      <td>0.73</td>
      <td>9.7</td>
      <td>7</td>
    </tr>
    <tr>
      <th>62</th>
      <td>7.5</td>
      <td>0.52</td>
      <td>0.16</td>
      <td>1.9</td>
      <td>0.085</td>
      <td>12.0</td>
      <td>35.0</td>
      <td>0.9968</td>
      <td>3.38</td>
      <td>0.62</td>
      <td>9.5</td>
      <td>7</td>
    </tr>
  </tbody>
</table>
</div>



### pH と不揮発酸濃度との関係を考える


```python
X = df1['pH'].values
Y = df1['fixed acidity'].values
```


```python
plt.title('Good wines')
plt.scatter(X, Y)
plt.xlabel('pH')
plt.ylabel('fixed acidity')
plt.grid()
plt.show()
```


![output_14_0.png](https://qiita-image-store.s3.amazonaws.com/0/211162/3198ca55-3ef8-dbbd-3b9e-7bdc1e40d50f.png)



# scikit-learn を用いた線形単回帰分析


```python
from sklearn import linear_model
clf = linear_model.LinearRegression()
X2 = [[x] for x in X]
clf.fit(X2, Y) # 予測モデルを作成

# 散布図
plt.scatter(X2, Y)
 
# 回帰直線
plt.title('Linear regression')
plt.plot(X2, clf.predict(X2))
plt.xlabel('pH')
plt.ylabel('fixed acidity')
plt.grid()
plt.show()

print("回帰係数= ", clf.coef_)
print("切片= ", clf.intercept_)
print("決定係数= ", clf.score(X2, Y))
```


![output_16_1.png](https://qiita-image-store.s3.amazonaws.com/0/211162/01bbf260-b6df-d1ce-03a9-3b672dad6d09.png)



    回帰係数=  [-9.985041]
    切片=  41.685825856217654
    決定係数=  0.5948141769596296


# 演習：上の結果を scikit-learn を使わず再現しよう

## 平均値


```python
# 平均値を求める関数を作ろう
def average(list):

```


```python
print(average(X))
print(average(Y))
```

    3.2888018433179704
    8.84700460829493


## 分散


```python
# 分散を求める関数を作ろう
def variance(list):

```


```python
print(variance(X))
print(variance(Y))
```

    0.023753403130242778
    3.9814772027437417


## 標準偏差


```python
# 標準偏差を求める関数を作ろう
import math
def standard_deviation(list):

```


```python
print(standard_deviation(X))
print(standard_deviation(Y))
```

    0.1541213908912153
    1.9953639273936326


## 共分散


```python
# 共分散 = 偏差積の平均 （偏差値、ではありません。偏差積、です）を作ろう
def covariance(list1, list2): 

```


```python
print(covariance(X, Y))
```

    -0.23717870415596012


## 相関係数


```python
# 相関係数 = 共分散を list1, list2 の標準偏差で割ったものを作ろう
def correlation(list1, list2):

```


```python
print(correlation(X, Y))
```

    -0.7712419704344603


## 回帰直線


```python
# 回帰直線の傾き＝相関係数＊（（yの標準偏差）／（xの標準偏差））を求める関数を作ろう
def a_fit(xlist, ylist):

```


```python
print(a_fit(X, Y))
```

    -9.985041000461308



```python
# y切片＝yの平均－（傾き＊xの平均）を求める関数を作ろう
def b_fit(xlist, ylist):

```


```python
print(b_fit(X, Y))
```

    41.68582585621759



```python
# 回帰直線の式
print("y = ax + b; (a, b) = ({}, {})".format(a_fit(X, Y), b_fit(X, Y)))
```

    y = ax + b; (a, b) = (-9.985041000461308, 41.68582585621759)



```python
# 回帰式を作ろう
def f(x):

```


```python
f(3.5)
```




    6.738182354603012



## 決定係数


```python
# 決定係数を求める関数を作ろう
def r2(X, Y):

```


```python
r2(X,Y)
```




    0.5948141769596293


```python
# 散布図と回帰直線を描いてみよう
```

![output_43_0.png](https://qiita-image-store.s3.amazonaws.com/0/211162/e76278bc-2429-2f07-13dd-4b33715ed39e.png)



    回帰係数=  -9.985041000461308
    切片=  41.68582585621759
    決定係数=  0.5948141769596293

# 補足：計算速度の改善の話

講義中に「回帰式を作ろう」の部分で

```python
def f(x):
    return a_fit(X, Y) * x + b_fit(X, Y)
```

と書いてしまいましたが、それは良くなかったですね。回帰係数と切片は定数なので、本来は一度だけ計算すれば済むものなのですが、これだと、f(x) を呼び出すたびに回帰係数と切片を再計算しなければならなくなるため、決定係数の計算が大変遅くなります。試しに品質の良い赤ワインのデータで時間を計測してみると

```python
%timeit r2(X,Y)
```

    328 ms ± 9.09 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)


もかかってしまいました。そこで、回帰係数と切片を再計算しなくて済むように改良すると

```python
afit = a_fit(X, Y)
bfit = b_fit(X, Y)
def f(x):
    return afit * x + bfit
```

```python
%timeit r2(X,Y)
```

    6.68 ms ± 297 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

だいぶ改善されました。

# 課題

以下の例について、scikit-learnを使わずに線形単回帰を行ない、散布図と回帰直線を描画し、回帰係数・切片・決定係数を表示してください。また、scikit-learnを使った結果と比較することにより検算してください。

* 品質の良い赤ワインの密度(density)から不揮発酸濃度(fixed acidity)を予測
* 品質の良くない赤ワインのpHから不揮発酸濃度(fixed acidity)を予測
* 品質の良くない赤ワインの密度 (density) から不揮発酸濃度 (fixed acidity) を予測
* 品質の良い白ワインのpHから不揮発酸濃度(fixed acidity)を予測
* 品質の良い白ワインの密度(density)から不揮発酸濃度(fixed acidity)を予測
* 品質の良くない白ワインのpHから不揮発酸濃度(fixed acidity)を予測
* 品質の良くない白ワインの密度(density)から不揮発酸濃度(fixed acidity)を予測

ただし、「品質が良くない」 = 「quality が 5 未満」とします。

以上の解析から、赤ワイン・白ワインの品質と、pH、密度、不揮発酸濃度の関係について考察してください。