PRML復活の呪文 part6 (2.3.4 - 2.3.6)

TL;DR

• 最尤推定による多変量ガウス分布の平均/分散パラメータの推定は1章でやった
  単変量ガウス分布と同様にできる（ただし式変形大変。特に分散の推定）
• ガウス推定による平均/分散パラメータの推定は事前分布に尤度と同じ関数の形である
  共役事前分布を使うと比較的簡単に導出できる

2.3.4 ガウス分布の最尤推定

1.2.4節「ガウス分布」で単変量ガウス分布の平均、分散の最尤推定を行った。本節では多変量ガウス分布の平均、分散の最尤推定を行う。
...といいたいところだが、分散の計算は複雑で演習問題に回っているので割愛。平均の最尤推定の導出のみ行う。

平均の最尤推定に必要な線形代数のテクニック

• 2次形式の微分¹
  $$
  \frac{ \partial }{ \partial x } (x^T A x) = ( A + A^T) x　\tag{A}
  $$
  （こういうものだと思おう。気になった方はググれば賢い人が証明を書いてくれています）
• 共分散行列$\Sigma$は対称行列
• （テキストp.83の下から9行目）対称行列の逆行列も対称行列。なので、$\Sigma^{-1} = ( \Sigma^{-1} )^T$。

導出

多変量ガウス分布は

$$
N( \mathbf{x} | \boldsymbol{\mu}, \Sigma ) = \frac{1}{ (2 \pi)^{D/2} } \frac{1}{| \Sigma |^{1/2}}
\exp \left\{ - \frac{1}{2} (\mathbf{x} - \boldsymbol{\mu} )^T \Sigma^{-1} (\mathbf{x} - \boldsymbol{\mu} ) \right\} \tag{2.43}
$$

なので、対数尤度をとろう。ガウス分布の平均、分散である$\mu, \Sigma$を固定したときにデータ$X$が得られる確率が尤度なので、尤度は$p(X | \mu, \Sigma)$と表記して

$$
\ln p(X | \mu, \Sigma) = - \frac{ND}{2} \ln (2 \pi) - \frac{N}{2} \ln | \Sigma |

• \frac{1}{2} \sum_{n=1}^N ( \mathbf{x} _n - \mu)^T \Sigma^{-1} ( \mathbf{x} _n - \mu) \tag{2.118}
  $$

これを$\mu$について偏微分しよう。$\mu$によらない式(2.118)右辺の第1項、2項は無視してよい。
上述テクニックの式(A)におけるxが$( \mathbf{x} _n - \mu) $であると考えると

\begin{align}

\frac{\partial}{\partial \mu} \ln p(X | \mu, \Sigma) &= - \frac{1}{2} \sum_{n=1}^N
\left\{ \Sigma^{-1} + ( \Sigma^{-1} )^T \right\} ( \mathbf{x} _n - \mu) \frac{ \partial ( \mathbf{x} _n - \mu) } { \partial \mu } \\

&= -\frac{1}{2} \sum_{n=1}^N ( 2 \Sigma^{-1} ) ( \mathbf{x} _n - \mu) (-1) \\
&= \sum_{n=1}^N \Sigma^{-1} ( \mathbf{x} _n - \mu) \tag{2.120}

\end{align}

対数尤度が最大になる$\mu$を求めるために 式(2.120)=0 の方程式を解く。

\begin{align}

\frac{\partial}{\partial \mu} \ln p(X | \mu, \Sigma) 
&= \sum_{n=1}^N \Sigma^{-1} ( \mathbf{x} _n - \mu) =  0 \\

\text{nによらない}\Sigma^{-1}を総和のシグマの外に出して、さらに両辺を\Sigma^{-1}で割って \\
\sum_{n=1}^N  ( \mathbf{x} _n - \mu) =  0 \\
\sum_{n=1}^N  \mathbf{x} _ n = \sum_{n=1}^N   \mu = N \mu \\
\therefore \mu_{ML} &= \frac{1}{N} \sum_{n=1}^N  \mathbf{x} _ n \tag{2.121}

\end{align}

結局、単変量だろうが多変量だろうが、観測値xの総和をデータ数Nで割れば平均の推定値になる、という直感どおりの結果が導出できた。

分散は計算が複雑で、テキストでも本文中では結果のみの掲載なので、ここでも結果だけ示しておく。
（観測値 -　平均）の2乗をNで割る、というこちらも直観通りの結果が導出できている。

$$
\Sigma_{ML} = \frac{1}{N} \sum_{n=1}^N ( \mathbf{x} _n - \mu _ {ML} ) ( \mathbf{x} _n - \mu _ {ML} )^T \tag{2.122}
$$

2.3.5 逐次推定

ガウス分布の平均の最尤推定値は式(2.121)で導出したように、N個のデータの総和（と割り算）で求められる。
では、Nがウン万といった膨大な数になりそうなときも、平均の推定のために過去すべてのデータを記憶する必要があるか？
答えはNoで結論を先に言うと、N個めのデータを観測したときの平均の推定値は、

• (N-1)個めのデータでの平均$\mu_{ML}^{(N-1)}$と
• N個めのデータ$\mathbf{x} _n$

で求められる。

\begin{align}

\mu_{ML}^{(N)} &= \frac{1}{N} \sum_{n=1}^N  \mathbf{x} _n \\
& \text{1からNまでの総和を1からN-1の総和とNの足し算に分割して} \\
&= \frac{1}{N} \mathbf{x} _N + \frac{1}{N} \sum_{n=1}^{N-1}  \mathbf{x} _n \\

& \text{ここで}　\mu_{ML}^{(N-1)} = \frac{1}{N-1} \sum_{n=1}^{N-1}  \mathbf{x} _n より
\sum_{n=1}^{N-1}  \mathbf{x} _n  = (N-1) \mu_{ML}^{(N-1)}であるから \\
&= \frac{1}{N} \mathbf{x} _N + \frac{N-1}{N} \mu_{ML}^{(N-1)} \\
&= \mu_{ML}^{(N-1)} + \frac{1}{N} ( \mathbf{x} _N - \mu_{ML}^{(N-1)} ) \tag{2.126}


\end{align}

N-1個のデータが観測された時点での平均$\mu_{ML}^{(N-1)}$から、N個めのデータが観測されたときはN個めと観測値とN-1での平均の誤差～～1/Nに比例する小さな「誤差信号」といえる量～～$( \mathbf{x} _ N - \mu_{ML}^{(N-1)} )$だけ$\mu_{ML}^{(N-1)}$を修正してN個めのデータの平均$\mu_{ML}^{(N)}$を補正していると解釈できる。

幾何的解釈は以下の図のとおり。

このような逐次推定の汎用的な定式化はRobbins-Monroアルゴリズムというらしいが、テキストの説明ではよく分からないので「そんな手法があるんだな～」ぐらいにしておく。

2.3.6 ガウス分布に対するベイズ推論

2.3.4節ではガウス分布の平均、分散というパラメータを最尤推定で点推定（この値だと1つの値で推定）した。
ここでは、ガウス分布の平均、分散というパラメータを事前分布を導入してベイズ推定にてパラメータの分布（事後分布）を行う。

• 平均が未知で分散が既知の場合
• 分散が未知で平均が既知の場合
• 平均も分散も未知の場合

の３パターンでベイズ推定を行うが、適当な事前分布を選ばないと数学的に解析が厳しくなる。
解析の簡単化のために、以前学習した共役性をもつ共役事前分布 (conjugate prior) を事前分布として選ぶ。
先に上記３パターンの共役事前分布を示しておく。

パターン	事前分布（単変量）	事前分布（多変量）
平均が未知で分散が既知の場合	ガウス分布	ガウス分布
分散が未知で平均が既知の場合	ガンマ分布	ウィシャート分布
平均も分散も未知の場合	ガウス-ガンマ分布	ガウス-ウィシャート分布

単変量ガウス分布で分散が既知、平均が未知の場合

ここではデータが生成されているガウス分布の分散$ \sigma^2 $は既知で、平均$ \mu $が未知であるとして、平均・分散の事後分布$\mu_N, \sigma^2_N$を既知の$\sigma^2, \sigma_0^2, \mu_0$で表すことを考える。

単変量のガウス分布のexpの中身の二次形式の係数から平均・分散を求める。
単変量ガウス分布は

$$
N(x | \mu ) = \frac{1}{ (2 \pi \sigma^2 )^{1/2} } \exp
\left\{ - \frac{1}{2 \sigma^2 } (x - \mu )^2 \right\} \tag{2.42}
$$

なのでexpの中身を式変形すると

$$

• \frac{1}{2} \cdot \frac{x^2}{\sigma^2} + \frac{\mu}{\sigma^2}x + (xによらない項)
  $$

となるので

• $ - \frac{1}{2} x^2 $の係数は分散の逆数になっている
• $ x $の係数は$ \frac{ \mu } { \sigma^2 } $になっている

ことを頭にいれておこう。

事後分布$ p(\mu | \sf{x} ) $は尤度関数$ p( \sf{x} | \mu )$と$\mu$に関する事前分布$ p(\mu) $の積で表せて

$$
p(\mu | \sf{x} ) \propto p( \sf{x} | \mu ) p(\mu) \tag{2.139}
$$

ここで尤度はガウス分布から生成させるデータが観測される確率の積なので

$$
p( \sf{x} | \mu ) = \prod_{n=1}^N p(x_n | \mu ) =

• \frac{1}{ (2 \pi \sigma^2)^{(N/2)} } \exp \left\{ - \frac{1}{2 \sigma^2} \sum_{n=1}^N (x_n - \mu )^2 \right\} \tag{2.137}
  $$

となる。これを（これから求めたい）$\mu$の関数としてみると、expの中身が$\mu$の2次形式になっている。
そこで、事前分布には尤度と同じ形の関数、すなわちexpの中身$\mu$の2次形式となっている関数、つまりガウス分布を仮定する（共役性。これにより式変形が簡単になる）

$$
p(\mu) = N( \mu | \mu_0, \sigma_0^2 ) \tag{2.138}
$$

（ベイズ推定には事前分布を考えないといけないが、事前分布であるガウス分布の平均$\mu_0$の設定に自信がないときは分散$\sigma_0^2$を大きくしておくことで、$\mu_0$の値は大きくばらつきうる、つまり$\mu_0$の設定にあまり自信がないことを表現できる。こうすることで事後分布の分散が大きくなる）

前述のとおり、事後分布のexpの中身を$\mu$についての2次形式に書き直し、係数を比較することで平均と分散を求める。

\begin{align}

p(\mu | \sf{x} ) &\propto　
\exp \left\{ - \frac{1}{2 \sigma^2} \sum_{n=1}^N (x_n - \mu )^2 \right\}
\times
\exp \left\{ - \frac{1}{2 \sigma_0^2 } ( \mu - \mu_0 )^2  \right\} \\

\text{exp}の中身 &= - \frac{1}{2 \sigma^2} \sum_n ( x_n - \mu )^2 - \frac{1}{2 \sigma_0^2} ( \mu - \mu_0)^2 

\end{align}

$ - \frac{\mu^2}{2} $の係数が分散の逆数になっているので

$$
\frac{1}{\sigma_N^2} = \frac{1}{\sigma_0^2} + \frac{N}{\sigma^2} \tag{2.142}
$$

$\mu$の係数が$ \frac{ \mu } { \sigma_N^2 } $になっているので

\begin{align}

\left( \frac{\mu_0}{\sigma_0^2} + \frac{1}{\sigma^2} \sum_n x_n \right) = 
\left( \frac{1}{\sigma_0^2} + \frac{N}{\sigma^2} \right) \mu_N \\

ここで \mu_{ML} = \frac{1}{N} \sum_n x_n \tag{2.143} とおくと

\end{align}

\begin{align}

\left( \frac{\mu_0}{\sigma_0^2} + \frac{N \mu_{ML} }{\sigma^2} \right) &= 
\left( \frac{1}{\sigma_0^2} + \frac{N}{\sigma^2} \right) \mu_N \\

\mu_N &= \left( \frac{\mu_0}{\sigma_0^2} + \frac{N \mu_{ML} }{\sigma^2} \right)
\left( \frac{1}{\sigma_0^2} + \frac{N}{\sigma^2} \right)^{-1} \\

&= \left( \frac{ \mu_0 \sigma^2 + N \mu_{ML} \sigma^2_0}{ \sigma^2 \sigma_0^2} \cdot 
\frac{ \sigma^2 \sigma_0^2 }{ \sigma^2 + N \sigma_0^2 }
\right) \\

&= \frac{ \sigma^2}{ N \sigma_0^2 + \sigma^2 } \mu_0 + 
\frac{ N \sigma_0^2 }{ N \sigma_0^2 + \sigma^2} \mu_{ML} \tag{2.141}

\end{align}

式(2.141), (2.142)の解釈は

• 観測データ数$N=0 $ならば、事後分布の平均は事前分布の平均に一致
• $ N \to \infty $では、事後分布の平均は最尤推定解に一致
• $ N \to \infty $で事後分布の分散は0に近づき、最尤推定解の近くで無限に尖った形となる。つまり、最尤推定により点推定（この値だと1つの値で推定）と等価になる。

多変量で分散が既知、平均が未知の場合

単変量と同じやり方で導出できるらしい。事前分布にガウス分布を想定するといいが、具体的な計算手順は割愛。

逐次更新式とベイズ推論の考え方

2.3.5節でN個のデータを観測したときの平均をN個めのデータ$ \mathbf{x}_n $とN-1個のデータの平均とで表現できることを見た。
この表現をベイズ的な考えで見てみよう。最終データ点$ \mathbf{x} _n $を分けて事後分布の式を書くと
以下になる。

カギかっこ（[ ]）の中身はN-1個のデータを観測したときの平均の事後分布となっているが、N個のデータを観測したときの平均の事後分布を求める際はこの部分を事前分布とみなせば、ベイズの定理によって平均の事後分布が計算できることが分かる。

この考え方は観測データが独立同分布（つまり、N-1個のデータの尤度が総乗でかける）に従うと仮定した際に、どんな問題に対しても使える汎用的なものである。

単変量ガウス分布で平均が既知、分散が未知の場合

ここではデータが生成されているガウス分布の平均$ \mu $は既知で、分散$ \sigma^2 $が未知であるとして、平均・分散の事後分布$\mu_N, \sigma^2_N$を既知の$\sigma_0^2, \mu_0, \mu$で表すことを考える。

ここからの式変形では分散の代わりに精度$ \lambda = 1 / \sigma^2$で考えたほうが簡単になるので、精度を使う。

精度で尤度を書くと

$$
p( \text{X} | \lambda ) = \prod_{n=1}^N N (x_n | \mu, \lambda^{-1} ) \propto
\lambda^{N/2} \exp \left\{ - \frac{\lambda}{2} \sum_{n=1}^N (x_n - \mu)^2 \right\}
\tag{2.145}
$$

これを（これから求めたい）$\lambda$の関数としてみると、($\lambda$のべき乗)×(expの中身が$\lambda$の一次式)となっている。
そこで、事前分布には尤度と同じ形の関数であるガンマ分布を共役事前分布として仮定する（これにより式変形が簡単になる）。

$$
\text{Gam} (\lambda | a, b) = \frac{1}{\Gamma(a)} b^a \lambda^{a-1} \exp ( -b \lambda ) \tag{2.146}
$$

式(2.146)のa, bはガンマ分布の形を決めるハイパーパラメータ。a, bの値を変えると下図のようにいろんな形の分布を表現でき、表現力が高い分布だと理解しておこう。

精度に関する事後分布$p( \lambda | \sf{x}) $は事前分布である式(2.146)と尤度である式(2.145)の積であるので

\begin{align}

p( \lambda | x ) &\propto \lambda^{a_0 - 1} \lambda^{N/2} \exp 
\left\{ - b_0 \lambda - \frac{\lambda}{2} \sum_{n=1}^N (x_n - \mu)^2 \right\}
\tag{2.149} \\

&= \lambda^{ ( a_0 + \frac{N}{2} ) - 1}  \exp 
\left\{ - ( b_0 + \frac{1}{2} \sum_{n=1}^N (x_n - \mu)^2 )\lambda \right\}

\end{align}

この事後分布は共役性から（事前分布に用いた）ガンマ分布$\text{Gam} (\lambda | a_N, b_N)$になっている。

式(2.146)と比べると事後分布のガンマ分布のパラメータ$a_N$は

\begin{align}

a_N &= a_0 + \frac{N}{2} \tag{2.150} \\
&= \frac{1}{2} ( N + 2 a_0 )

\end{align}

となっており、データを観測する前、すなわち$N=0$の時点であたかも$2 a_0$個のデータを観測しているような「事前知識」を与えられていると解釈できる。

式(2.146)と比べると事後分布のガンマ分布のパラメータ$b_N$は

\begin{align}

b_N &= b_0 + \frac{1}{2} \sum_{n=1}^N (x_n - \mu)^2 \\

& ここで \sigma_{ML}^2 = \frac{1}{N} \sum_n ( x_n - \mu)^2、
すなわち\sum_n ( x_n - \mu)^2 = N \sigma_{ML}^2 であることを使って \\

&= b_0 + \frac{N}{2} \sigma_{ML}^2 \tag{2.151} \\

&= \frac{ \sigma_{ML}^2 }{2} \left( N + \frac{2}{ \sigma_{ML}^2 } b_0 \right)
\end{align}

となっており、データを観測する前、すなわち$N=0$の時点であたかも$\frac{2}{ \sigma_{ML}^2 } b_0$個のデータを観測しているようにみえる。$a_N$の解釈で$2 a_0$個のデータが事前に観測されていると解釈したので

\begin{align}

\frac{2}{ \sigma_{ML}^2 } b_0 = 2 a_0 \\
\therefore \sigma_{ML}^2 = \frac{b_0}{a_0}

\end{align}

となる。
まとめると、分散が$\frac{b_0}{a_0}$なるデータが事前に$2 a_0$個観測されていると解釈できる。

###分散も平均も未知の場合

上表に記載したように単変量の場合は事前分布にガウス-ガンマ分布、多変量の場合はガウス-ウィシャート分布を使うと計算が楽になるらしい（共役性）。
ただし、テキストに細かい導出は書いておらず、おそらく計算は大変なので平均もしくは分散が既知の場合と同じような感じで計算できそうだな、ぐらいに思っておこう。

1. Yuuki Saitohさんのスライド(p.31)を参考にさせていただきました。https://www.slideshare.net/YuukiSaitoh/prml2-3 ↩