ベイズ線形回帰 | マサムネの部屋

線形回帰をベイズ統計学の範疇で説明します。積分も微分も必要なくて、行列の計算だけで面白い結果を得る事が出来ます。

ベイズの定理
線形回帰と確率の関係付け
ベイズの定理による計算
行列の計算
まとめ

ベイズの定理

ベイズの定理の復習をします。詳しい話は以下の記事をご覧ください。

、ベイズの定理について解説し、具体例の計算をします。良くあるサイコロの問題と、正規分布を事前分布に仮定した時の事後分布のパラメーターを求めます。

確率密度関数や確率質量関数の事を確率と呼んでしまいます。
確率$p(x,y)$に対して、$p(x), p(x|y)$は以下で定義されました。
$$ \begin{eqnarray}
p(x) &=& \int_x p(x,y) dy \\
p(x|y) &=& \frac{p(x,y)}{p(y) }
\end{eqnarray}$$
ベイズの定理は次の式です。
$$ \begin{eqnarray}
p(x|y) &=& \frac{p(y|x)p(x )}{p(y) }
\end{eqnarray}$$

線形回帰と確率の関係付け

ベイズの定理で計算するために、正規分布と多次元正規分布を使います。基本的な性質は以下の記事にまとめてあります。

正規分布の性質

正規分布の導入をします。確率になっていること、平均の値、分散の値を具体的に計算します。

多次元正規分布の性質

多次元正規分布の定義の平均値や分散を計算します。最尤推定によって得られた平均値や分散が、不偏推定量になっているか確かめます。その結果をもとに、多次元正規分布が正規分布の拡張になっている事を確かめます。

線形回帰は、以下のモデルの誤差関数$ \epsilon $ を最小にするモデルでした。
$$\begin{eqnarray}
y=\vec{w} \cdot \vec{x}+ \epsilon
\end{eqnarray}$$
確率とつなげる為に、$ \epsilon = \mathcal{N} ( \epsilon | 0, \lambda ^{-1} ) $ を仮定しましょう。ただし、平均$ \mu $ 、分散$ \lambda ^{-1} $の正規分布を$ \mathcal{N} ( x|\mu , \lambda ^{-1} ) $で表しました。すると、 yは$ \vec{x} , \vec{w} $を与えたとき、以下の確率に従います。
$$\begin{eqnarray}
y=\mathcal{N}(y |\vec{w} \cdot \vec{x} , \lambda ^{-1} )
\end{eqnarray}$$
目標は、データxを与えたときの予測値yの確率$ p(y|x ) $ です。計算を進めるため、$ \vec{w} $が従う確率分布を多次元正規分布に仮定します。
$$\begin{eqnarray}
p( \vec{w} )=\mathcal{N}(\vec{w} |\vec{m}, \Lambda ^{-1} )
\end{eqnarray}$$
$ \vec{m}, \Lambda ^{-1} $はハイパーパラメーターで、適当なベクトルや行列を入れます。ゼロベクトルと単位行列とかで良いです。
これで最低限の準備が整いました。計算していきましょう。

ベイズの定理による計算

データ$ X = \{ \vec{x_1} , \cdots , \vec{x_n} \} , Y=\{ y_1 , \cdots , y_n \} $が得られたとして、$ p(\vec{w} |X,\vec{y} ) $を求めます。
$$\begin{eqnarray}
p(\vec{w} |X,\vec{y} ) =\frac{ p(\vec{w} ) \prod p(y_i |x_i , \vec{w} )} {p(Y|X) }
\end{eqnarray}$$
$ \vec{w} $に関係があるのは分子のみです。定番ですが、$ \log $をとって $ \vec{w} $ で整理しましょう。
$$\begin{eqnarray}
& \log& p(\vec{w} |X,\vec{y} ) \\
&=&-\frac{1}{2} \left\{ \vec{w}^{T}
(\lambda \sum \vec{x_i} \vec{x_i} ^{T} +\Lambda )
-2\vec{w} ^{T} (\lambda \sum y_i \vec{x_i} +\Lambda \vec{m} \right\} +const
\end{eqnarray}$$
データX,Yからパラメーター$ \vec{w}$が以下のように予測できます。

$$\begin{eqnarray}
p(\vec{w} |X,Y ) &=& \mathcal{N} (\vec{w} | \hat{ \vec{m} }, \hat{\Lambda } ^{-1}) \\
\hat{\Lambda } &=& \lambda \sum \vec{x_i} \vec{x_i} ^{T} + \Lambda \\
\hat{\vec{m} } &=& \hat{\Lambda} ^{-1} \left( \lambda \sum y_i \vec{x_i} +\Lambda \vec{m} \right)
\end{eqnarray}$$

新しくデータ$ x_{\ast} $を得たとして、その予測値を$ y_{\ast} $としましょう。ベイズの定理から
$$\begin{eqnarray}
p(\vec{w} |x_{\ast} ,y_{\ast} ) = \frac{ p(\vec{w} ) p(y_{\ast}| \vec{x} , \vec{w} )}{p(y_{\ast}|x_{\ast} ) }
\end{eqnarray}$$
です。これを変形して、
$$\begin{eqnarray}
p(y_{\ast}| x_{\ast} ) \propto \frac{ p(y_{\ast}| \vec{x} , \vec{w} )}{ p(\vec{w} |x_{\ast} ,y_{\ast} ) }
\end{eqnarray}$$
logを取って、$ y_{\ast} $についてまとめましょう。正規分布が出てくるはずです。
$$\begin{eqnarray}
\log p(y_{\ast}| x_{\ast} ) \propto
-\frac{1}{2} \left\{
(\lambda – \lambda ^2 \vec{x_{\ast} }^{T} \hat{\Lambda }^{-1} \vec{x_{\ast}} )y_{\ast} ^2
-2 ( \lambda \vec{x_{\ast}} ^{T} \hat{\Lambda }^{-1} \Lambda \vec{m} ) y_{\ast}
\right\}
\end{eqnarray}$$
この式から、$ y_{\ast} $の従う確率は以下になります。

$$\begin{eqnarray}
p(y_{\ast}| x_{\ast} ) &=& \mathcal{N}\left( y_{\ast} |\vec{m} _{\ast} , \lambda _{\ast} \right) \\
\lambda _{\ast} &=& \lambda – \lambda ^2 \vec{x_{\ast} }^{T} \hat{\Lambda }^{-1} \vec{x_{\ast}} \\
\vec{m} _{\ast} &=& \lambda _{\ast} ^{-1} ( \lambda \vec{x_{\ast}} ^{T} \hat{\Lambda }^{-1} \Lambda \vec{m} )
\end{eqnarray}$$

この段階でも実装しようと思えば出来ますが、式が非常に汚いです。$ \hat{\Lambda }^{-1} $に関する項をもう少し綺麗な形にしましょう。その為には逆行列に関する公式が必要 ¹ です。

行列の計算

良く出てくる行列に関しては、逆行列が知られています。

[ウッドベリーの公式]
$$\begin{eqnarray}
(A+UBV)^{-1} = A^{-1} – A^{-1} U( B^{-1} +VA^{-1} U )VA^{-1}
\end{eqnarray}$$
証明は、$ A+UBV $にかけてみると単位行列になることを確認すればよいです。

回帰分析の話に戻って、$x_{\ast} ^{T} \hat{\Lambda }^{-1} x_{\ast} $を計算しましょう。初めに、$ \hat{\Lambda } = \Lambda + \vec{x_{\ast}} \lambda \vec{x_{\ast} } ^{T} $にウッドベリーの公式を当てはめると、逆行列が計算できます。
$$\begin{eqnarray}
( \Lambda + \vec{x_{\ast}} \lambda \vec{x_{\ast} } ^{T} )^{-1} &=& \Lambda^{-1} – \Lambda ^{-1} \vec{x_{\ast}} ^{T} ( \lambda ^{-1} + \vec{x_{\ast} } ^{T} \Lambda ^{-1} \vec{ x_{\ast} } ) \vec{x_{\ast} } ^{T} \Lambda ^{-1} \\
&=& \Lambda^{-1}
– \frac{1}{ \lambda ^{-1} + \vec{x_{\ast} } ^{T} \Lambda ^{-1} \vec{ x_{\ast} } }
\Lambda ^{-1} \vec{x_{\ast}} \vec{x_{\ast} } ^{T} \Lambda ^{-1}
\end{eqnarray}$$
計算の見通しを良くするために、$ \mathcal{L} = \vec{x_{\ast} } ^{T} \Lambda ^{-1} \vec{ x_{\ast} } $と置きましょう。この置き換えをしておいて、次のような計算が出来ます。
$$\begin{eqnarray}
x_{\ast} ^{T} \hat{\Lambda }^{-1} x_{\ast} &=&
\mathcal{L} – \frac{\mathcal{L} ^2} { \lambda ^{-1} + \mathcal{L} } \\
&=& \frac{ \lambda ^{-1} \mathcal{L} }{\lambda ^{-1} + \mathcal{L} }
\end{eqnarray}$$
この計算結果で、汚かった$y_{\ast} $を計算してみましょう。

$$\begin{eqnarray}
\lambda _{\ast} &=& \lambda – \frac{ \lambda \mathcal{L} }{\lambda ^{-1} + \mathcal{L} } \\
&=& \frac{1}{\lambda ^{-1} +\mathcal{L} } \\

\vec{m} _{\ast}
&=&
(\lambda ^{-1} +\mathcal{L} ) ( \lambda \vec{x_{\ast}} ^{T} \hat{\Lambda }^{-1} \Lambda \vec{m} ) \\
&=&
\lambda (\lambda ^{-1} +\mathcal{L} ) \left( 1 – \frac{\mathcal{L}}{\lambda ^{-1} + \mathcal{L} } \right) \vec{x_{\ast}} ^{T} \vec{m} \\
&=&
\vec{x_{\ast}} ^{T} \vec{m}
\end{eqnarray}$$
この計算から、予測式は次のようにまとめる事が出来ます。
$$\begin{eqnarray}
p(y_{\ast}|\vec{x_{\ast}} )= \mathcal{N} \left( y_{\ast} | \vec{x_{\ast}} ^{T} \vec{m} , \lambda ^{-1} +\vec{x_{\ast}} ^{T} \Lambda ^{-1} \vec{x_{\ast}} \right)
\end{eqnarray}$$
事前にデータを取っていて、パラメーターが推定できている時は、$ \vec{m} ,\Lambda $ の代わりに $ \hat{\vec{m}} ,\hat{\Lambda} $を使います。

予測式の中身を見てみましょう。
$ \lambda , \vec{m} , \Lambda $ ハイパーパラメーターです。$ \lambda $ は誤差関数の分散なので、データ全体の分散を表すハイパーパラメーターです。$ \vec{m} , \Lambda $ は、パラメーター $ \vec{w} $の平均値と分散でした。予測値の平均値は、大体$ \vec{w} ^{T} \vec{x_{\ast}} $になると言っています。分散はについては誤差関数の分はそのまま受け継ぎ、計量を$ \vec{w} $の分散とした新しいデータのノルム分だけ増加させます。つまり、今までのデータと比べて大きく異なっていれば、分散を大きくするという事です。
言葉で書くと納得のいく式ですが、計算は大変でした。機械学習は計算が全部コンピュータなので、どうしてその結果が出てきたか良く分からないということがあります。一方、ベイズ統計学は結果の式を見ると何が起きているか分かることが多いです。