廿TT

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

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

Louis の公式

Louis (1982) はEMアルゴリズムを用いたときの観測情報行列に関する, 以下の関係式を導いた.

 \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 次導関数行列の負値である.

さらに未観測のデータ x が, 多項分布(カテゴリカル分布)に従うインジケータの場合,

 \displaystyle I(\theta; y) =S_c(\hat x | \theta)S_{c}^{T}(\hat x | \theta)

となると述べている.

ここで  \hat x = E[x | y; \theta ] である.

なぜそうなるのかはわからない. だれか教えてください.

混合二項分布の場合

完全な観測に基づく対数尤度は

\displaystyle l^C (\boldsymbol{\mu,\xi}) = \sum^{n}_{i=1}\sum^{K}_{j=1} z_{ij} \log f (y_i; n_i, p_j) +\\
 \displaystyle \quad \sum^{n}_{i=1}\sum^{K}_{j=1} z_{ij} \log \xi _j

である. ここで f(y) は二項分布の確率関数.

一回微分すると,

\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 \frac{\partial}{\partial p_j}l^C = \sum_{i=1}^{n} z_{ij}\left(\frac{y_i}{p_j} - \frac{n_i-y_i}{1-p_j} \right)

シミュレーション

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

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

 \xi=(0.2,0.3,0.5), p = (0.1,0.5,0.7) とし, サンプルサイズを100, 200, 500, と変化させ, 各2000回のシミュレーションをしたときの結果は以下のようになった.

n \xi_1 \xi_2 p_1 p_2 p_3
100 0.94 0.94 0.94 0.94 0.93
200 0.93 0.94 0.95 0.95 0.96
500 0.95 0.95 0.94 0.95 0.96

おおむね, 名目上の数字と近い値になることが確認できる.

R のコード

doEM_mix_binom <- function(y,n,xi_inits,p_inits,K,maxit=1000,tol=1e-6){
  z_update <- function(y,n,xi,p,K){
    tmp <-sapply(1:K,function(j)xi[j]*dbinom(y,n,p[j]))
    tmp/apply(tmp,1,sum)
  }
  xi_update <- function(z){
    apply(z,2,mean)
  }
  p_update <- function(y,n,z){
    apply(z*y,2,mean)/apply(z*n,2,mean)
  }
  xi_cur <- xi_inits; p_cur <- p_inits
  for(i in 1:maxit){
    z <- z_update(y,n,xi_cur,p_cur,K)
    xi_new <- xi_update(z) 
    p_new <- p_update(y,n,z)
    # Stop iteration if the difference between the current and new estimates is less than a tolerance level
    if( all(abs(xi_cur - xi_new) < tol) & all(abs(p_cur - p_new) < tol) ){
      break
    }
    # Otherwise continue iteration
    xi_cur <- xi_new; p_cur <- p_new
  }
  
  list(z=z, xi=xi_new, p=p_new, iter=i)
}

Sfunc <-function(j,y,n,z,p.est,xi.est,K){
  if(j %in% 1:(K-1)){
    z[,j]/xi.est[j]-z[,K]/xi.est[K]
  }else{
    z[,j-(K-1)]*(y/p.est[j-(K-1)]-(n-y)/(1-p.est[j-(K-1)]))
  }
}
infomat <- function(z,y,n,p.est,xi.est,K){
  mat <-matrix(NA, K*2-1, K*2-1)
  for(i in 1:(K*2-1)){
    for(j in 1:(K*2-1)){
      mat[i,j] <- sum(Sfunc(i,y,n,z,p.est,xi.est,K)*
                        Sfunc(j,y,n,z,p.est,xi.est,K))
    }
  }
  mat
}
# N=100
K=3
true_xi = c(0.2,0.3,0.5)
true_p = c(0.1,0.5,0.7)
q <- qnorm((1-0.95)/2,lower.tail = FALSE)
simfunc <- function(N){
cover_xi <- matrix(NA,2000,K-1)
cover_p <- matrix(NA,2000,K)
for(i in 1:2000){
true_z <- sample.int(K,N,TRUE,true_xi)
n = rpois(N,50)
y <- rbinom(N,n,true_p[true_z])
est <-doEM_mix_binom(y,n,true_xi,true_p,K)
ord <-order(est$p)
z <- est$z[,ord]
p.est <-est$p[ord]
xi.est <- est$xi[ord]
mat <-infomat(z,y,n,p.est,xi.est,K)
se <-sqrt(diag(solve(mat)))
cover_xi[i,] <- xi.est[-K] - q*se[1:2] < true_xi[-K] & true_xi[-K] < xi.est[-K] + q*se[1:2]
cover_p[i,] <- p.est - q*se[K:(K*2-1)] < true_p & true_p < p.est + q*se[K:(K*2-1)]
}
c(apply(cover_xi,2,mean),apply(cover_p,2,mean))
}

result_sim <-sapply(c(100,200,500),simfunc)

result_mat <-cbind(c(100,200,500),t(result_sim))
colnames(result_mat) <- c("n","xi_1","xi_2","p_1","p_2","p_3")
print(round(result_mat,2))

参考文献・関連エントリ

Louis, T. A. (1982). Finding the observed information matrix when using the EM algorithm. Journal of the Royal Statistical Society. Series B (Methodological), 226-233.

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