読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

多項分布のパラメータ推定:モンテカルロ EM アルゴリズムの例題

R

『R によるモンテカルロ法入門』p.177 の例題をやります.
『R によるモンテカルロ法入門』は書き方に省略が多くて意味がわからないところがけっこうあるので『計算統計学の方法』pp.82-pp.88 もあわせて参照しています.

Rによるモンテカルロ法入門

Rによるモンテカルロ法入門

計算統計学の方法―ブートストラップ・EMアルゴリズム・MCMC (シリーズ予測と発見の科学 5)

計算統計学の方法―ブートストラップ・EMアルゴリズム・MCMC (シリーズ予測と発見の科学 5)

どっちか買うとしたら『計算統計学の方法』がおすすめです。

本文

観測値 \boldsymbol{x} = \left(x_1,x_2,x_3,x_4\right)=(125,18,20,34) が得られたとます.
対応する確率変数 \boldsymbol{X} は, 各カテゴリの発現確率が \displaystyle \left(\frac{1}{2}+\frac{\theta}{4},\frac{1}{4}(1-\theta),\frac{1}{4}(1-\theta),\frac{\theta}{4} \right) の多項分布に従います.

第一カテゴリの観測 x_1 はふたつのカテゴリからの観測 z_1, z_2 の和 x_1=z_1+z_2 だと考えます. 確率変数 \boldsymbol{Z}=(Z_1,Z_2,Z_3,Z_4,Z_5) は発現確率 \displaystyle \left(\frac{1}{2},\frac{\theta}{4},\frac{1}{4}(1-\theta),\frac{1}{4}(1-\theta),\frac{\theta}{4} \right) の多項分布に従うとします.

\boldsymbol{Z} の確率関数は,

 \displaystyle p(\boldsymbol{z}|\theta) = \frac{n!}{z_1!z_2!z_3!z_4!z_5!}\{\frac{1}{2}\}^{z_1}\{\frac{1}{4}\theta\}^{z_2}\\
 \displaystyle ~~~ \times \{\frac{1}{4}(1-\theta)\}^{z_3}\{\frac{1}{4}(1-\theta)\}^{z_4}\{\frac{1}{4}\theta\}^{z_5}

です.

完全な観測 z_1,z_2,z_3,z_4,z_5 が得られたとして, \boldsymbol{Z}対数尤度は

\displaystyle l(\theta;\boldsymbol{z}) = \log \frac{n!}{z_1!z_2!z_3!z_4!z_5!} \\
~~~ - z_1 \log 2 + z_2( \log \theta - \log 4) \\
~~~ +z_3 (-\log 4 +\log(1-\theta))\\
~~~ + z_4 (-\log 4 +\log(1-\theta)) \\
~~~ +z_5  (-\log 4 +\log\theta)\\
\propto (z_2+z_5)\log \theta + (z_3+z_4)\log(1-\theta) .

この対数尤度の条件付き期待値は

Q(\theta;\theta_0)\\
= E_{\theta_0}[(Z_2+Z_5)\log \theta + (Z_3+Z_4)\log(1-\theta)|\boldsymbol{X}=\boldsymbol{x} ] \\
=  (E_{\theta_0}[Z_2|\boldsymbol{X}=\boldsymbol{x} ]+x_4)\log \theta + (x_2+x_3)\log(1-\theta).

ここで, 観測 \boldsymbol{X}=\boldsymbol{x} が得られたという条件のもとでの Z_2 の条件付き分布は, サイズ X_1, 確率

\displaystyle \left(\frac{1}{4}\theta_0 \right) \big/ \left(\frac{1}{2} + \frac{1}{4}\theta_0 \right)=\frac{\theta_0}{2+\theta_0}

の二項分布になります.

よって

\displaystyle E_{\theta_0}[Z_2|\boldsymbol{X}=\boldsymbol{x} ]=x_1\frac{\theta_0}{2+\theta_0}.

Q(\theta;\theta_0) \theta に関して最大化するには, \theta微分して 0 と置いて解けばよいのでパラメータの更新式は,

 \displaystyle \left( \frac{\theta_0}{2+\theta_0}x_1 + x_4 \right) \big/ \left(\frac{\theta_0}{2+\theta_0}x_1+x_2+x_3 + x_4  \right)

となります.

モンテカルロ EM では期待値  E_{\theta_0}[Z_2|\boldsymbol{X}=\boldsymbol{x} ] を以下の経験平均で置き換えます.

 \bar Z_2 = \frac{1}{m} \sum^{m}_{i=1}Z_{2i}.

m \bar Z_2 はサイズ m x_1, 確率 \theta_0/(2+\theta_0) の二項分布からシミュレーションします.

練習問題

上記の設定で, \theta_0 =0.5 を初期値として EM 系列をプロットします.

また MCEM 系列を1000回シミュレーションして, 分散のレンジをプロットします.

y <- c(125,18,20,34)
EM <- function(y,theta0){
  theta_history <- theta0
  for(i in 1:100){
    z <- y[1]*(theta0/(2+theta0))
    den <- z+y[4]
    num <- z+y[2]+y[3]+y[4]
    theta1 <- den/num
    theta_history <- c(theta_history,theta1)
    if(abs(theta1-theta0)<1e-6)break
    theta0 <- theta1
  }
  list(theta=theta1,theta_history=theta_history,iterations=i)
}
res1 <-EM(y,0.5)
it <-res1$iterations
m <- 100
theta0 <- 0.5
theta_history <- matrix(,it+1,1000)
theta_history[1,] <- theta0
set.seed(1234)
for(i in 1:it){
  z <- rbinom(1000,size=m*y[1],prob=theta0/(2+theta0))/m
  den <- z+y[4]
  num <- z+y[2]+y[3]+y[4]
  theta1 <- den/num
  theta_history[i+1,] <- theta1
  theta0 <- theta1
}
upper=apply(theta_history,1,max)
lower=apply(theta_history,1,min)
ggplot() +
  geom_line(aes(x=0:it,y=res1$theta_history))+
  geom_ribbon(aes(x=0:it,ymax=upper,ymin=lower),alpha=0.3)+
  xlab("iterations")+ylab("MCEM sequences")+
  theme_grey(20)

f:id:abrahamcow:20160320034836p:plain