単回帰分析は、機械学習の１つの手法である、直線近似の手法です。統計学や実務のデータサイエンス分野の教科書でも、一番最初に登場するくらい基本的な分析手法です。
単回帰分析を実際に行う場合は、Pythonやエクセル、Rなどの数値計算ソフトを用いることで、一瞬で単回帰分析自体は行うことができます。一方で、単回帰分析を解析的に追っていくと、偏微分や関数の最適化などの話も登場し、少しとっつきにくい側面もあります。
一方で、単回帰分析に登場するような解析的な手法を理解できると、今後機械学習で登場する数式の理解度が飛躍的に上昇します。また今回の記事では取り扱いませんが、重回帰分析の導出などを学ぶと、線形代数の理解や知識が深まります。
今回の記事では、単回帰分析の理論的な側面や導出に重点をおいて、できるだけわかりやすく解説し、単回帰分析の手法やそこに登場する数学の概念を理解できるようにします。
また、実際の身長と体重のデータセットであるDavisデータセットを用い、体重から身長を単回帰分析する例題を通して、実務や研究で行うような単回帰分析の手法的な側面についても解説していきます。
単回帰分析の先は、重回帰分析や線形回帰モデルなどより複雑なモデルが登場しますが、最も基本的なモデルである、単回帰分析についてまずは理解をしていきましょう。
- 単回帰分析の概要を理解する
- 単回帰分析の導出について理解する
- Pythonで実際のデータセットを用いて単回帰分析を行う
単回帰分析の概要を理解する
まず最初に単回帰分析について簡単にイメージを持ちましょう。単回帰分析はこういうものかとわかってしまえば、本当に簡単です。
単回帰分析とは、1つの説明変数から1つの目的変数を予測する手法である。
例えば、身長・体重の例では、体重から身長を予測することや、家賃の例では、部屋の広さから、家賃を予測するといったものが単回帰分析です。
この場合、体重や部屋の広さといった項目が説明変数となり、身長や家賃が目的変数になります。一般的に説明変数が既知であり、目的変数の値を予測したい場合に、用います。
一方この記事では重回帰分析という言葉も登場していますが、重回帰分析においては、単回帰分析において1つだった説明変数が、2つ以上になった場合です。つまり、重回帰分析は単回帰分析の説明変数を拡張した場合に過ぎないので、まずは単回帰分析を理解することが重要です。単回帰分析自体は、非常に簡単な概念ですが、実際に計算が少し複雑です。
しかし、これを理解するのは、今後、さまざまな統計的手法・機械学習、ディープラーニング等を理解していく上で本当に必須なので、しっかり理解していきましょう。
ここで、先ほどと体重から身長を回帰する問題について考えます。
今、あるクラスメイトの体重と身長のデータが次のように得られているとします。（これは、後程説明するDavisデータセットを利用しています。）単回帰分析では、これらのデータから、次の図のような直線を得ることが目的です。
この直線は、一次関数で表現できます。
説明変数を$x$, 目的変数を$y$とした時、単回帰分析の回帰式は次の式で表現される。
\begin{equation}
y =wx + b
\end{equation}ここで、$w$と$b$はモデルパラメータである。($w, b \in \mathbb{R}$ )
(1)の回帰式に登場する、$w$と$b$はモデルパラーメータで、$w$はweight、$b$はbiasの頭文字になっています。単回帰分析の目的は、既にあるデータ（図中の緑色の点）から、この2つのモデルパラメータ$w$と$b$を見つけ、適切な回帰直線を見つけることです。
単回帰分析の導出を理解する
ここまででざっくりと単回帰分析の概要について話してきたので、次に、単回帰分析の数学的・解析的側面について、数式を交えながら説明していきます。（ここが正念場だと思います）
先ほど、単回帰分析のゴールは、データから適切なモデルパラメータを見つけることと書きました。では、どのように適切なモデルパラメータを見つけるのでしょうか。単回帰分析では、最小二乗法(least squares)を利用することで、最適なモデルパラメータを決定します。
単回帰分析における最小二乗法
最小二乗法とは、モデルからの予測値と実際の値の差（これを誤差(Error)と言います）の二乗の和(二乗和誤差, sum-of-squares error)を最小にするモデルパラメータを決定する方式です。
まず、誤差についてですが、今回学習用の$n$個のデータを$\mathcal{D}$とします。そして、$\mathcal{D}$のうち、m番目のデータ$(x_m, y_m)$に注目すると、誤差は次のように表現されます。
\begin{equation}
\begin{split}
(誤差) &= \hat{y}- y_m \\
&= (wx_m + b) - y_m
\end{split}
\end{equation}ここで、モデルの予測値を$\hat{y}$で表現しました。(2)の誤差の表現だと、符号がついてしまいますが、誤差の大きさ（絶対値)に注目したいので、通常は誤差を二乗して取り扱います。これが二乗誤差です。
\begin{equation}
\begin{split}
(二乗誤差) &= (\hat{y_m}- y_m)^2 \\
&=
\biggr \{ (wx_m + b) - y_m \biggr\}^2
\end{split}
\end{equation}$n$個のデータがある$\mathcal{D}$のうち、m番目のデータの誤差に注目すると、(3)になります。なので、$n$個のデータ全ての渡って足し合わせた、二乗和誤差が、最小化したい目的関数になります。
単回帰分析では、二乗和誤差を最小化するモデルパラメータが、最適なモデルパラメータになります。
一般的に、最適なモデルパラメータを決定する際に、最小化する関数のことを、目的関数(
objective function)と言います。単回帰分析の目的関数$$は次のようになります。
\begin{equation}
\begin{split}
\mathcal{L} &=
\sum_{n=1}^{N} (\hat{y_n}- y_n)^2 \\
&=
\sum_{n=1}^{N} \biggr \{ (wx_n + b) - y_n \biggr\}^2
\end{split}
\end{equation}単回帰分析で最適なパラメータを見つける
単回帰分析においては、目的関数$\mathcal{L}$を最小化させることが、最適なパラメータとなります。
では、どのように目的関数を最小化すれば良いでしょうか。ここで、(4)式について、整理すると、
\begin{equation}
\begin{split}
\mathcal{L} &=
\sum_{n=1}^{N} \biggr \{ (wx_n + b) - y_n \biggr\}^2 \\
&= \sum_{n=1}^{N} (w^2x_n^2 + b^2 + 2wbx_n -2wy_nx_n -2by_n +y_n^2)
\end{split}
\end{equation}
となります。(5)式をパラメータ$b$と$w$について注目すると、それぞれ$b$と$w$の、下に凸の2次関数になっていることがわかります。なぜなら、(5)式において$\sum$の中身が、$b$, $w$に対して、下に凸の2次関数であり、それらをいくつ足しても、下に凸な2次関数であることに変わりはないからです。
下に凸な2次関数の場合、1次微分の値がゼロになるときに関数が最小化されます。つまり、
\begin{equation}
\begin{split}
\frac{\partial}{\partial w} \mathcal{L} =0 \\
\frac{\partial}{\partial b} \mathcal{L} =0 \\
\end{split}
\end{equation}(6)式を満たす、$w$と$b$が求めるパラメータとなります。
(6)を計算すると、
\begin{equation}
\begin{split}
\frac{\partial}{\partial w} \mathcal{L} &= 2wx_n^2 + 2bx_n -2 y_nx_n \\
&= 2 \sum_{n=1}^{N} ( w {x_n}^2 + bx_n- y_nx_n ) = 0 \\ \\
\frac{\partial}{\partial b} \mathcal{L} &= 2b + 2wx_n -2 y_n \\
&= 2 \sum_{n=1}^{N} (b + wx_n -y_n) =0
\end{split}
\end{equation}を満たす、$w$と$b$が求めるべき、パラメータとなります。これをさらに整理すると、
\begin{equation}
\begin{split}
\sum_{n=1}^{N} w{x_n}^2 + b\sum_{n=1}^{N} x_n &= \sum_{n=1}^{N} x_n y_n \\
bN + w\sum_{n=1}^{N} x_n &= \sum_{n=1}^{N} y_n
\end{split}
\end{equation}これらは、$w$と$b$の連立方程式になっています。パラメータ2つで式が2本の連立方程式なので、解くことはできます。(8)式の下の式より、
\begin{equation}
b = \frac{1}{N} \biggl \{ \sum_{n=1}^{N} y_n - w \sum_{n=1}^{N} x_n \biggr \}
\end{equation}これを、(8)式の上の式に代入すると、
\begin{split}
\sum_{n=1}^{N} w{x_n}^2 + \frac{1}{N} \biggl \{ \sum_{n=1}^{N} y_n - w \sum_{n=1}^{N} x_n \biggr \}\sum_{n=1}^{N} x_n &= \sum_{n=1}^{N} x_n y_n
\end{split}
これを成立すると、
\begin{split}
w \biggl \{
\sum_{n=1}^{N} {x_n}^2 -\frac{1}{N} \biggl (\sum_{n=1}^{N} x_n \biggr )^2
\biggr \} = \sum_{n=1}^{N} x_n y_n - \frac{1}{N} \sum_{n=1}^{N} y_n \sum_{n=1}^{N} x_n
\end{split}となるので、
\begin{equation}
\begin{split}
w = \frac
{ N \sum_{n=1}^{N} x_n y_n - \sum_{n=1}^{N} y_n \sum_{n=1}^{N} x_n}
{ N \sum_{n=1}^{N} {x_n}^2 - (\sum_{n=1}^{N} x_n)^2 }
\end{split}
\end{equation}となります。また、これを(9)に代入すると、
\begin{equation}
\begin{split}
b = \frac
{ \sum_{n=1}^{N} x_n^2 \sum_{n=1}^{N} y_n - \sum_{n=1}^{N} x_ny_n \sum_{n=1}^{N} x_n}
{ N \sum_{n=1}^{N} {x_n}^2 - (\sum_{n=1}^{N} x_n)^2 }
\end{split}
\end{equation}となります。
(8)の連立方程式を代入法を用いて解きましたが、実際には、行列計算でも解くことができます。(8)の連立方程式を行列で表現すると、
\begin{equation}
\begin{split}
\begin{pmatrix}
N & \sum_{n=1}^{N} x_n \\
\sum_{n=1}^{N} x_n & \sum_{n=1}^{N} x_n^2
\end{pmatrix}
\begin{pmatrix}
b \\
w
\end{pmatrix}
=
\begin{pmatrix}
\sum_{n=1}^{N} y_n \\
\sum_{n=1}^{N} x_n y_n
\end{pmatrix}
\end{split}
\end{equation}このように表現できます。これを、単回帰分析の正規方程式(normal equation)といいます。正規方程式を解くことで、単回帰のパラメータを求めることができます。
この正規方程式を特には、(12)式の左辺の行列に注目します。この2次の正方行列の逆行列は、次で表されます。
\begin{equation}
\begin{pmatrix}
N & \sum_{n=1}^{N} x_n \\
\sum_{n=1}^{N} x_n & \sum_{n=1}^{N} x_n^2
\end{pmatrix}^{-1} =
\frac{1}{N\sum_{n=1}^{N} x_n^2 - (\sum_{n=1}^{N} x_n)^2 }
\begin{pmatrix}
\sum_{n=1}^{N} x_n^2 & - \sum_{n=1}^{N} x_n \\
- \sum_{n=1}^{N} x_n & N
\end{pmatrix}
\end{equation}ここでは、2次の正方行列の逆行列の公式を使いました。この公式の詳しい説明は、下記の記事で扱っています。
(13)で得ることができた逆行列を(12)式に左からかけると、
\begin{equation}
\begin{split}
\begin{pmatrix}
b \\
w
\end{pmatrix}
&=
\frac{1}{N\sum_{n=1}^{N} x_n^2 - (\sum_{n=1}^{N} x_n)^2 }
\begin{pmatrix}
\sum_{n=1}^{N} x_n^2 & - \sum_{n=1}^{N} x_n \\
- \sum_{n=1}^{N} x_n & N
\end{pmatrix}
\begin{pmatrix}
\sum_{n=1}^{N} y_n \\
\sum_{n=1}^{N} x_n y_n
\end{pmatrix} \\
&=
\frac{1}{N\sum_{n=1}^{N} x_n^2 - (\sum_{n=1}^{N} x_n)^2 }
\begin{pmatrix}
\sum_{n=1}^{N} x_n^2 \sum_{n=1}^{N} y_n - \sum_{n=1}^{N} x_ny_n \sum_{n=1}^{N} x_n \\
N \sum_{n=1}^{N} x_n y_n - \sum_{n=1}^{N} y_n \sum_{n=1}^{N} x_n
\end{pmatrix} \\
\end{split}
\end{equation}となり、(10)、(11)と同じ結果を得ることができました。
Pythonで実際のデータセットを用いて単回帰分析を行う
ここまで、数式をいじりながら、解析的に単回帰分析の最適なパラメータの計算をしてきました。どうでしたか？単純な単回帰であっても、難しいと思った人も多いと思います。
実際に実務で単回帰分析をする際は、sckit-learnなどのライブラリが良い感じにやってくれるので、ここまで苦労することはありませんが、まず解析的に数式を理解することは、今後新しいアルゴリズムや論文等を読む際に必ず必要な力になるので、できる限り理解できると良いでしょう。
さて、ここからは、実際にPythonを用いて、実際のデータセットで単回帰分析をいきます。
今回扱うデータセットは、記事の冒頭で取り扱った、身長と体重の2つの項目があるDavisデータセットを利用し、体重データから身長を予測するという単回帰回を考えることにします。
データセットの準備
ここからはPythonを用いてデータセットを準備していきます。
まず、次のコードで、今回用いる必要なライブラリと、Davisデータセットをダウンロードして、pandasのDataFrameに読み込みます。
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
df = pd.read_csv('https://vincentarelbundock.github.io/Rdatasets/csv/carData/Davis.csv', index_col=False)
df = df[['weight', 'height']]
mask = ((df['weight'] < 140) & (df['height'] > 120)) # 正常なデータの範囲を指定することで外れ値を除外
df = df[mask]
data = df.to_numpy()Davisデータセットは、外れ値が含まれるので、今回はそれを除外しています。Davisデータセットをとり扱う方法については、下記の記事で詳しく解説しているため、今回は省略します。
データを確認してみます。
fig, ax = plt.subplots(dpi=100, figsize=(5,5))
ax.set_facecolor('w')
ax.set_xlabel('weight')
ax.set_ylabel('height')
ax.set_aspect('equal')
ax.set_xlim([0, 200])
ax.set_ylim([0, 200])
ax.scatter(data[:, 0], data[:, 1], s=10, marker="o", facecolor='#35966e')
plt.grid(color='gray', linestyle='dotted', linewidth=1)
plt.show()
このような分布になることがわかります。ここで実際に単回帰分析を行い、直線の方程式を求めてみます。
scikit-learnで単回帰分析
単回帰分析はscikit-learnで行う場合、非常に簡単です。次のコードで、データに対して直線をフィッティングすることができます。
from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(data[:, 0].reshape(-1, 1), data[:, 1].reshape(-1, 1))得られた、直接の傾きと切片は、それぞれ傾きが、model.coef_ 、切片がmodel.intercept_で得ることができています。
print("y = {} + {}".format(model.coef_[0][0] , model.intercept_[0]))
# => y = 0.5168935754352395 + 136.83660744117836結果はこのようにありました。傾きが、0.5168935754352395ということなので、体重が1kg増加するごとに、身長が0.5cm程度増加するという結果が、単回帰分析によって得ることができました。
得られた直線を図示すると次のようになります。
fig, ax = plt.subplots(dpi=100, figsize=(5,5))
ax.set_facecolor('w')
ax.set_xlabel('weight')
ax.set_ylabel('height')
ax.set_aspect('equal')
ax.set_xlim([0, 200])
ax.set_ylim([0, 200])
xs = np.linspace(0, 200, 10)
pridict = model.coef_[0] * xs + model.intercept_[0]
ax.scatter(data[:, 0], data[:, 1], s=10, marker="o", facecolor='#35966e')
ax.plot(xs, pridict, linestyle="solid")
plt.grid(color='gray', linestyle='dotted', linewidth=1)
plt.show()
うまくデータの傾向を捉えることができていますね。