2. Pythonで綴る多変量解析　1-2. 単回帰分析(アルゴリズム)

基本的な数値計算に必要なNumpyとPandasだけを使って、単回帰式を求めていきます。

####**⑴ ライブラリをインポートする**

```python
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
```

####**⑵ データを取り込み、中身を確認する**

```python
df = pd.read_csv("https://raw.githubusercontent.com/karaage0703/machine-learning-study/master/data/karaage_data.csv")
print(df.head())
```
![002_002_001.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/3d24fb18-c140-74a6-730a-449f11a492a8.png)

##**<font color="Gray">最小二乗法</font>**

単回帰分析は、回帰式に含まれる２つの定数、すなわち回帰係数$a$と切片$b$を求めることが目標でした。
その際、より精度の高い単回帰式を得るためには、全体としての誤差、すなわち残差 $y - \hat{y}$ ができるだけ小さくなるように定数$a$、$b$が決定されねばなりません。
このの**残差の定義**を考えます。
![002_002_002.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/9ae97b16-40d7-6355-2f54-b36d1342e19f.png)
$e_{1}+e_{2}+e_{3}…+e_{n}$
となりそうですが、これは誤りです。
実測値は、回帰直線に対し、正負それぞれの側に偏在しています。つまりプラスとマイナスとが互いに打ち消し合ってしまい、合計すると 0 になってしまうのです。
![002_002_003.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/cae041e1-ddc0-e174-c6c2-64254a40fadd.png)
そこで、個体ごとの残差を二乗することで正負をなくし、単に隔たりの大きさ（距離）として扱えるようにします。
$Q = {e_{1}}^{2}+{e_{2}}^{2}+{e_{3}}^{2}…+{e_{n}}^{2}$
回帰直線からの距離の総量として$Q$が定義されました。
この $Q$ が最も小さくなることが、回帰直線の傾き$a$の決め手となり、それが得られれば$y$軸との交点$b$もおのずと求まります。
この方法を**最小二乗法**といいます。


最小二乗法にもとづいて、単回帰式を解いていきます。


####**⑶ 変数x, yそれぞれの平均値を算出する**
```python
mean_x = df['x'].mean()
mean_y = df['y'].mean()
```
![002_002_004.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/c9d083fd-6e5c-5da7-de3d-2106711a2411.png)


####**⑷ 変数x, yそれぞれの偏差を算出する**
偏差とは、各個体の数値と平均値との差です。変数$x$については$x_{i}-\bar{x}$、変数$y$については$y-\bar{y}$を算出します。各変数ともデータの個数分を計算することになります。

```python
# xの偏差
dev_x = []
for i in df['x']:
    dx = i - mean_x
    dev_x.append(dx)
# yの偏差
dev_y = []
for j in df['y']:
    dy = j - mean_y
    dev_y.append(dy)
```
![002_002_005.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/218f60b5-b58e-6a2c-fbe2-3a3b4e245d21.png)


####**⑸ 変数xの分散を算出する**
⑷で求めた偏差をつかって分散を計算します。分散とは、偏差の二乗の平均をとったもの、すなわち偏差ごとに二乗して、その総和を (データの個数-1) の数値で除算します。

```python
# 偏差平方和
ssdev_x = 0
for i in dev_x:
    d = i ** 2
    ssdev_x += d
# 分散
var_x = ssdev_x / (len(df) - 1)
```
![002_002_006.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/78719353-880d-bf7d-98a8-5c8e9d0930ff.png)



####**⑹ 共分散を算出する**
共分散$s_{xy}$とは、２変数の関係の強さを表わす指標の一つで、次式で定義されます。
$s_{xy} = \frac{1}{n - 1} \displaystyle \sum_{i = 1}^n 
{(x_i - \overline{x})(y_{i} - \overline{y})}$
データの個体ごとの組を考えます。$(x_{1}, y_{1}), (x_{2}, y_{2}), ..., (x_{n}, y_{n})$という$n$組があるとき、組ごとに$x$の偏差と$y$の偏差を乗算し、それらの総和を (データ個数-1) の数値で除算します。

```python
# 偏差積和
spdev = 0
for i,j in zip(df['x'], df['y']):
    spdev += (i - mean_x) * (j - mean_y)
# 共分散
cov = spdev / (len(df) - 1)
```
![002_002_007.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/bd20cb74-0539-74ad-b77a-20b6810ed238.png)


####**⑺ 回帰係数aを算出する**
最小二乗法によって回帰係数を求める公式を示します。
$a = \frac{S_{xy}}{Sx^2}$
⑹で求めた共分散$S_{xy}$を、⑸で求めた変数$x$の分散$Sx^2$で除算することで回帰係数$a$が求められます。

```python
a = cov / var_x
```
![002_002_008.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/db277fcb-f24b-a2d1-8f88-bc379d997071.png)


####**⑻ 切片bを算出する**
単回帰式 $y = ax + b$ を変形し、$b = y -ax$として、⑶で求めた平均値$\bar{x}, \bar{y}$と、⑺で求めた回帰係数$a$を代入します。

```python
b = mean_y - (a * mean_x)
```
![002_002_009.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/b7147dd4-6236-f726-d292-b534cb9ffdcf.png)


**以上、最小二乗法の公式により単回帰式が求められました。**
**先に、機械学習ライブラリ scikit-learn を利用して得られた計算結果と一致しています。そこで、決定係数の確認もまた自前で計算してみます。**

####**⑼ 決定係数を算出し、回帰式の精度を確認する**
回帰式をつかって予測値のデータを作成し、その分散を求めます。
それが、実測値$y$の分散に対して何割を占めるか、すなわちどの程度もとの変量$y$を説明できているか。

```python
# 予測値zのデータ作成
df['z'] = (a * df['x']) + b
print(df)

# 予測値zの分散
ssdev_z = 0
for i in df['z']:
    j = (i - df['z'].mean())**2
    ssdev_z += j
var_z = ssdev_z / (len(df) - 1)
print("予測値の分散:", var_z)

# 実測値yの分散
ssdev_y = 0
for i in dev_y:
    j = i ** 2
    ssdev_y += j
var_y = ssdev_y / (len(df) - 1)
print("実測値yの分散:", var_y)

# 決定係数
R = var_z / var_y
print("決定係数R:", R)
```
![002_002_010.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/52f6a2a5-bfd5-154a-2248-7193381b30ee.png)
決定係数もまた先の scikit-learn による計算結果と一致することが確認できました。



####**⑽ 散布図に回帰直線を合わせて図示する**
```python
plt.plot(x, y, "o") # 散布図
plt.plot(x, z, "r") # 回帰直線
plt.show()
```
![002_002_011.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/502109/638c3f01-d969-7473-8d80-0c4d31d7c56e.png)
---
ここまで、単回帰分析のアルゴリズムを学びました。
しかし、現実の社会では、ある現象をただ１つの要因で説明できるようなケースはまずありません。ある現象の背景には、さまざまな要因が大なり小なり同時にからみ合っているものです。
次に、３つ以上の変数をあつかう重回帰分析を学びます。