2012年5月6日日曜日

『計算統計学の方法 』読書ノート;EMアルゴリズムの練習問題(入門編)

追記
グラフを入れた (6/1)


この本↓ の事例を勉強のために実装してみる. R を使用.


小西 貞則, 越智 義道, 大森 裕浩(2008)『シリーズ〈予測と発見の科学〉 5 計算統計学の方法 ―ブートストラップ,EMアルゴリズム,MCMC―』(朝倉書店)


「6.1 1変量正規分布の場合」より.


これは, 欠測データがある場合に, 正規分布のパラメータ $\mu$ と $\sigma ^2$ を推定する問題である.
以下の様なデータがあるとする.
54.62, 57.98, 49.41, 51.77, 46.72, 43.09, NA, 51.68, 40.33, 50.78, 43.63, NA, 44.89, 52.34, 47.28, 45.73, 49.02, 49.52, 46.41, 45.44, 46.50, 53.57, NA
NA の部分は何らかの理由で値が分からない. この NA の部分を欠測データと呼ぶ.
こういう場合には, EM アルゴリズムが使える.

準備
標記の都合上, 上のデータを並べかえる. そして, 前半の l 個は観測の行われたデータ, 後半の m 個が欠測データであることにする. また, 完全な観測に対応する確率変数ベクトルを $\boldsymbol{X}=(X_1, \ldots,X_{l},X_{l+1}, \ldots, X_n)^T$, 実際に観測の行われたデータに対応する確率変数ベクトルを $\boldsymbol{Y}=(Y_1, \ldots, Y_l)^T$, 欠測データに対応する確率変数ベクトルを $\boldsymbol{Z}=(Z_1,\ldots,Z_m)^T$として,その実現値を小文字で表すことにする(T は転置).
こんな風に.
$x_1$ $x_2$ $x_3$ ...... $x_{18}$ $x_{19}$ $x_{20}$ $x_{21}$ $x_{22}$ $x_{23}$
$y_1$ $y_2$ $y_3$ ...... $y_{18}$ $y_{19}$ $y_{20}$ $z_{1}$ $z_{2}$ $z_{3}$
54.62 57.98 49.41 ...... 45.44 46.50 53.57 ? ? ?

$X_i$ ($i =1, \ldots ,n$) は正規分布に従うと仮定する.
このとき, (仮に欠測部分のデータが観測されたとして)完全データに対する尤度関数 $L^C$ は, \begin{align} L^C (\boldsymbol{\theta},\boldsymbol{x}) &= \prod ^{n}_{i=1} f(x_i) \\ &= \prod ^{n}_{i=1} \frac{1}{\sqrt{2 \pi \sigma ^2}} \exp\left(- \frac{(x_i - \mu)^2}{2 \sigma ^2} \right) \end{align} 対数尤度関数 $l^C$ は \begin{align} l^C (\boldsymbol{\theta},\boldsymbol{x}) &=\log \left( \prod ^{n}_{i=1} \frac{1}{\sqrt{2 \pi \sigma ^2}} \exp\left(- \frac{(x_i - \mu)^2}{2 \sigma ^2} \right)\right) \\ &=\sum^{n}_{i=1}\left( -\frac{1}{2} \log (2 \pi \sigma^2) - \frac{(x_i - \mu)^2}{2 \sigma ^2} \right) \\ &= -\frac{n}{2} \log (2 \pi \sigma^2) - \frac{1}{2 \sigma ^2}\sum^{n}_{i=1} (x_i - \mu)^2 \\ &= -\frac{n}{2} \log (2 \pi \sigma^2) - \frac{1}{2 \sigma ^2}\left\{ \sum^{l}_{i=1} (x_i - \mu)^2 + \sum^{n}_{i=l+1} (x_i - \mu)^2 \right\} \\ &= -\frac{n}{2} \log (2 \pi \sigma^2) - \frac{1}{2 \sigma ^2}\left\{ \sum^{l}_{i=1} (y_i - \mu)^2 + \sum^{m}_{j=1} (z_j - \mu)^2 \right\} \end{align} である.ただし,ここで $\boldsymbol{\theta}= (\mu, \sigma^2)$ である.

E ステップ
E ステップでは,$Q(\boldsymbol{\theta}, \boldsymbol{\theta}^{(k)}) = E_{\boldsymbol{\theta}^{(k)}}[l^C (\boldsymbol{\theta},\boldsymbol{X})|\boldsymbol{Y}=\boldsymbol{y} ]$ を計算する.
\begin{align} Q(\boldsymbol{\theta}, \boldsymbol{\theta}^{(k)})&= E_{\boldsymbol{\theta}^{(k)}}[l^c (\boldsymbol{\theta},\boldsymbol{X})|\boldsymbol{Y}=\boldsymbol{y} ] \\ &= E_{\boldsymbol{\theta}^{(k)}}[ -\frac{n}{2} \log (2 \pi \sigma^2) - \frac{1}{2 \sigma ^2}\left\{ \sum^{l}_{i=1} (Y_i - \mu)^2 + \sum^{m}_{j=1} (Z_j - \mu)^2 |\boldsymbol{Y}=\boldsymbol{y} ]\right\}\\ &= -\frac{n}{2} \log (2 \pi \sigma^2) - \frac{1}{2 \sigma ^2}\left\{ \sum^{l}_{i=1} E_{\boldsymbol{\theta}^{(k)}}[ (Y_i - \mu)^2|\boldsymbol{Y}=\boldsymbol{y}] + \sum^{m}_{j=1} E_{\boldsymbol{\theta}^{(k)}}[ (Z_j - \mu)^2 |\boldsymbol{Y}=\boldsymbol{y} ]\right\} \\ &= -\frac{n}{2} \log (2 \pi \sigma^2) - \frac{1}{2 \sigma ^2}\left\{ \sum^{l}_{i=1} (y_i - \mu)^2 + \sum^{m}_{j=1} E_{\boldsymbol{\theta}^{(k)}}[ (Z_j - \mu)^2]\right\} \end{align} ところで, $E_{\boldsymbol{\theta}^{(k)}}[Z]=\mu^{(k)}$, $E_{\boldsymbol{\theta}^{(k)}}[(Z - \mu^{k})^2]=(\sigma ^2)^{(k)}$ であるから, \begin{align} E_{\boldsymbol{\theta}^{(k)}}[(Z_j -\mu)^2]&= E_{\boldsymbol{\theta}^{(k)}}[(Z_j -\mu^{(k)}+\mu^{(k)}-\mu)^2] \\ &= E_{\boldsymbol{\theta}^{(k)}}[(Z_j -\mu^{(k)})^2 +2(Z_j -\mu^{(k)})(\mu^{(k)} -\mu)+ (\mu^{(k)} -\mu)^2] \\ &= (\sigma ^2)^{(k)} +2E_{\boldsymbol{\theta}^{(k)}}[(Z_j -\mu^{(k)})(\mu^{(k)} -\mu)]+ (\mu^{(k)} -\mu)^2 \\ &= (\sigma ^2)^{(k)} +2E_{\boldsymbol{\theta}^{(k)}}[(Z_j \mu^{(k)} -Z_j \mu -(\mu^{(k)})^2 -\mu\mu^{(k)})]+ (\mu^{(k)} -\mu)^2 \\ &= (\sigma ^2)^{(k)} +2( (\mu^{(k)})^2 - \mu^{(k)} \mu -(\mu^{(k)})^2 -\mu\mu^{(k)})+ (\mu^{(k)} -\mu)^2 \\ &= (\sigma ^2)^{(k)} + (\mu^{(k)} -\mu)^2 \end{align} よって, \[ Q(\boldsymbol{\theta}, \boldsymbol{\theta}^{(k)})= -\frac{n}{2} \log (2 \pi \sigma^2) - \frac{1}{2 \sigma ^2}\left\{ \sum^{l}_{i=1} (y_i - \mu)^2 + m (\mu^{(k)} - \mu)^2 + m(\sigma^2)^{(k)} \right\} \]
M ステップ
M ステップでは,$Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(k)})$ を最大化する $\boldsymbol{\theta}$ を求める.したがって, 以下の方程式をとけば良い. \begin{align} \boldsymbol{0}&= \frac{\partial}{\partial \boldsymbol{\theta}} Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(k)})\\ &= \left( \begin{array}{c} \frac{\partial}{\partial \mu} Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(k)}) \\ \frac{\partial}{\partial \sigma ^2} Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(k)}) \end{array} \right) \\ &= \left( \begin{array}{c} \frac{1}{\sigma^2}\left\{\sum^{m}_{i=1}(y_i - \mu) + m(\mu^{(k)} -\mu) \right\} \\ -\frac{n}{2\sigma^2}+\frac{1}{2(\sigma^2)^2} \left\{ \sum^m_{i=1}(y_i -\mu)^2+m(\mu^{(k)} -\mu)^2 + m(\sigma)^{(k)} \right\} \end{array} \right) \end{align} よって, \begin{align} \boldsymbol{\theta}^{(k+1)} &= \left( \begin{array}{c} \mu ^{(k+1)} \\ (\sigma ^2) ^{(k+1)} \end{array} \right) \\ &= \left( \begin{array}{c} \frac{1}{n}\left\{ \sum^{l}_{i=1}y_i + m \mu ^{(k)}\right\} \\ \frac{1}{n}\left\{ \sum^{l}_{i=1}(y_i - \mu ^{(k+1)})^2+m(\mu ^{(k)}- \mu ^{(k+1)} )^2 +m(\sigma ^2)^{(k)}\right\} \end{array} \right) \end{align} これが EM アルゴリズムの更新式となる.

実装
#データ
x <-c(54.62, 57.98, 49.41, 51.77, 46.72, 43.09, NA, 51.68, 40.33, 50.78, 43.63, NA, 44.89, 52.34, 47.28, 45.73, 49.02, 49.52, 46.41, 45.44, 46.50, 53.57, NA)

y <- x[!is.na(x)]
m <- sum(is.na(x))
n=length(x)

mu1=0; sigma1=1 #初期値
for(k in 1:100){
  mu2= (sum(y)+m*mu1)/n ;sigma2= (sum((y-mu2)^2)+m*(mu1-mu2)^2+m*(sigma1))/n #更新式
  if(abs(mu1-mu2) + abs(sigma1 -sigma2) < 1e-05) break #収束判定
  mu1=mu2
  sigma1=sigma2
}
結果,
> k
[1] 11
> mu2
[1] 48.5355
> sigma2
[1] 17.93601
$\mu=48.5355$, $\sigma^2 =17.9360$ となり, (k=11 で break しているということは,10回目と11回目がほぼ同じ値だったということなので)10回の反復で収束したことが分かる.実はこれらは欠測データを除いたデータの標本平均,標本分散と変わらない.
> mean(y)
[1] 48.5355
#R の関数 var はデフォルトだと不偏分散を求めるので…
> var(y)*(length(y)-1)/length(y)
[1] 17.936
ついでに,欠測観測数 m が 0,1,2,3 と変化するとき,収束状況がどのようになるか見てみよう.
theta =matrix(0, nrow=10, ncol=8)
for(m in 0:3){
  mu1=0; sigma1=1  #初期値
  for(k in 1:10){
    mu2= (sum(y)+m*mu1)/n
    sigma2= (sum((y-mu2)^2)+m*(mu1-mu2)^2+m*(sigma1))/n
    mu1=mu2
    sigma1=sigma2
    theta[k,(2*m+1):(2*m+2)]=c(mu2,sigma2)
  }
}
#print(xtable(theta),type="html")
結果,



m=0 m=1 m=2 m=3
k $\mu^{(k)}$ $(\sigma^2)^{(k)}$ $\mu^{(k)}$ $(\sigma^2)^{(k)}$ $\mu^{(k)}$ $(\sigma^2)^{(k)}$ $\mu^{(k)}$ $(\sigma^2)^{(k)}$
1 48.54 17.94 46.22 123.96 44.12 211.08 42.20 282.91
2 48.54 17.94 48.43 23.23 48.13 37.10 47.71 57.04
3 48.54 17.94 48.53 18.19 48.50 19.69 48.43 23.11
4 48.54 17.94 48.54 17.95 48.53 18.10 48.52 18.61
5 48.54 17.94 48.54 17.94 48.54 17.95 48.53 18.02
6 48.54 17.94 48.54 17.94 48.54 17.94 48.54 17.95
7 48.54 17.94 48.54 17.94 48.54 17.94 48.54 17.94
8 48.54 17.94 48.54 17.94 48.54 17.94 48.54 17.94
9 48.54 17.94 48.54 17.94 48.54 17.94 48.54 17.94
10 48.54 17.94 48.54 17.94 48.54 17.94 48.54 17.94

となり, m が大きくなるにつれて,収束が遅くなる傾向が見て取れる.
また,m = 0 のときは一回で収束している.これは m=0 のときの更新式が, \begin{align} \boldsymbol{\theta}^{(k+1)} &= \left( \begin{array}{c} \mu ^{(k+1)} \\ (\sigma ^2) ^{(k+1)} \end{array} \right) \\ &= \left( \begin{array}{c} \frac{1}{n} \sum^{l}_{i=1}y_i\\ \frac{1}{n}\left\{ \sum^{l}_{i=1}(y_i - \mu ^{(k+1)})^2 \right\} \end{array} \right) \\ &= \left( \begin{array}{c} \bar{y}\\ s^2_y \end{array} \right) \end{align} となるのだから,考えてみれば当たり前のことであった.

その他
更新式をみると,EM アルゴリズムとは,欠測 $\boldsymbol{z}$ の値を, $\boldsymbol{Z}$ の期待値 $E_{\boldsymbol{\theta}^{(k)}}[\boldsymbol{Z}|\boldsymbol{Y}=\boldsymbol{y}]$ で置き換えて計算する方法ではないことがわかる.
$\mu$ についての更新式 \[\mu^{(k+1)} = \frac{1}{n}\left\{ \sum^{l}_{i=1}y_i + m \mu ^{(k)}\right\}\] は, そのように考えたときと一致している. つまり $\mu^{(k+1)}$ と \[\frac{1}{n}\left\{ \sum^{l}_{i=1}y_i + \sum^{m}_{j=1} E_{\boldsymbol{\theta}^{(k)}}[Z_j|\boldsymbol{Y}=\boldsymbol{y}] \right\}\] は等しい. しかし, $\sigma^2$ についての更新式 \[(\sigma^2)^{(k+1)} = \frac{1}{n}\left\{ \sum^{l}_{i=1}(y_i - \mu ^{(k+1)})^2+m(\mu ^{(k)}- \mu ^{(k+1)} )^2 +m(\sigma ^2)^{(k)}\right\} \] は, \[\frac{1}{n}\left\{ \sum^{l}_{i=1}(y_i - \mu ^{(k+1)})^2 + \sum^{m}_{j=1}(E_{\boldsymbol{\theta}^{(k)}}[Z_j|\boldsymbol{Y}=\boldsymbol{y}] - \mu ^{(k+1)})^2\right\}\] とは別物である.

0 件のコメント:

コメントを投稿