線形回帰における最尤推定・MAP推定・ベイズ推定

線形回帰を確率的に捉える。

はじめに

線形回帰では、 $m$ 次元特徴ベクトル $\bm x\in\R^m$ とパラメータ $\bm\beta\in\R^m$ の内積によってターゲット $y$ を予測する。

\begin{align} \hat y = \bm x^\top\bm\beta \end{align}

特徴ベクトルと対応するターゲットが $n$ 組得られたとする。これらを

$X\in\R^{n\times m}$
$\bm y\in\R^n$

と表しておく。

線形回帰では、予測値と正解の二乗和誤差が最小になるようなパラメータ $\bm\beta$ を求める。これは目的関数を微分して $=0$ とおくことで解析的に求められる。

\begin{align} J(\bm w) &= \| \bm y - X\bm\beta \|^2 \\ \frac{\partial J}{\partial \bm\beta} &= 2X^\top X\bm\beta - 2X^\top\bm y \\ \hat{\bm\beta} &= (X^\top X)^{-1}X^\top\bm y \end{align}

さてこの最小二乗法によって求めたパラメータは、正解と予測値の差が最小になるようなパラメータではなく、尤度が最大になるようなパラメータとも解釈できる。誤差が正規分布に従うと仮定した上で最尤推定を行うと、途中で二乗和誤差の最小化が出てくるためだ。

このように、線形回帰は、最尤推定をはじめとした確率分布のパラメータ推定手法から考えることが出来る。本稿では、いくつかのパラメータ推定手法から線形回帰を捉えてみる。

最尤推定に基づく線形回帰

最尤推定

尤度が最大になるパラメータ $\theta$ の推定。

\begin{align} \hat\theta = \argmax_\theta p(X;\theta) \end{align}

多くの場合で計算しやすい対数尤度 $\log p(X;\theta)$ を扱う。

線形回帰モデルのパラメータ $\bm\beta$ に対して最尤推定を行なってみよう。次の条件付き分布を考える。

\begin{align} p(y|\bm x;\bm\beta) \end{align}

ここで、ターゲット $y$ が以下のように発生すると仮定する。

\begin{align} y &= \bm x^\top\bm\beta + \epsilon \\ \epsilon &\sim \mathcal N(0, \sigma^2) \end{align}

$\epsilon$ は誤差で、平均0、分散 $\sigma^2$ の正規分布に従うと仮定した。すると先の条件付き分布を次のようにモデル化できる。

\begin{align} p(y|\bm x; \bm\beta) = \mathcal N(y; \bm x^\top\bm\beta,\sigma^2) \end{align}

この下で対数尤度を計算してみると

\begin{align} \ln p(\bm y| X; \bm\beta) &= \sum_i \ln p(y_i| \bm x_i; \bm\beta) \\ &= \sum_i \ln \mathcal N(y_i; \bm x_i^\top\bm\beta, \sigma^2) \\ &= \sum_i \ln \left( \frac{1}{\sqrt{2\pi\sigma^2}} \exp \left( -\frac{(y_i - \bm x_i^\top \bm\beta)^2}{2\sigma^2} \right) \right) \\ &= -\frac{1}{2\sigma^2} \sum_i (y_i - \bm x_i^\top \bm\beta)^2 + \text{const} \\ \end{align}

負の二乗和誤差が出てきた。最尤推定ではこれを最大化する、つまり二乗和誤差を最小化する。ここからの流れは同じ。線形回帰における最小二乗法によるパラメータ推定は、誤差に正規分布を仮定した最尤推定であると言える。

MAP推定に基づく線形回帰

MAP推定

事後確率が最大になるパラメータ $\theta$ の推定。

\begin{align} \hat\theta &= \argmax_\theta p(\theta|X) \\ &= \argmax_\theta p(X;\theta)p(\theta) \end{align}

最尤推定ではあるパラメータ $\theta^*$ における尤度 $p(X;\theta^*)$ のみを考慮していたが、MAP推定では、「そもそもそのパラメータが得られる確率はどんなもんなん？」というところまで考慮する。主観に基づくパラメータそのものの妥当性 $p(\theta^*)$ で尤度 $p(X;\theta^*)$ に重み付けする。

MAP推定でも同じように線形回帰を考えられる。

\begin{align} \hat{\bm\beta} &= \argmax_{\bm\beta} p(\bm\beta| X, \bm y) \\ &= \argmax_{\bm\beta} p(\bm y| X; \bm\beta)p(\bm\beta) \end{align}

尤度 $p(\bm y| X; \bm\beta)$ は最尤推定の章で示した通り、誤差に正規分布を仮定してモデル化する。事前分布 $p(\bm\beta)$ は自身で適当なものを設定する。

例として、事前分布に平均 $0$ 、分散 $1/\lambda$ をもつ独立な正規分布を仮定してみよう。

\begin{align} p(\bm\beta) = \mathcal N(\bm\beta; \bm 0, \frac{1}{\lambda} I) \end{align}

対数をとる。

\begin{align} \ln p(\bm\beta) &= \ln \mathcal N(\bm\beta; \bm 0, \frac{1}{\lambda} I) \\ &= \ln \left( \left( \frac{\lambda}{2\pi} \right)^{d/2} \exp\left( -\frac{1}{2}\lambda \| \bm\beta \|^2 \right) \right) \\ &= -\frac{1}{2}\lambda \|\bm\beta\|^2 + \text{const} \end{align}

材料が揃ったので、対数事後分布を見てみる。

\begin{align} \ln p(\bm\beta| X, \bm y) &= \ln p(\bm y| X; \bm\beta) + \ln p(\bm\beta) \\ &= -\frac{1}{2} \| \bm y - X\bm\beta \|^2 - \frac{1}{2}\lambda \| \bm\beta\|^2 + \text{const} \end{align}

見やすくするために条件付き分布の分散 $\sigma^2$ は1としている。

MAP推定ではこれを最大化するわけだが、さてこの式、どこかで見たことがないだろうか？実はこれはRidge回帰の目的関数と一致する。Ridge回帰は二乗和誤差にパラメータのL2ノルムを正則化項として加えた以下の目的関数を最小化する。これは上記の事後分布の最大化と同じ意味になる。

\begin{align} J(\bm\beta) = \| \bm y - X\bm\beta \|^2 + \lambda \|\bm\beta \|^2 \end{align}

つまりRidge回帰は、事前分布に平均0の独立な正規分布を仮定したMAP推定と見られる。この事前分布が正則化項として働く。この正則化項は、機械学習的には「パラメータが0から離れたときにペナルティを与えるもの」と解釈できるが、ベイズ的には「そもそもパラメータって0付近から発生するよね」という立場を反映させたものとなる。事前分布の分散 $1/\lambda$ が小さいほど= $\lambda$ が大きいほど強い正則化が行われるようになるが、これは分散が小さくなることで「パラメータって0付近から発生するよね」という立場が強くなるためである。

L2ノルムで正則化を行うRidge回帰だけでなく、L1ノルムで正則化を行うLasso回帰もMAP推定の特殊な形として解釈できる。事前分布としてラプラス分布を仮定してMAP推定を進めるとLasso回帰の目的関数が出現する。

\begin{align} p(\bm\beta) &= \prod_j \frac{1}{2b} \exp (-\lambda |\beta_j|) \\ \ln p(\bm\beta| X, \bm y) &= \ln p(\bm y| X; \bm\beta) + \ln p(\bm\beta) \\ &= -\frac{1}{2} \| \bm y - X\bm\beta \|^2 - \lambda \| \bm\beta \|_1 + \text{const} \end{align}

ベイズ推定に基づく線形回帰

ベイズ推定

事後分布そのものの推定。

\begin{align} p(\theta|X) &= \frac{p(X|\theta) p(\theta)}{p(X)} \end{align}

MAP推定ではこれが最大となる $\theta$ を求めていたが、ベイズ推定では分布そのものを求める。

線形回帰におけるパラメータ $\bm\beta$ についてベイズ推定を行ってみる。

\begin{align} p(\bm\beta| X, \bm y) &= \frac{p(\bm y| X; \bm\beta)p(\bm\beta)}{p(\bm y| X)} \end{align}

尤度 $p(y| \bm x; \bm\beta)$ はこれまで通り正規分布でモデル化する。そして事前分布 $p(\bm\beta)$ にも共役な分布である正規分布を仮定する。

\begin{align} p(y| \bm x; \bm\beta) &= \mathcal N(y; \bm x^\top\bm\beta, \sigma_y^2) \\ p(\bm\beta) &= \mathcal N(\bm\beta; \bm 0, \sigma_\beta^2 I) \end{align}

すると、事後分布を以下の正規分布で表せる。（証明は付録）

\begin{align} p(\bm\beta| X, \bm y) &= \mathcal N(\bm\beta; \bm\mu_{\bm\beta}, \Sigma_{\bm\beta}) \\ \Sigma_{\bm\beta} &= \left( \frac{1}{\sigma_y^2}X^\top X - \frac{1}{\sigma_{\bm\beta}^2}I \right)^{-1} \\ \bm\mu_{\bm\beta} &= \Sigma_{\bm\beta} \left( \frac{1}{\sigma_y^2} X^\top\bm y \right) \end{align}

これでベイズ推定が完了し、事後分布が求められた。

線形回帰におけるベイズ推定には優れた側面がある。事後分布を求める事で、新たに得られた特徴ベクトル $\bm x_*$ に対応するターゲット $y_*$ が分布として求められるようになる。

\begin{align} p(y_*| X, \bm y, \bm x_*) &= \int p(y_*| \bm x_*; \bm\beta) p(\bm\beta| X, \bm y) \, d\bm\beta \end{align}

これは先と同じ仮定でモデル化すると次の正規分布で表せる。（証明は付録）

\begin{align} p(y_*| X, \bm y, \bm x_*) &= \mathcal N(y_*; \mu_*, \sigma_*^2) \\ \mu_* &= \bm x_*^\top \bm\mu_{\bm\beta} \\ \sigma_*^2 &= \bm x_*^\top \Sigma_{\bm\beta} \bm x_* + \sigma_y^2 \end{align}

実際に求めてみよう。適当なサンプルを用意する。

Loading image...

ここに三次式を仮定してベイズ推定を行う。これまで同様、条件付き分布、事前分布には適当な正規分布を仮定した。推定した事後分布を用いて予測分布を求めると図のようになる。

Loading image...

中心の太線は期待値、色がついている範囲は標準偏差 $\sigma_*$ を表している。

図を見ると、ある部分では分散が大きく、またある部分では分散が小さいことが分かる。分散の小さい部分の予測は比較的信頼できる一方、そうでない部分の予測には不確実性が残ることになる。ベイズ推定ではこのような不確実性の比較ができるため、より多くの情報を考慮した意思決定が可能になる。

👀

点推定でも条件付き分布 $p(y| \bm x; \bm\beta)$ を定義しているので予測値が分布として求められるが、ベイズ推定と異なり、分散が目的変数 $\bm x$ に依らないため、相対的な不確実性の高さが分からない。

おわりに

線形回帰における各種推定手法をまとめた。機械学習モデルを確率モデルとして解釈するための基礎的な理解が身についたと思う。

参考文献

須山敦志. ベイズ深層学習. 講談社, 2019.
後藤正幸, 小林学. 入門パターン認識と機械学習. コロナ社, 2014.

付録

事後分布の証明

以下を示す。

\begin{align} p(\bm\beta| X, \bm y) &= \mathcal N(\bm\beta; \bm\mu_{\bm\beta}, \Sigma_{\bm\beta}) \end{align}

条件付き分布、事前分布が以下とすると

\begin{align} p(y| \bm x; \bm\beta) &= \mathcal N(y; \bm x^\top\bm\beta, \sigma_y^2) \\ p(\bm\beta) &= \mathcal N(\bm\beta; \bm 0, \sigma_\beta^2 I) \end{align}

対数事後分布はこうなる。

\begin{align} &\ln p(\bm\beta| X, \bm y) \\ &= \ln p(\bm y| X; \bm\beta) + \ln p(\bm\beta) + C_1 \\ &= -\frac{1}{2\sigma_y^2} \| \bm y - X\bm\beta \|^2 - \frac{1}{2\sigma_{\bm\beta}^2} \| \bm\beta \|^2 + C_2 \\ &= -\frac{1}{2\sigma_y^2} (\bm y - X\bm\beta)^\top(\bm y - X\bm\beta) - \frac{1}{2\sigma_{\bm\beta}^2} \bm\beta^\top \bm\beta + C_2 \\ &= -\frac{1}{2\sigma_y^2} (\bm y^\top\bm y - 2\bm y^\top X\bm\beta + \bm\beta^\top X^\top X\bm\beta) - \frac{1}{2\sigma_{\bm\beta}^2} \bm\beta^\top \bm\beta + C_2 \\ &= -\frac{1}{2\sigma_y^2}\bm\beta^\top X^\top X\bm\beta - \frac{1}{2\sigma_{\bm\beta}^2} \bm\beta^\top \bm\beta + \frac{1}{\sigma_y^2} \bm y^\top X\bm\beta + C_3 \\ &= -\frac{1}{2}\bm\beta^\top \left( \frac{1}{\sigma_y^2}X^\top X - \frac{1}{\sigma_{\bm\beta}^2}I \right) \bm\beta + \bm\beta^\top \left( \frac{1}{\sigma_y^2} X^\top\bm y \right) + C_3 \\ \end{align}

ここで、多変量正規分布は

\begin{align} \ln \mathcal N(\bm x; \bm\mu, \Sigma) &= - \frac{1}{2} (\bm x - \bm\mu)^\top \Sigma^{-1} (\bm x - \bm\mu) + C_4 \\ &= - \frac{1}{2} \bm x^\top\Sigma^{-1}\bm x - \bm x^\top\Sigma^{-1}\bm\mu + C_5 \\ \end{align}

なので、先の事後分布は以下の多変量正規分布とみなせる。

\begin{align} p(\bm\beta| X, \bm y) &= \mathcal N(\bm\beta; \bm\mu_{\bm\beta}, \Sigma_{\bm\beta}) \\ \Sigma_{\bm\beta} &= \left( \frac{1}{\sigma_y^2}X^\top X - \frac{1}{\sigma_{\bm\beta}^2}I \right)^{-1} \\ \bm\mu_{\bm\beta} &= \Sigma_{\bm\beta} \left( \frac{1}{\sigma_y^2} X^\top\bm y \right) \end{align}

予測分布の証明

以下を示す。

\begin{align} p(y_*| X, \bm y, \bm x_*) &= \mathcal N(\bm y_*; \mu_*, \sigma_*^2) \\ \mu_* &= \bm x_*^\top \bm\mu_{\bm\beta} \\ \sigma_*^2 &= \bm x_*^\top \Sigma_{\bm\beta} \bm x_* + \sigma_y^2 \end{align}

といっても、多変量正規分布の線形変換に関する定理を使うだけ。

定理

$\bm x \in \R^d$ が多変量正規分布 $\mathcal N(\bm\mu, \Sigma)$ に従うとき、 $\bm x$ の線形変換 $\bm y = A\bm x + \bm b$ は、平均 $A\bm\mu + \bm b$ 、共分散 $A\Sigma A^\top$ の多変量正規分布に従う。

\begin{align} \bm y \sim \mathcal N(A\bm\mu + \bm b, A\Sigma A^\top) \end{align}

$\bm\beta \sim \mathcal N(\bm\mu_{\bm\beta}, \Sigma_{\bm\beta})$ なので、 $A=\bm x_*^\top, \bm b=\bm 0$ とすると、定理より以下が成り立つ。

\begin{align} \bm x_*^\top\bm\beta \sim \mathcal N(\bm x_*^\top\bm\mu, \bm x_*^\top\Sigma_{\bm\beta}\bm x_*) \end{align}

そしてここにノイズ $\epsilon \sim \mathcal N(0, \sigma_y^2)$ を足した $y_* = \bm x_*^\top\bm\beta + \epsilon$ も当然正規分布となり、平均・分散をそれぞれ足して

\begin{align} y_* &\sim \mathcal N(\bm x_*^\top\bm\mu, \bm x_*^\top \Sigma_{\bm\beta} \bm x_* + \sigma_y^2) \end{align}

が成り立つ。