1. この記事は何？
機械学習や統計学でよく扱われる線形回帰において、なぜ「二乗誤差」を最小化するのか、という疑問に対し確率的な視点から考えると納得感が出るというお話について考えてみます。
2.想定読者
そもそも最小二乗法って？という方もいらっしゃるかもしれません。すみませんが、今回の記事は回帰分析とはどういうもので最小二乗法がどういうものか、ということについてはある程度理解されている方を読者として想定しています。興味のある方は下記の記事などがわかりやすいと思います。
3. 何が問題なのか？
-
線形回帰分析では、データの散布状態への「当てはまりの良い直線」を探す際に、各点から回帰直線までの垂直距離の二乗和を最小化します。いわゆる最小二乗法です。
-
しかし、だいたいの統計学入門書に書いてあることは「差異をそのまま合計するとゼロになってしまうから二乗する」という話までです。それだけの理由だったら必ずしも二乗でなくても良いはずです。なぜ絶対値や他のべき乗ではなく二乗なのか？と思ったことはないでしょうか。この点について、書籍にきちんと答えが書いてあることは少ないです。記載があっても「絶対値を使うより計算がしやすいから」というような雑な説明で終わっていたりします。
-
もちろん、中には説明してくれている親切な書籍もあります。たとえば杉山聡『本質を捉えたデータ分析のための分析モデル入門』では計算しやすいという理由の他にもいくつかの理由を紹介してくれています。
- 最良線形不偏推定量(BLUE)だから
- 最尤推定、正規分布と相性がよいから
- 最小二乗法を情報幾何の観点から考えるとベクトルから平面に垂直に下ろした距離を最小化することと等価だから
-
…というような説明ですが、どの説明もレベルが高すぎてちょっと何言ってるかわからないので困ってしまいます。
-
もう少し分かる説明ないかな…ということで、いろいろ調べた中で、とても納得感があったのがArtem Kirsanovさんという方が作られた動画です。
- この動画の説明に感動したので、忘れないうちに自分なりに理解したことをまとめておこうと思いました。それが本記事となります。
4. 確率的な視点の導入
- データは「理想的な線形モデル + ノイズ」によって生成されると考える。
- 「理想的な線形モデル」を式にすると $ y_{ideal} = W^TX $ となります。この式では右辺のウェイト$W$、 入力値$X$をベクトルとして考えてまとめて表現しています。 $W^T$ は $X$ の各値にかけるウェイトである$W$を転置したものです。$y_{ideal}$は内積の結果なのでスカラー値です。
- 実際に観測される値 $y$ は理想的な線形モデルで予想される値 $y_{ideal}$から上下にぶれることになります。このノイズを ε とします。理想の式にノイズとして ε を加えると、観測される値 $y$ = モデル理想値 $(W^T X)$ + ノイズ $(ε)$ となります。
- このノイズ (ε) は、測定誤差や考慮外の要因など、多数の小さな独立した要因の合計であると仮定できます。
- 中心極限定理により、多数の独立な確率変数の和はガウス分布（正規分布）に従うため、ここでもノイズ $ε$ はガウス分布に従うと仮定できます。
ガウス分布（正規分布）はよく見かけるこれ↓です。σは標準偏差です。
5. 最小二乗法導出の流れはこうなる
これから行う説明の概要を先に示します。
- あるモデル（重み$W$）が与えられたとき、特定のデータ点 $(X, y)$ が観測される確率は、ちょうど $ε = y - W^TX$ となるようなノイズがガウス分布から生じる確率に等しい。
- するとデータセット全体の観測確率は、各データ点ごとの生じる確率の積で表される。
- 最適なモデル $W$ を見つけることは、「理想的な線形モデル」をめがけてこの「データセット全体での観測確率」を最大化すること（最尤推定）と同じだと言える。
- ここで計算を容易にするための便法として確率の対数を取る（対数尤度）。対数を取っても最大値を与える $W$ は変わらないはずなので計算が楽な方が良い。
- 対数尤度を最大化する計算を進めると、最終的に「二乗誤差の合計を最小化」するという、よく知られた最小二乗法の目的関数が現れる。
- つまり最小二乗法は、ノイズがガウス分布に従うという仮定のもとで、データの尤度を最大化した結果として自然に導かれる。
うむ、なるほどわからん、と思うかもしれません。概要を説明しただけなので無理もありません。そこでこれがどういうことなのかもう少し具体的に見て行きましょう。
6. 最小二乗法の導出についてさらにくわしく
観測データ $(X, y)$ がどのように生成されるかを考えます。
6-1. 線形関係
データの背後には、入力 $X$ と出力 $y$ の間に理想的に線形な関係があると仮定します。これは、ある「真の」係数＝重みベクトル $W$ を使って理想的y: $y_{ideal}= W^TX$ と表せるのでした。（$W^TX$ は $W$ と $X$ の内積）
6-2. ガウスノイズ
しかし、現実のデータ $y$ は、この理想的な値 $y_{ideal}$ に、様々な要因（測定誤差、モデル化できない要素など）によるノイズ ε が加わったものとして観測されます。
\displaylines{
y =y_{ideal} + ε \\
y_{ideal}= W^TX \\
y = y_{ideal} + ε = W^TX + ε\\
↓\\
y = W^TX + ε
}
6-3. ノイズの分布
このノイズ ε が平均 0、分散 σ² のガウス分布（正規分布）に従うと仮定します。これは、ノイズが多数の独立した小さな要因の合計である場合に中心極限定理から妥当とされる仮定です。ガウス分布の確率密度関数は有名ですが、次のようになります。
\displaylines{
σは標準偏差、e^x = exp(x)として\\
P(ε) = {\frac{1}{σ\sqrt{2π}}}\times{exp(-\frac{ε^2}{2σ^2})}
}
6-4. 尤度の最大化
今、私たちの目標は、手元にあるデータセット $D = {(X_1, y_1), ..., (X_n, y_n)}$ （n個のデータ点）を最もよく説明する（データにフィットする）モデルのパラメータ $W$ を見つけることです。いわゆる最尤推定です。
(1) 単一データ点の尤度
あるモデル係数 $W$ が与えられたとき、特定のデータ点 $(X_i, y_i)$ が観測される確率（尤度）を考えます。$y_i = W^TX_i + ε_i$ より、このデータ点が観測されるためには、ノイズが $ε_i = y_i - W^TX_i$ という値を取る必要があります。
この $ε_i$ がガウス分布に従うと仮定しているため、特定のノイズ値 $ε_i=y_i−W^TX_i$ が発生する確率密度はガウス分布の式で表現できるはずです。このことを用いて入力値$X_i$、モデル$W$、分散$σ^2$ が与えられたときに$y_i$ が観測される条件付き確率（尤度）を表すと、次のようになります。
\displaylines{
P(y_i | X_i, W, σ^2) = P(ε_i = y_i - W^TX_i) \\
= (\frac{1}{σ\sqrt{2π}}) \times exp(-\frac{(y_i - W^TX_i)² }{2σ²}) ←ガウス分布確率密度の式より
}
これは、モデル $W$ のもとでデータ点 $(Xᵢ, yᵢ)$ が観測される「もっともらしさ」を表します。
(2) データセット全体の尤度
通常、データセットの各点は互いに独立に生成されたと仮定します。そのため、モデル$W$、分散$σ^2$ が与えられたときにデータセット全体 $D$ が観測される確率（同時尤度）は、各データ点の尤度の積で表されます。
\displaylines{
P(D | W, σ²) = \prod_{i=1}^nP(yᵢ | Xᵢ, W, σ²) \\
上記(1)の式より\\
\prod_{i=1}^nP(yᵢ | Xᵢ, W, σ²) = \prod_{i=1}^n[(\frac{1}{σ\sqrt{2π}}) \times exp(-\frac{(y_i - W^TX_i)² }{2σ²})]\\
すなわち\\
P(D | W, σ²) = \prod_{i=1}^n[(\frac{1}{σ\sqrt{2π}}) \times exp(-\frac{(y_i - W^TX_i)² }{2σ²})]
}
ここで $\prod_{i=1}^n$ は 1 から n までのすべての i についての積を表します。
(3) 対数尤度の導入
この尤度 $P(D | W, σ²)$ を最大化する $W$ を見つけたいのですが積の形は扱いにくいという問題があります。そこで対数を取ります。式はこうなります。
log(P(D|W,σ²)) = log[\prod_{i=1}^n[(\frac{1}{σ\sqrt{2π}}) \times exp(-\frac{(y_i - W^TX_i)² }{2σ²})]
対数関数は単調増加なので尤度 $P(D|W, σ²)$ を最大化する$W$は、対数尤度 $log(P(D|W, σ²))$ を最大化する $W$ と同じになるはずです。
いったん対数にしてしまえば、学校で習った以下のような対数の性質を使って式を簡単にできます。
\displaylines{
log(ab) = log(a) + log(b) ←中の掛け算はlogの足し算になる\\
log(e^x) = log_ee^x = log_eexp(x)= xlog_ee = x ←logとexpは逆関数の関係にある\\
log\frac1a = log(a^{-1})=-log(a)←分数は-1乗の指数になり、指数はlogの外に出せる
}
これらを利用すると、対数尤度の式の変形はこうなります。
まずlogの中の$\prod$はlogの足し算にできるので$\sum$に書き換えられます。
log(P(D|W,σ²)) = \sum_{i=1}^n log [(\frac{1}{σ\sqrt{2π}}) \times exp(-\frac{(y_i - W^TX_i)² }{2σ²})]
logの中の掛け算をlogの足し算に変形します。
log(P(D|W,σ²)) = \sum_{i=1}^n [log (\frac{1}{σ\sqrt{2π}}) + log(exp(-\frac{(y_i - W^TX_i)² }{2σ²}))]
さらに分数を-1乗に変形したうえで指数をlogの外に出し、log(exp(...)) の部分は打ち消し合って中の式だけが残ります
log(P(D|W,σ²)) = \sum_{i=1}^n [-log{(σ\sqrt{2π})} -\frac{(y_i - W^TX_i)² }{2σ²}]
$\sum$を分割します。
log(P(D|W,σ²)) = \sum_{i=1}^n [-log(σ\sqrt{2π})] - \sum_{i=1}^n[\frac{(y_i - W^TX_i)² }{2σ²}]
右辺の最初の項の$Σ$の中には i がいないので$\sum$は単純な$n$倍になり、第2項で $\frac1{2σ^2}$は定数なので$\sum$の外に出せます。
log P(D | W, σ²) = -n \times log(σ\sqrt{2π}) - \frac1{2σ²} \sum_{i=1}^n(y_i - W^TX_i)²
(4) 最小二乗法への接続
ところで、私たちはこの対数尤度 $log(P(D | W, σ²))$ を 最大化 する $W$ を探していたのでした。最後にできあがった式をもう一度見てみます。
log(P(D | W, σ²)) = -n \times log(σ\sqrt{2π}) - \frac1{2σ²} \sum_{i=1}^n(y_i - W^TX_i)^2
この式は以下のように評価できます。
log(P(D | W, σ²)) = (Wに依存しない定数項) - (正の定数) \times \sum_{i=1}^n(y_i - W^TX_i)^2
- 第1項 $-n \times log(σ\sqrt{2π})$ は $W$ に依存しない定数なので、$W$ の最適化には影響しません。
- 第2項の係数 $- \frac1{2σ^2} $ は負の定数です。
こう考えると、対数尤度$log(P(D | W, σ²))$を最大化することはこの式の中の $\sum_{i=1}^n(y_i - W^TX_i)^2$ の部分を 最小化 することと等価であるとわかります。
したがって、対数尤度を最大化する問題は結局のところ、以下の式を最小化するという問題に帰着します。
\sum_{i=1}^n(y_i - W^TX_i)^2
7.結論
ここまでで明らかなように、$\sum_{i=1}^n(y_i - W^TX_i)^2 $ は、モデルの理想値 $y_{ideal(i)}=W^TX_i$ と実際の観測値 $y_i$ との誤差 $ε_i$の二乗和に他なりません。
\sum_{i=1}^n(y_i - W^TX_i)^2 = \sum_{i=1}^n(y_i - y_{ideal(i)})^2 = \sum_{i=1}^nε_i^2
本記事のポイントは、この最小二乗法というアプローチが単に「計算がしやすいから」などという理由ではなく、「観測データのノイズはガウス分布に従う」という確率的な仮定のもとでデータの尤度を最大化する（＝最もデータをもっともらしく説明するモデルを選ぶ）という、より根本的な原理から自然に導かれるということです。Kirsanov氏の動画で説明されていたのはこの導出プロセスでした。また杉山聡氏のいう「最尤推定、正規分布と相性がよい」というのもこのことを指しているのだと思われます。動画ではさらにこの後正則化の話題にも触れていましたが、長くなりましたので割愛しますこのテーマについては別途、続編となる記事を作成しました。よろしければご笑覧ください。
以上です。ここまで読んでいただきありがとうございました。
この記事がどなたかのお役に立てば幸いです。