【回帰分析の数理】#2 線形回帰のパラメタを決定する正規方程式

はじめに

データマイニング・機械学習における最もシンプルな分析手法である線形回帰について説明します。前回の記事では、線形回帰のうち説明変数が1つの場合である単回帰分析を例に、尤度関数について説明しました。

【回帰分析の数理】#1 線形回帰で学ぶ最尤推定
はじめに データマイニング・機械学習における最もシンプルな分析手法である線形回帰について説明します。この記事では線形回帰を例に、データマイニングの基礎概念である尤度関数について学びます。データマイニングとは、データから意思決定を高度化...

今回の記事では、重回帰分析や多項式によるフィッティングに拡張した議論を行います。前回と同じく、データセットはこちらを使用します。

Advertising Dataset
Kaggle is the world’s largest data science community with powerful tools and resources to help you achieve your data science goals.

残差プロット

前回はデータセットのうち、売り上げをTV広告費という1変数のみを用いてモデリングすることを試みました。

グラフを観察すると、売り上げである”Sales”の値が高いほど、赤い直線からのズレが大きくなっているように見えます。次のような残差プロット(Residual Plot)を描画して、残差のパターンを確認してみます。

横軸のPredicted Valuesがモデルにより推定されたSalesの値、縦軸が残差\(\epsilon^{(i)}\)に相当します。Predicted Valuesが大きいほど、残差の絶対値が大きくなっていくパターンが見えます。

しかし、残差は\(\epsilon^{(i)} \sim \mathbb{N}(0, \sigma^2)\)とモデリングしているので、上記のようなパターンが見られるということは、TV広告費によるモデリングだけでは表現しきれないパターンがあったことを示唆しています。

重回帰分析

実は、データセットにはTV広告費以外にも、Radio広告費・Newspaper広告費という変数が格納されています。ここで、データに対する変数の命名について次のように定めておきましょう。

interceptは、前回の記事で説明した\(\theta^{\mathrm{T}}x\)という表記を可能にするための切片です。続いて、TV, Radio, Newspaperという3つの説明変数が格納されています。この3つの説明変数を用いて、目的変数であるSalesの売り上げをモデリングすることを目指します。一般的な線形の重回帰分析では、次の式のようにモデリングができます。

$$
\large
\begin{align}
y^{(i)} &=  h_\theta(x^{(i)} ) + \epsilon^{(i)} \\
&= \theta_0 x_0^{(i)} + \theta_1 x_1^{(i)} + \theta_2 x_2^{(i)} + \theta_3 x_3^{(i)} + \epsilon^{(i)}\\
&= \theta^\mathrm{T} x^{(i)}  + \epsilon^{(i)}
\end{align}
$$

すなわち、式の形としては前回の記事と全く同じように扱えます。尤度関数の最大化は、同じように最小二乗法と等しくなります。また、最小化するコスト関数である残差の二乗和を\(J(\theta)\)とおきます。

$$
\large
\begin{align}
\hat{\theta} &= \arg\max_\theta l(\theta) \\
&= \arg\min_\theta \sum_{i=1}^N (y^{(i)}-\theta^\mathrm{T} x^{(i)})^2 \\
J(\theta) &= \sum_{i=1}^N (y^{(i)}-\theta^\mathrm{T} x^{(i)})^2
\end{align}
$$

補足
偏微分する際の利便性から、残差の二乗和の二分の一をコスト関数とする流儀もあります。

正規方程式によるモデルパラメタの決定

さて、上記のような\(\theta\)は、デザイン行列(design matrix)\(X\)と目的変数のベクトル\(\vec{y}\)を用いて、次のように解析的に求めることができます(正規方程式)。デザイン行列を用いると二乗誤差の合計を簡潔に書き下すことができます。

$$
\large
\begin{align}
\nabla_\theta J(\theta) &= \nabla_\theta \sum_{i=1}^N (y^{(i)}-\theta^\mathrm{T} x^{(i)})^2 \\
&= \nabla_\theta (X\theta-\vec{y})^\mathrm{T}(X\theta-\vec{y}) \\
&= \nabla_\theta (\theta^\mathrm{T}(X^\mathrm{T}X)\theta – 2(X^\mathrm{T}\vec{y})^\mathrm{T}\theta) \\
&= 2(X^\mathrm{T}X\theta – X^\mathrm{T}\vec{y})
\end{align}
$$

補足
上記の式変形におけるベクトル微分は、以下の二つの性質を利用すると簡潔に行うことができます。ただし、\(x, b\)はベクトル、\(A\)は対称な行列、すなわち\(A^\mathrm{T}=A\)を満たします。

$$
\large
\begin{align}
\nabla_x b^\mathrm{T}x &= b \\
\nabla_x x^\mathrm{T} Ax &= 2Ax
\end{align}
$$

特に、2式目は線形代数で二次形式と呼ばれるものです。
\((X^\mathrm{T}X)^\mathrm{T} = X^\mathrm{T}(X^\mathrm{T})^\mathrm{T} = X^\mathrm{T}X\)ゆえ、\(A=X^\mathrm{T}X\)、\(x=\theta\)として対称行列であることを利用します。

\(J(\theta)\)が下向きの凸関数(convex function)であることを考慮すれば、偏微分値がゼロになる点が最小値をとります。よって次が成り立ちます。

$$
\large
X^\mathrm{T}X\theta = X^\mathrm{T}\vec{y}
$$

よってモデルパラメタは次のように決定されます。

$$
\large
\theta = (X^\mathrm{T}X)^{-1}X^\mathrm{T}\vec{y}
$$

補足
なお、これは\(X^\mathrm{T}X\)が逆行列を持つ(正則である)ことを前提にしています。ごく稀に正則でない場合は、擬似逆行列を用いるというテクニックがあります。実際にPython等を用いてこのようなアルゴリズムを実装する際に「逆行列が存在しない」というエラーが起きることがありますが、そのほとんどはコーディングのミスによるゼロ除算(Zero Division Error)やオーバーフロー(Overflow)によるものです。

結果の確認

上記のようにモデルパラメタ\(\theta\)を決定して、残差プロットを再び描画してみます。

左側が一変数のみを用いた場合、右側が三変数を用いた場合の残差プロットです。右側では、Predicted Valuesが大きくなるほど残差の絶対値が大きくなる、というパターンを緩和できていることが分かります。

まとめ

今回は、重回帰分析への線形回帰の拡張を行いました。まだまだ分析の余地が残るデータであり、ここからの変数選択や特徴量の作成などがデータマイニングの力量が問われる部分です。次回は多項式など一般の線形回帰への拡張を行います。

【回帰分析の数理】#3 線形回帰の多項式回帰への拡張
はじめに データマイニング・機械学習における最もシンプルな分析手法である線形回帰について説明します。前回の記事では、説明変数が一つである単回帰分析から、説明変数が複数である重回帰分析への拡張を行いました。 今回の記事では...

参考文献

1. Stanford cs229 Lecture Notes Linear Regression (Andrew Ng)
2. Introduction to Statistical Learning (Gareth et al.)

 

SD/T
Written by
SD/T