線形回帰をベイズ統計学の範疇で説明します。積分も微分も必要なくて、行列の計算だけで面白い結果を得る事が出来ます。
ベイズの定理
ベイズの定理の復習をします。詳しい話は以下の記事をご覧ください。
確率密度関数や確率質量関数の事を確率と呼んでしまいます。
確率\(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} \)に関する項をもう少し綺麗な形にしましょう。その為には逆行列に関する 公式が必要 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} \)の分散とした 新しいデータのノルム分だけ増加させます。つまり、今までのデータと比べて大きく異なっていれば、分散を大きくするという事です。
言葉で書くと納得のいく式ですが、計算は大変でした。機械学習は計算が全部コンピュータなので、どうしてその結果が出てきたか良く分からないということがあります。一方、ベイズ統計学は結果の式を見ると何が起きているか分かることが多いです。
まとめ
・線形回帰をベイズ統計学の立場で書き直すことが出来る。
・事前分布を正規分布に取ることで、計算が出来る。
・ベイズ統計学で話を進めると可読性の高い結果を得る事が出来る。
・近いうちに数値実験を追記します。