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

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

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

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

df = pd.read_csv("https://raw.githubusercontent.com/karaage0703/machine-learning-study/master/data/karaage_data.csv")
print(df.head())

最小二乗法

単回帰分析は、回帰式に含まれる２つの定数、すなわち回帰係数$a$と切片$b$を求めることが目標でした。
その際、より精度の高い単回帰式を得るためには、全体としての誤差、すなわち残差 $y - \hat{y}$ ができるだけ小さくなるように定数$a$、$b$が決定されねばなりません。
このの残差の定義を考えます。

$e_{1}+e_{2}+e_{3}…+e_{n}$
となりそうですが、これは誤りです。
実測値は、回帰直線に対し、正負それぞれの側に偏在しています。つまりプラスとマイナスとが互いに打ち消し合ってしまい、合計すると 0 になってしまうのです。

そこで、個体ごとの残差を二乗することで正負をなくし、単に隔たりの大きさ（距離）として扱えるようにします。
$Q = {e_{1}}^{2}+{e_{2}}^{2}+{e_{3}}^{2}…+{e_{n}}^{2}$
回帰直線からの距離の総量として$Q$が定義されました。
この $Q$ が最も小さくなることが、回帰直線の傾き$a$の決め手となり、それが得られれば$y$軸との交点$b$もおのずと求まります。
この方法を最小二乗法といいます。

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

⑶ 変数x, yそれぞれの平均値を算出する

mean_x = df['x'].mean()
mean_y = df['y'].mean()

⑷ 変数x, yそれぞれの偏差を算出する

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

# 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)

⑸ 変数xの分散を算出する

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

# 偏差平方和
ssdev_x = 0
for i in dev_x:
    d = i ** 2
    ssdev_x += d
# 分散
var_x = ssdev_x / (len(df) - 1)

⑹ 共分散を算出する

共分散$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) の数値で除算します。

# 偏差積和
spdev = 0
for i,j in zip(df['x'], df['y']):
    spdev += (i - mean_x) * (j - mean_y)
# 共分散
cov = spdev / (len(df) - 1)

⑺ 回帰係数aを算出する

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

a = cov / var_x

⑻ 切片bを算出する

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

b = mean_y - (a * mean_x)

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

⑼ 決定係数を算出し、回帰式の精度を確認する

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

# 予測値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)

決定係数もまた先の scikit-learn による計算結果と一致することが確認できました。

⑽ 散布図に回帰直線を合わせて図示する

plt.plot(x, y, "o") # 散布図
plt.plot(x, z, "r") # 回帰直線
plt.show()

ここまで、単回帰分析のアルゴリズムを学びました。
しかし、現実の社会では、ある現象をただ１つの要因で説明できるようなケースはまずありません。ある現象の背景には、さまざまな要因が大なり小なり同時にからみ合っているものです。
次に、３つ以上の変数をあつかう重回帰分析を学びます。