こんにちは、テックタヌキです。
普段はデータ分析のお仕事をしていて、このブログでは学んだことのアウトプットをしています。
本記事では、前回の記事に引き続き、ベイズ推定について学んだことをまとめます。
何か間違っているところがあれば指摘してもらえると喜びます。
前回のおさらい
前回は（分散既知の）１次元ガウス分布の平均値のパラメータ推定を行いました。
共役事前分布を与えた場合、解析的に求めることができましたが、
共役事前分布ではない標準コーシー分布で与えた場合には、周辺尤度であるp(X)=\int P(X|\mu)p(\mu) d\mu = \int_{-\infty}^{\infty} \exp(- \sum_n(x_n- \mu)^2/2)/\pi(1+\mu^2) d\muの積分が解析的に求めることができず、モンテカルロ法によって求めました。
本記事では、モンテカルロ法ではなく、変分ベイズ法によって求める方法について解説していきます。
変分ベイズ
結論から説明します。以後未知パラメータをまとめて、Zと書きます。
変分ベイズの根幹アイデアは、事後分布p(Z|X)を求めることはあきらめて、近似分布q(Z)をもとめることです。
ZがZ=(z_1, \cdots, z_K)のように複数の未知パラメータからなる場合は、パラメータ間に独立性の仮定を入れて
p(Z|X) \simeq q(Z_1)q(Z_2) \cdots q(Z_G)
とします。 ここで、Gはパラメータをグループ分けしたときのグループ数です。パラメータのグループごとに独立性を仮定します。
すべてのパラメータ間に独立性を仮定する場合は、G=Kです。
話を元に戻します。近似分布との近さを測る指標としては、以下のKLダイバージェンスを用います。
KL(q(Z)||p(Z|X)) = E_{Z \sim q(Z)} \Big[\log \frac{q(Z)}{p(Z|X)}\Big] = \int q(Z) \log \frac{q(Z)}{p(Z|X)} dZ
このKLダイバージェンスが最小化されるように、q(Z)を求めることが目標ですが、
KLダイバージェンスの中に事後分布が含まれているので、計算ができません。
しかし、実はこのKLダイバージェンスの最小化がほかの言葉で置き換えられるということを以下に示していきます。
対数周辺尤度の計算
p(X)を周辺尤度といい、p(X)=\int p(X, Z)dzで(X, Z)の同時分布から周辺化できます。
この周辺尤度に対して対数を取った\log p(Z)を対数周辺尤度といい、q(Z)を使って、イェンセンの不等式を使って数式をいじってみます。
\begin{aligned}
\log p(X) &= \log \int P(X, Z) dZ \\
&= \log \int q(Z)
\frac{P(X, Z)}{q(Z)} dZ \\
& \geq \int q(Z) \log\frac{P(X, Z)}{q(Z)} dZ =: L(q(Z))
\end{aligned}
最後の不等式が、イェンセンの不等式です。\logが上に凸で、q(Z)を足し合わせると1になることに注意すればこの演算の正当性はわかるかと思います。
この下から抑えたものを、ここではL(q(Z))と表します。
ここで、天下り的ですが、L(q(Z))とKL(q(Z)||p(Z|X))を足し合わせたものを計算してみます。
\begin{aligned}
L(q(Z)) + KL(q(Z)||p(Z|X)) &= \int q(Z) \log\frac{P(X, Z)}{q(Z)} dZ + \int q(Z) \log \frac{q(Z)}{p(Z|X)} dZ \\
&= \int q(Z) \log \Big(\frac{P(X, Z)}{q(Z)}\frac{q(Z)}{p(Z|X)}\Big) dZ\\
&= \int q(Z) \log \Big(\frac{P(X, Z)}{p(Z|X)}\Big) dZ\\
&= \int q(Z) \log \Big(\frac{P(X, Z)}{p(Z,X)/p(X)} \Big) dZ \\
&= \int q(Z) \log p(X) dZ \\
&= \log p(X) \int q(Z) dZ = \log p(X)
\end{aligned}
すなわち、L(q(Z))とKL(q(Z)||p(Z|X))を足し合わせたものは、対数周辺尤度に一致するということです。ここで、対数周辺尤度は一定なので、特に
L(q(Z))+KL(q(Z)||p(Z|X))= const.ということがいえました。
KLダイバージェンスの最小化の言いかえ
もともとの目的は、KL(q(Z)||p(Z|X))を最小化することでした。
ここで、L(q(Z))+KL(q(Z)||p(Z|X))= const.なので、KL(q(Z)||p(Z|X))を最小化するということは、L(q(Z))を最大化することに言い換えられます。
L(q(Z))の計算を進めてみます。
\begin{aligned}
L(q(Z))&:=\int q(Z) \log\frac{p(X, Z)}{q(Z)} dZ \\
&= \int q(Z) \log\frac{p(X|Z)p(Z)}{q(Z)} dZ \\
&= \int q(Z) \log p(X|Z) + q(Z) \log \frac{p(Z)}{q(Z)} dZ \\
&= \int q(Z) \log p(X|Z) dZ- \int q(Z) \log \frac{q(Z)}{p(Z)} dZ \\
&= E_{Z \sim q(Z)} [\log p(X|Z)] - KL(q(Z) || p(Z))
\end{aligned}
すなわち、L(q(Z))を対数尤度のqによる期待値とqと事前分布とのKLダイバージェンスに分解することができました。
つまり、L(q(Z))を最大化するということは、対数尤度の期待値を最大化し、q(Z)を事前分布に近づけることです。
また、KL(q(Z)||p(Z|X))の最小化では、事後分布を含んでおり最小化しにくい形になっていますが、L(q(Z))の中には計算できるものしか含まれていないので最大化しやすい形になっています。
変分ベイズ自由エネルギーについて
L(q(Z))を-1倍したものを変分ベイズ自由エネルギーといいます。
すなわち、L(q(Z))の最大化は、変分ベイズ自由エネルギーの最小化になります。
変分ベイズ自由エネルギーは関数を引数とした関数なので、汎関数と呼ばれるものになります。
汎関数の最適化には、通常変分法とよばれるものを使い、フレシェ微分の意味での微分が0の点を探索します。
フレシェ微分をするには定義から、関数空間とそれに付随したノルムが必要になります。
しかし、探す範囲の関数空間を何らかのパラメトリックな確率分布全体と設定してしまえば、探索範囲を有限次元のユークリッド空間に帰着できます。
（分散既知）の1次元ガウス分布の平均の場合の変分ベイズ
具体的なケースとして、分散既知の1次元ガウス分布の平均を変分ベイズによって推定してみましょう。平均の事後分布を以下のように近似できると設定します。
事前分布には、標準コーシー分布を設定します。
ここまでの結論から、KL(q(\mu)||p(\mu|X))の最小化は、L(q(\mu))の最大化に帰着でき、以下のようになります。
L(q) = E_{\mu \sim q(\mu)} [\log p(X|\mu)] - KL(q(\mu) || p(\mu))
q(\mu)の探索範囲を正規分布全体とすると以下のようになります。
q(\mu|a, b) = \frac{1}{(2\pi b)^{1/2}} \exp(-(\mu-a)^2/(2b))
これを使って、L(q)の第1項と2項を計算してみると、以下のようになります。（正規分布の１次２次モーメントに注意する）
\begin{aligned}
&E_{\mu\sim q(\mu)}[\log p(X|\mu)] \\
&= \int_{-\infty}^{\infty} \frac{1}{(2\pi b)^{1/2}} \exp(-(\mu-a)^2/(2b)) \times (-\log (2\pi) -\frac{1}{2}\sum_n (x_n - \mu)^2) d\mu \\
&= - \log (2\pi ) - \frac{1}{2}\int_{-\infty}^{\infty} (N\mu^2 -2 (\sum_n x_n)\mu + \sum_n x_n^2) \frac{1}{(2\pi b)^{1/2}}\exp(-(\mu-a)^2/(2b)) d\mu \\
&= - \log (2 \pi)-\frac{1}{2} \Big( N(a^2 +b) - 2(\sum_n x_n)a + \sum_n x_n^2\Big)
\end{aligned}
\begin{aligned}
&- KL(q(\mu) || p(\mu)) \\
&= - \int q(\mu) \log \frac{q(\mu)}{p(\mu)} d\mu \\
&= - \int \frac{1}{(2\pi b)^{1/2}} \exp(-(\mu-a)^2/(2b)) \log \Big( \pi (1+ \mu^2) \frac{1}{(2\pi b)^{1/2}} \exp(-(\mu-a)^2/(2b))\Big) d \mu \\
&= - \log \pi + \log ( 2 \pi b)^{1/2} \\
&+ \frac{1}{2b} ((a^2 + b) -2a^2+ a^2) - \int_{-\infty}^{\infty} \frac{\log (1 + \mu^2)}{(2\pi b)^{1/2}}\exp(-(\mu-a)^2/(2b)) d \mu \\
&= - \log \pi + \log ( 2 \pi b)^{1/2} + \frac{1}{2}
- \int_{-\infty}^{\infty} \frac{\log (1 + \mu^2)}{(2\pi b)^{1/2}}\exp(-(\mu-a)^2/(2b)) d \mu \\
\end{aligned}
第１項目は解析的に計算できましたが、第２項目には対数の期待値が出てきてしまいました。
このままでは、解析的に計算できないので、苦肉の策としてテイラー展開をします。
log(1 + \mu^2) = \mu^2 -\frac{1}{2} \mu^4 + \frac{1}{3}\mu^6 \cdots
\log(1+\mu^2) = \mu^2とすると、
-KL(q(\mu) || p(\mu)) = - \log \pi + \log ( 2 \pi b)^{1/2} + \frac{1}{2}
-(a^2 + b)
よって、
\begin{aligned}
&L(q;a, b) \\
&= - \log (2 \pi)-\frac{1}{2} \Big( N(a^2 +b) - 2(\sum_n x_n)a + \sum_n x_n^2\Big) \\
&- \log \pi + \log ( 2 \pi b)^{1/2} + \frac{1}{2}-(a^2 + b) \\
&= -\frac{1}{2} \Big( N(a^2 +b) - 2(\sum_n x_n)a \Big) + \log ( 2 \pi b)^{1/2} -(a^2+b) + const.\\
&= -\frac{N+2}{2}(a^2 +b) + (\sum_n x_n)a + \frac{1}{2} \log b + const.
\end{aligned}
となります。a, bについて最大化したいので、以下で勾配が0となる点をもとめてみます。
\begin{aligned}
\frac{\partial L}{\partial a} &= -(N+2)a + \sum_n x_n = 0 \\
\frac{\partial L}{\partial b} &= -\frac{N+2}{2} + \frac{1}{2b} = 0
\end{aligned}
を解くと、
\hat{a} = \frac{\sum_n x_n}{N+2}, \hat{b} = \frac{1}{N+2}
結論として、事後分布は以下のように近似することができました。
p(\mu|X) \sim \mathcal{N}(\mu|\frac{\sum_n x_n}{N+2}, \frac{1}{N+2})
事前分布に、共役事前分布である標準正規分布を与えた場合では、以下の事後分布が得られていましたので、それと比較すると平均と分散をやや小さく推定することがわかりました。
p(\mu|X) = \mathcal{N}(\mu|\frac{\sum_n x_n}{N+1}, \frac{1}{N+1})
最後に
変分ベイズの原理について解説を行いました。
変分ベイズ法をもっとも簡単な例で適用してみましたが、推定するパラメータが複数になる場合は独立性の仮定をうまく利用することで、推定することが可能になるようです。気になる方はご自身でやってみてください。
Discussion