【Python実装】重回帰分析の理論と実装を理解する

多変量解析の中で、重回帰分析は非常によく用いられている手法です。この記事を読んでいる人の中でも、その名前くらいは耳にしたことも多いのではないでしょうか。

重回帰分析は、説明変数が一般に多次元あるような場合の回帰問題を解く手法で、説明変数が1つである単回帰分析を単純に多次元に拡張したような手法となっています。（単回帰分析について扱っている記事はこちらをご参照ください）

実際のでデータ分析で重回帰分析は非常によく利用されます。実際にPythonを用いる場合は、ほんの10行足らずのコードで分析自体はできるのですが、重回帰分析は、実際にコードを書いて実行できるのと同じくらい、手法について数式レベルで理解していることが重要です。

重回帰分析の場合、多数のデータに対して多次元の説明変数を用いることから、その導出に行列表現が頻出しますが、どの式変形も機械学習やデータサイエンス分野で頻出のため、このあたりの数式を理解することは、大学初頭レベルの数学の復習に大変役に立ちます。

今回の記事では、重回帰分析について、その解析的な解の導出を数式を交えて実際に証明します。また、記事の後半では、実際に得られた解とデータセットを用いて、重回帰分析をPythonで実施します。

本記事の内容を理解することで、重回帰分析についてかなり理解度を上げることができるだけでなく、数式についても理解できるようになると思います。

本記事の内容

重回帰分析の概要

重回帰分析の概要について紹介します。まず、重回帰分析を学ぶには、単回帰分析について理解していない人は、まずそちらを学習することをお勧めします。

重回帰分析の理論に入っていきます。重回帰分析では、多次元の説明変数を扱うので、$D$次元の説明変数$\bm{x} = { x_1, x_2, \cdots, x_D }$と表現します。

また学習のデータセットの数は$N$個とし、全部のデータの集合は、$\mathcal{D} = { (\bm{x_1}, y_1), (\bm{x_2}, y_2), \cdots , (\bm{x_N}, y_N)}$と表現します。

重回帰分析の回帰式は次のようになります。

重回帰分析の定義式

\begin{equation} \begin{split} \hat{y} &= w_0 + w_1x_1 + w_2 x_2 + \cdots + x_Dw_D \\ &= \bm{w}^T \bm{x} \end{split} \end{equation}

ここで、$\bm{w}$は重みパラメータのベクトルで、$\bm{w} = \{w_0, w_1, \cdots, w_d \}$である。

重回帰分析のゴールは、この重みパラメータ$\bm{w} = \{w_0, w_1, \cdots, w_d \}$をデータ$\mathcal{D} $から求めることになります。

ここで、(1)を全てのデータに対して当てはめることを考えると

\begin{equation} \begin{split} \bm{\hat{y}} = \begin{pmatrix} \hat{y_1} \\ \hat{y_2} \\ \vdots \\ \hat{y_D} \\ \end{pmatrix} & = \begin{pmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1D} \\ 1 & x_{21} & x_{22} & \cdots &x_{2D} \\ \vdots & \vdots & \vdots & \vdots & \ddots \\ 1 & x_{n1} & x_{n2} & \cdots & x_{nD} \\ \end{pmatrix} \begin{pmatrix} w_0 \\ w_1 \\ \vdots \\ w_d \\ \end{pmatrix} \\ &= \bm{X} \bm{w} \end{split} \end{equation}

とかけることがわかります。ここで、$x_{nd}$は、n番目のデータセットのd番目の説明変数を意味しています。重回帰分析の導出で、(2)の式変形を見ると、最初はギョッとしますが、実際に手を動かしてみると、単純な話なので、ぜひ式変形をしてみてください。

(2)で登場する$\bm{X}$は、行にデータを並べ、列方向に説明変数を並べています。このような行列はデザイン行列や計画行列と呼ばれ、機械学習や統計分野でよく登場する形なので、覚えておいた方が良いでしょう。デザイン行列については、下記の記事で解説しています。

ここまでが、重回帰分析における準備段階です。

では、続いて、$\bm{w} = \{w_0, w_1, \cdots, w_d \}$を決定するにはどうしらいいでしょうか？重回帰分析では、単回帰分析と同様に、モデルの予測値$\hat{y_n}$と正解値$y_n$の二乗和誤差を最小化する、最小二乗法によってパラメータを決定します。

二乗和誤差を$E$とすると、二乗和誤差は次のようにかけます。

\begin{equation} \begin{split} E &= \sum_{n=1}^{N} (\hat{y_n} - y_n)^2 = \sum_{n=1}^{N}(y_n - \bm{w^T}\bm{x})^2 \\ &= (\bm{y} - \bm{X} \bm{w} )^T (\bm{y} - \bm{X} \bm{w} ) \\ &= (\bm{y}^T - (\bm{X} \bm{w})^T ) (\bm{y} - \bm{X} \bm{w} ) \\ &= \bm{y}^T \bm{y} -\bm{y}^T \bm{X} \bm{w} - (\bm{X} \bm{w})^T\bm{y} + (\bm{X}\bm{w})^T\bm{X} \bm{w} \\ &= \bm{y}^T \bm{y} - \bm{w}^T \bm{X}^T \bm{y} - \bm{w}^T \bm{X}^T \bm{y} + \bm{w}^T \bm{X}^T\bm{X} \bm{w} \\ &= \bm{y}^T \bm{y} - 2 \bm{w}^T \bm{X}^T \bm{y} + \bm{w}^T \bm{X}^T\bm{X} \bm{w} \end{split} \end{equation}

となります。ここで、誤差$E$をパラメータ$\bm{w}$で微分することを考えます。

(3)の二乗和誤差は、いろいろと行列の計算をしていますが、最終的に得られる値はスカラー値です。そのため、スカラー値をベクトルで偏微分する必要性があります。

(3)の二乗和誤差を$\bm{w}$で偏微分すると、

\begin{equation} \begin{split} \frac{\partial}{\partial \bm{w}} E &= -2\bm{X^T}\bm{y} +2(\bm{X^T}\bm{X})\bm{w} =0 \end{split} \end{equation}

よって、

\begin{equation} \begin{split} (\bm{X^T}\bm{X})\bm{w} = \bm{X^T}\bm{y} \end{split} \end{equation}

となります。これは、重回帰モデルの正規方程式と呼ばれています。

ここで、$\bm{X^T}\bm{X}$が逆行列$(\bm{X^T}\bm{X})^{-1}$を持つとき、(5)の正規方程式の両辺に左から$(\bm{X^T}\bm{X})^{-1}$をかけることで、

\begin{equation} \bm{w} = (\bm{X^T}\bm{X})^{-1}\bm{X^T}\bm{y} \end{equation}

が重回帰分析の係数ベクトルとなります。

ここまでで、重回帰分析のパラメータの導出まで扱いました。ここからはPythonと実際のデータセットを用いて、実際に重回帰分析をやっていきます。