廿TT

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

EMアルゴリズムを用いたときの観測情報行列の求め方(混合正規分布)

Louis の公式

EMアルゴリズムの弱点は, 推定量の分散の評価を直接的には得られないことにあった.

しかし, 完全なデータが得られたという条件の下で対数尤度を用いて情報行列を計算する方法がいくつか考案されている.

Louis (1982) は, 以下の関係式を示した.

 \displaystyle I(\theta; y) = E[B_c(x; \theta) | y; \theta ]_{\theta=\hat \theta} - \\
\displaystyle \quad E[S_c(x | \theta)S_{c}^{T}(x | \theta) | y; \theta ]_{\theta=\hat \theta}

ここで  S_c(x; \theta) は完全データの対数尤度の 1 次導関数ベクトル,  B_c(x; \theta) は 2 次導関数行列の負値である.

例題:混合正規分布

混合正規分布のパラメータをEMアルゴリズムで推定する場合の更新式は 混合正規分布のパラメータ推定(あるいは EM アルゴリズム練習問題) - 廿TT に書いた.

今回は, 正規分布の分散は 1 で既知とする. 混合数は K=2とする.

上のエントリと同様の記法を用いて, 完全観測に基づく対数尤度 l^C (\boldsymbol{\theta ,x})
\displaystyle l^C (\boldsymbol{\mu,\xi}) = \sum^{n}_{i=1}\sum^{K}_{j=1} z_{ij} \log \phi (y_i-\mu_j) +\\
 \displaystyle \quad \sum^{n}_{i=1}\sum^{K}_{j=1} z_{ij} \log \xi _j
となる. ここで\phi(x) は標準正規分布の密度関数.

\displaystyle \frac{\partial}{\partial \xi_j}l^C = \sum_{i=1}^{n} \left(\frac{z_{ij}}{\xi_j} - \frac{z_{iK}}{\xi_K} \right)

これを 0 と置いて解くと,  \hat \pi_j = \sum z_i/n なので,

\displaystyle \frac{\partial^2}{\partial \xi_j}l^C = \sum_{i=1}^{n} \left(-\frac{z_{ij}}{\xi_j^2} - \frac{z_{iK}}{\xi_K^2} \right)\\
\displaystyle = -\frac{n}{\hat \xi_j \hat \xi_K}

また,

\displaystyle \frac{\partial}{\partial \mu_j}l^C = \sum_{i=1}^{n}z_{ij}\{-(y_i-\mu_j)\}

\displaystyle \frac{\partial^2}{\partial \mu_j}l^C = -\sum_{i=1}^{n}z_{ij}=-n\hat \xi_j

二次の導関数は, z_{ij} について線形の式になっているので, E[B_c(x; \theta) | y; \theta ]_{\theta=\hat \theta} を求めるには単に z_{ij} をその条件付き期待値で置き換えればよい.

E[z_{ij}|y]= w_{ij} と書くことにする.

E[S_c(x | \theta)S_{c}^{T}(x | \theta) | y; \theta ]_{\theta=\hat \theta} を順番に計算していく.

E[S_c(x | \theta)S_{c}^{T}(x | \theta) | y; \theta ]_{\theta=\hat \theta}l,m 成分を  S_{lm} と書くことにする.

S_{11} は,

\displaystyle \frac{\partial}{\partial \xi_j}l^C = \sum_{i=1}^{n} \left(\frac{z_{ij}}{\xi_j} - \frac{z_{iK}}{\xi_K}\right) \\
\displaystyle \quad = \sum_{i=1}^{n} \left( \frac{z_{i1}(1-\xi_1) - (1-z_{i1})\xi_1}{\xi_1 \xi_K}\right)\\
\displaystyle \quad = \sum_{i=1}^{n} \left( \frac{z_{i1} - \xi_1}{\xi_1 \xi_K}\right)

より,

 \displaystyle S_{11} = V(\sum_{i=1}^{n} \left( \frac{z_{i1} - \xi_1}{\xi_1 \xi_K}\right))\\
=\displaystyle \sum_{i=1}^{n} \left( \frac{w_{i1}(1-w_{i1})}{\xi_1^2\xi_K^2}\right)

さらに,

 \displaystyle S_{22} =V( \sum_{i=1}^{n}z_{ij}\{-(y-\mu_j)\})\\
=\displaystyle \sum_{i=1}^{n}w_{i1}(1-w_{i1})(x-\mu_1)^2

 \displaystyle S_{33} = \sum_{i=1}^{n}w_{i1}(1-w_{i1})(x-\mu_2)^2

であり, また多項分布の共分散より,

 \displaystyle S_{12} = {\rm Cov}(\sum_{i=1}^{n} \left( \frac{z_{i1} - \xi_1}{\xi_1 \xi_K}\right), \sum_{i=1}^{n}z_{i2}\{-(y-\mu_1)\}) \\
 \displaystyle = \sum_{i=1}^{n}w_{i1}w_{i2}(y-\mu_1)/(\xi_1 \xi_K)

 \displaystyle S_{23} = {\rm Cov}( \sum_{i=1}^{n}z_{i2}\{-(y-\mu_1)\}), \sum_{i=1}^{n}z_{i3}\{-(y-\mu_2)\}) \\
 \displaystyle = -\sum_{i=1}^{n}w_{i1}w_{i2}(y-\mu_1)(y-\mu_2)

と各要素が求まる.

シミュレーション

観測情報行列から推定量の標準誤差が近似的に求まる.

これをもとに95%信頼区間を作り, 真の値が信頼区間に入る割合(coverage probability)を求める.

\xi=0.4, \mu_1=0, \mu_2=2 とし, サンプルサイズを100, 500, 1000, と変化させ, 各2000回のシミュレーションをしたときの結果は以下のようになった.

n \xi \mu_1 \mu_2
100 0.92 0.93 0.94
200 0.94 0.94 0.94
500 0.94 0.94 0.95

サンプルサイズが小さいときはちょっと甘めに出てしまうようだ.

R のコード

doEM <- function(y,xi0,mu0,K,maxit=1000,tol=1e-6){
  xi_cur <- xi0
  mu_cur <- mu0
for(i in 1:maxit){
  tmp <- sapply(1:K,function(j)xi_cur[j]*dnorm(y,mu_cur[j]))
  w <- tmp/apply(tmp,1,sum)
  xi_new <- apply(w,2,mean)
  mu_new <- apply(w*y,2,mean)/xi_new
  if(all(abs(c(xi_new, mu_new) - 
             c(xi_cur, mu_cur)) < tol)){break}
  xi_cur <- xi_new
  mu_cur <- mu_new
}
  list(xi=xi_new, mu=mu_new, w=w, iter=i)
}

Infomat <- function(y,w,xi,mu,K){
  n <- length(y)
  Ic <-diag(c(n/prod(xi),n*xi))
  wprod <- apply(w,1,prod)
  ydifmu <- sapply(mu,function(x)(y-x))
  Im11 <- sum(wprod)/prod(xi^2)
  Im22_33 <- apply(wprod*sapply(mu,function(x)(y-x)^2),2,sum)
  Im23 <- -sum(wprod*apply(ydifmu,1,prod))
  Im12_13 <- apply(wprod*(ydifmu),2,sum)/prod(xi)
  Im <-diag(c(Im11,Im22_33))
  Im[2,3] <- Im23
  Im[3,2] <- Im23
  Im[1,2] <- Im12_13[1]
  Im[2,1] <- Im12_13[1]
  Im[1,3] <- -Im12_13[2]
  Im[3,1] <- -Im12_13[2]
  Ic-Im
}

sim_mix_norm <- function(B=10000,level=0.95,N,true_xi,true_mu,K){
  q <- qnorm((1-level)/2)
  mu_boot <- xi_boot <- matrix(NA,B,K)
    cover_mu1 <- cover_mu2 <- cover_xi <- numeric(B)
#  pb <- txtProgressBar(min = 1, max = B, style = 3)
  for(i in 1:B){
    Z <- sample.int(K,N,replace = TRUE,prob = true_xi)
    y <-rnorm(N,true_mu[Z])
    est <-doEM(y,true_xi,true_mu,K)
    ord <-order(est$mu)
    mu_boot[i,] <- est$mu[ord]
    xi_boot[i,] <- est$xi[ord]
    se <- sqrt(diag(solve(Infomat(y,est$w[,ord],est$xi[ord],est$mu[ord],K))))
    xi_lower <-est$xi[ord][1] + q*se[1]
    xi_upper <-est$xi[ord][1] - q*se[1]
    mu1_lower <-est$mu[ord][1] + q*se[2]
    mu1_upper <-est$mu[ord][1] - q*se[2]
    mu2_lower <-est$mu[ord][2] + q*se[3]
    mu2_upper <-est$mu[ord][2] - q*se[3]
    cover_xi[i] <- xi_lower <= true_xi[1] & true_xi[1] <= xi_upper
    cover_mu1[i] <- mu1_lower <= true_mu[1] & true_mu[1] <= mu1_upper
    cover_mu2[i] <- mu2_lower <= true_mu[2] & true_mu[2] <= mu2_upper
#    setTxtProgressBar(pb, i)
  }
  list(xi_boot=xi_boot,mu_boot=mu_boot,cover_xi=cover_xi,
       cover_mu1=cover_mu1,cover_mu2=cover_mu2)
}
wrap_sim <-function(N){
res = sim_mix_norm(B=2000,N=N,true_xi = c(0.4,0.6),true_mu = c(0,2),K=2)
cbind(
mean(res$cover_xi),
mean(res$cover_mu1),
mean(res$cover_mu2))
}

library(parallel)
result_sim <-mclapply(c(100,200,500),wrap_sim,mc.cores = 3)
names(result_sim)<-c("xi","mu1","mu2")
print(cbind(n=c(100,200,500),
      sapply(result_sim,function(x)round(x,2))))

参考文献・関連エントリ

http://www.markirwin.net/stat221/Refs/louis1982.pdf

The EM Algorithm and Extensions (Wiley Series in Probability and Statistics)

The EM Algorithm and Extensions (Wiley Series in Probability and Statistics)

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

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

http://ebsa.ism.ac.jp/ebooks/sites/default/files/ebook/1881/pdf/vol3_ch9.pdf

EMアルゴリズムを用いたときの観測情報行列の求め方(混合二項分布) - 廿TT