真の\(\mu\)、\(\sigma\)があるものとしてデータ\(\boldsymbol{x}=(x_1,x_2,\cdots,x_n)\)が最も起こりやすくなるように
\(\mu\)、\(\sigma\)を計算するのが最尤推定だった。
最尤推定においては、\(\mu_{ML}\)、\(\sigma_{ML}\)が先にあって、その結果として\(\boldsymbol{x}\)があると考えている。
逆に、データ\(\boldsymbol{x}=(x_1,x_2,\cdots,x_n)\)があって、
その結果として\(\mu\)、\(\sigma\)が確率的な分布をもって存在すると考える。
\(\mu\)、\(\sigma\)はデータ\(\boldsymbol{x}=(x_1,x_2,\cdots,x_n)\)次第で確率的に決まる。
この分布を事後確率分布と呼んでいて、
事後確率分布を最大にする\(\mu_{MAP}\)、\(\sigma_{MAP}\)を求めるのが最大事後確率推定(MAP推定)。
まずは、事後確率分布と事前確率分布の関係をまとめてみて、
最大化すべき事後確率分布とは何なのか、詳細に追ってみる。
最尤推定
\(N\)個のデータ\(\boldsymbol{x}=(x_1,x_2,\cdots,x_n)\)が観測されたとして、
観測値にはガウスノイズが含まれる(=観測値は正規分布に従う)。
\(\boldsymbol{x}\)が互いに独立であるならば、同時確率は周辺確率、つまり各々の確率の積で表せた。
\(\boldsymbol{x}\)は既知のデータ、統計量である\(\mu\)、\(\sigma\)は変数であるから、
同時確率は\(\mu\)、\(\sigma\)を変数とする関数として表せた。
\begin{eqnarray}
p(x | \mu, \sigma^2) = \prod_{n=1}^{N}N(x_n|\mu,\sigma^2)
\end{eqnarray}
既にある観測値\(\boldsymbol{x}\)を最も高い確率で実現するための統計量\(\mu\)、\(\sigma\)が
一つ存在するという仮定のもと\(\boldsymbol{x}\)を推定した。
\(p(x|\mu,\sigma^2)\)は、\(\mu\)、\(\sigma\)がとある値であるときに、
\(\boldsymbol{x}\)が発生する確率として捉えることができる(条件付き確率)。
事後確率分布
逆に、\(\boldsymbol{x}\)が得られたときに\(\mu,\sigma^2\)が得られる確率\(p(\mu,\sigma^2|x)\)を考える。
\(\boldsymbol{x}\)自体が確率的に生起するデータであり、\(\mu,\sigma^2\)が\(\boldsymbol{x}\)次第でブレる。
そんな状況。
ベイズの定理を使って因果を入れ替えると以下の通り。
\begin{eqnarray}
p(\mu, \sigma^2 | x) = \frac{p(\boldsymbol{x|}\mu,\sigma^2)p(\mu,\sigma^2)}{p(\boldsymbol{x})}
\end{eqnarray}
分子の\(p(\boldsymbol{x|}\mu,\sigma^2)\)は尤度。\(p(\mu,\sigma^2)\)は事前分布。
データ\(\boldsymbol{x}\)が事前分布\(p(\mu,\sigma^2)\)に従うと仮定したもので、
最初からそれがわかっていれば苦労しない訳だから、\(p(\mu,\sigma^2)\)はわからない。
\(p(\boldsymbol{x})\)はデータ\(\boldsymbol{x}\)が発生する確率。
観測することから得られる確率。\(\mu,\sigma^2\)に対して定数。
この\(p(\mu, \sigma^2 | x)\)を最大化する\(\mu, \sigma^2\)を求めようというのが、
最大事後確率推定。
\(p(\mu, \sigma^2 | x)\)の式において、分母の\(p(\boldsymbol{x})\)は\(\mu, \sigma^2\)に対して定数だから、
\(p(\mu, \sigma^2 | x)\)を最大化する\(\mu, \sigma^2\)は、分母を払った\(p(\boldsymbol{x|}\mu,\sigma^2)p(\mu,\sigma^2)\)も最大化する。
だから、代わりに\(p(\boldsymbol{x|}\mu,\sigma^2)p(\mu,\sigma^2)\)の最大化を考える。
\begin{eqnarray}
p(\mu, \sigma^2) p(x | \mu, \sigma^2) &=& N(x_n|\mu,\sigma^2) \prod_{n=1}^{N}N(x_n|\mu,\sigma^2) \\
&=& N(x_n|\mu,\sigma^2) p(x | \mu_{ML}, \sigma^2_{ML})
\end{eqnarray}
尤度関数\(p(\boldsymbol{x|}\mu,\sigma^2)\)の最尤解\(\mu_{ML}\)と\(\sigma^2_{ML}\)の求め方は前回記事。
データさえあれば標本平均、標本分散を計算して終了。
\begin{eqnarray}
\mu_{ML} &=& \frac{1}{N}\sum_{i=1}^{N}x_i \\
\sigma^2_{ML} &=& \frac{1}{N} \sum_{i=1}^N (x_i -\mu_{ML})^2
\end{eqnarray}
尤度\(\prod_{n=1}^{N}N(x_n|\mu,\sigma^2)\)は計算により値が決まる。
\begin{eqnarray}
\prod_{n=1}^{N}N(x_n|\mu,\sigma^2) &=& \prod_{j=1}^{n} \frac{1}{\sqrt{2\pi} \sigma_{ML}} exp \left( – \frac{1}{2} \left( \frac{x_j -\mu_{ML}}{\sigma_{ML}} \right)^2 \right)
\end{eqnarray}
共役事前確率分布
さて…
事後確率と事前確率の関係をわかりやすく書くと以下の通り。
だが、事前確率\(p(\mu,\sigma^2|x)\)は分からない。
\begin{eqnarray}
p(\mu,\sigma^2|x) = \frac{p(x|\mu,\sigma^2)}{p(x)} p(\mu,\sigma^2)
\end{eqnarray}
事後確率分布\(p(\mu,\sigma^2|x)\)が事前確率分布\(p(\mu,\sigma^2)\)と同じ確率分布であれば、
事前確率分布と尤度の計算から事後確率分布を求められる。
そういう操作をするために、事後確率分布と同じ形の事前確率分布を仮に作っておいて、
まず事後確率分布と仮の事前確率分布を結ぶ式を立てる。
仮に立てた事前確率分布を、共役事前確率分布と言うらしい。
いきなり事前確率分布を\(\mu\)、\(\sigma\)だけをもつ別の確率分布に置き換えることはできないから、
確率分布を制御する超パラメータ\(\alpha,\beta,\gamma,\cdots\)を使って表現する。
正規分布の共役事前分布は逆ガンマ分布を使う。
逆ガンマ分布は4個の超パラメータ\(\alpha,\beta,\gamma,\delta\)を持つ。
\begin{eqnarray}
p(\mu,\sigma^2) = \frac{\sqrt{\gamma}}{\sigma \sqrt{2\pi}} \frac{\beta^\alpha}{\Gamma[\alpha]} \left( \frac{1}{\sigma^2} \right)^{\alpha+1} exp \left( – \frac{2\beta + \gamma (\delta – \mu)^2 }{2\sigma^2} \right)
\end{eqnarray}
なので、ベイズの式は超パラメータ\(\alpha,\beta,\gamma,\delta\)を使って以下のように書ける。
左辺の事後確率分布と右辺の共役事前確率分布は同じ分布であって、
右辺には定数\(1/p(x)\)と尤度\( \prod_{n=1}^{N}N(x_n|\mu,\sigma^2)\)をかけた形になっている。
\begin{eqnarray}
p(\mu,\sigma^2|x) &=& \frac{p(x|\mu,\sigma^2)}{p(x)} p(\mu,\sigma^2) \\
&=& \frac{1}{p(x)} \prod_{n=1}^{N}N(x_n|\mu,\sigma^2) \frac{\sqrt{\gamma}}{\sigma \sqrt{2\pi}} \frac{\beta^\alpha}{\Gamma[\alpha]} \left( \frac{1}{\sigma^2} \right)^{\alpha+1} exp \left( – \frac{2\beta + \gamma (\delta – \mu)^2 }{2\sigma^2} \right)
\end{eqnarray}
左辺の事後確率を最大化する\(\mu,\sigma\)を見つけるわけだけども、
それは、尤度と共役事前確率で表された右辺を最大化する\(\mu,\sigma\)を見つけるのと同じ。
対数をとると積の部分が和になって計算しやすいから対数をとる。
やり方は\(\mu\)、\(\sigma\)で偏微分したものが0とする。
\begin{eqnarray}
\frac{\partial}{\partial \mu} p(\mu,\sigma^2|x) &=& \frac{\partial}{\partial \mu} \frac{p(x|\mu,\sigma^2)}{p(x)} p(\mu,\sigma^2) \\
&=& \frac{\partial}{\partial \mu} \frac{1}{p(x)} \prod_{n=1}^{N}N(x_n|\mu,\sigma^2) \frac{\sqrt{\gamma}}{\sigma \sqrt{2\pi}} \frac{\beta^\alpha}{\Gamma[\alpha]} \left( \frac{1}{\sigma^2} \right)^{\alpha+1} exp \left( – \frac{2\beta + \gamma (\delta – \mu)^2 }{2\sigma^2} \right) \\
&=& 0
\end{eqnarray}
ちょっと計算しようとすると目眩がするけど、超パラメータ\(\alpha,\beta,\gamma,\delta\)はそのまま残して、
\(\mu,\sigma\)を\(\alpha,\beta,\gamma,\delta\)を含む式で表す。
\begin{eqnarray}
\mu_{MAP} &=& \frac{\sum_{i=1}^{N}x_i + \gamma \delta}{N+\gamma} \\
&=& \frac{N\bar{x} + \gamma \delta}{N+\gamma}
\end{eqnarray}
同様に、
\begin{eqnarray}
\frac{\partial}{\partial \sigma^2} p(\mu,\sigma^2|x) &=& \frac{\partial}{\partial \sigma^2} \frac{p(x|\mu,\sigma^2)}{p(x)} p(\mu,\sigma^2) \\
&=& \frac{\partial}{\partial \sigma^2} \frac{1}{p(x)} \prod_{n=1}^{N}N(x_n|\mu,\sigma^2) \frac{\sqrt{\gamma}}{\sigma \sqrt{2\pi}} \frac{\beta^\alpha}{\Gamma[\alpha]} \left( \frac{1}{\sigma^2} \right)^{\alpha+1} exp \left( – \frac{2\beta + \gamma (\delta – \mu)^2 }{2\sigma^2} \right) \\
&=& 0
\end{eqnarray}
計算できないけど
\begin{eqnarray}
\sigma^2_{MAP} &=& \frac{\sum_{i=1}^{N} (x_i-\mu)^2 + 2\beta + \gamma (\delta-\mu)^2}{N+3+2\alpha} \\
&=& \frac{N\sigma^2 + 2\beta + \gamma (\delta-\mu)^2}{N+3+2\alpha}
\end{eqnarray}
最尤推定による解と並べてみる。
平均、分散ともに、最尤解の分母、分子をパラメータ\(\alpha,\beta,\gamma,\delta\)で調整したものになっている。
\begin{eqnarray}
\mu_{ML} &=& \bar{x} \\
\mu_{MAP} &=& \frac{N\bar{x}+\gamma \delta}{N+\delta} \\
\sigma_{ML}^2 &=& \bar{\sigma^2} \\
\sigma_{MAP}^2 &=& \frac{N\sigma^2 + 2\beta + \gamma (\delta-\mu)^2}{N+3+2\alpha}
\end{eqnarray}
