EMアルゴリズム

はじめに

千葉大学・株式会社Nospareの川久保です．今回は，EMアルゴリズムについて解説します．EMアルゴリズムは，欠測が含まれたデータに対する観測データの尤度を最大化するアルゴリズムです．

手法の説明

問題設定

$y$を観測データのベクトル，$z$を欠測している確率変数とし，パラメータベクトル$\theta$を用いて$y$の周辺密度関数が$f(y\mid \theta)$とモデル化されているとします．パラメータ$\theta$の推定法として，対数尤度関数$\log f(y\mid \theta)$の最大化，すなわち最尤法を考えます．ここで$f(y \mid \theta)$が直接モデリングされておらず，$y$と$z$の同時密度$f(y,z \mid \theta) = f(y \mid z,\theta)f(z \mid \theta)$がモデリングされていて，$y$の周辺密度を求める際の積分評価
$$
f(y \mid \theta) = \int f(y,z \mid \theta) dz
$$
が困難である場合があります．そのような場合に利用すると便利なのが，今から解説するEMアルゴリズムです．

EMアルゴリズムの手順

EMアルゴリズムはexpectation-maximizationアルゴリズムの略で，期待値計算を行うE-stepと，最大化を行うM-stepを交互に，収束するまで繰り返し計算を行うというアルゴリズムです．$(k+1)$回目のE-stepとM-stepはそれぞれ以下の計算を行います．

• E-step: $k$回目の繰り返しで得られたパラメータ推定値を$\theta^{(k)}$とし，以下の期待値計算を行います．

\begin{align}
Q(\theta, \theta^{(k)}) &= E_{\theta^{(k)}} \left[ \log f(y,z \mid \theta) \mid y \right] \\
&= \int \log \{ f(y,z \mid \theta) \} \ f(z \mid y, \theta^{(k)}) dz.
\end{align}

ただし$E_{\theta^{(k)}}[\cdot \mid y]$は，パラメータ値が$\theta^{(k)}$の$y$を所与とした$z$の条件付分布についての期待値演算を表します．

• M-step: E-stepで評価した$Q(\theta,\theta^{(k)})$を$\theta$について最大化し，最大となる$\theta$の値を$\theta^{(k+1)}$とします．

E-stepとM-stepの反復を収束条件を満たすまで（$\theta^{(k)}$の値がほとんど動かなくなるまで）繰り返し，最終的に得られた$\theta^{(k)}$を$\theta$の推定値とするのがEMアルゴリズムです．ここで$f(y \mid \theta)$が単峰な尤度関数である場合，EMアルゴリズムで得られる$\theta^{(k)}$は，$f(y \mid \theta)$を最大化する最尤推定値に収束することが示されます．

EMアルゴリズムの各繰り返しでは，最尤法の本来の目的関数である周辺対数尤度$\log f(y \mid \theta)$を最大化するかわりに，$\log f(y,z \mid \theta)$（完全データの対数尤度, complete data likelihood）を最大化します．ただし$\log f(y,z \mid \theta)$には欠測している確率変数$z$が含まれるので，$z$の値を条件付期待値で補って代入しているイメージです．厳密には$E[z \mid y]$を完全データの対数尤度に代入しているわけではなく，完全データの対数尤度を$z \mid y$の分布で期待値をとっています．

EMアルゴリズムの適用

EMアルゴリズムは欠測データ解析だけでなく，潜在変数を導入することで完全データの尤度を書くことが簡単になるようなモデルにおける最尤法にも用いられます．代表的な例は混合分布の推定で，それぞれの観測が混合分布のどのcomponentからの標本なのかを示す潜在変数ベクトル$z$を導入し，この$z$を欠測データとみなしてEMアルゴリズムを適用します．混合分布でのEMアルゴリズムは多くの解説書で紹介されているので，この記事では線形混合モデルの未知パラメータ推定への適用例を紹介したいと思います．

線形混合モデルへの適用

以前の記事で紹介した以下のような線形混合モデルの一般形を考えます．

$$
y = X\beta + Zb + \varepsilon, \quad b \sim \mathrm{N}_m(0,G), \quad \varepsilon \sim \mathrm{N}_n(0,R), \tag{1}
$$
このモデルおける未知パラメータを，変量効果$b$を欠測データとみなしたEMアルゴリズムで推定します．ここで簡単化のために，$G$や$R$に含まれる分散成分を全て既知とし，未知パラメータは$\beta$だけだとします．もちろん分散成分が未知の場合も，回帰係数と分散成分とをあわせてEMアルゴリズムで推定できます．最大化の目的関数の基礎となる完全データの対数尤度は，

\begin{align}
\log f(y,b \mid \beta) =& \ \log f(y \mid b,\beta) + \log f(b) \\
=& -{n \over 2}\log(2\pi) -{1 \over 2}\log |R| - {1 \over 2} (y - X\beta - Zb)^\top R^{-1}(y - X\beta - Zb) \\
& \ -{m \over 2}\log(2\pi) - {1 \over 2}\log|G| - {1 \over 2}b^\top G^{-1}b
\end{align}

と書けます．$(k+1)$回目のE-stepでは$Q(\beta, \beta^{(k)}) = E_{\beta^{(k)}}[\log f(y,b \mid \beta) \mid y]$を求め，M-stepは，$Q(\beta, \beta^{(k)})$を$\beta$について最大化します．$Q(\beta, \beta^{(k)})$の1階導関数は，
$$
{d \over d\beta}Q(\beta, \beta^{(k)}) = -(X^\top X)\beta + X^\top (y - ZE_{\beta^{(k)}}[b \mid y])
$$
なので，実際にはE-stepでは$b \mid y$の分布の一次モーメントだけがわかれば十分です．正規分布の性質から，
$$
E_{\beta^{(k)}}[b \mid y] = GZ^\top \Sigma^{-1}(y - X\beta^{(k)}) =: \tilde{b}^{(k)}
$$
とわかります．ただし，$\Sigma$は$y$の周辺分布の分散で，$\Sigma = ZGZ^\top + R$です．よって$(k+1)$回目のM-stepで求まる最適解は，
$$
\beta^{(k+1)} = (X^\top X)^{-1} X^\top (y - Z\tilde{b}^{(k)})
$$
です．

数値実験

(1)式のモデルの特殊形として，以下のモデルを考えます．
$$
y_{ij} = \beta_0 + x_{ij} \beta_1 + b_i + \varepsilon_{ij}, \quad b_i \sim \mathrm{N}(0,\tau^2), \quad \varepsilon_{ij} \sim \mathrm{N}(0,\sigma^2), \quad (i=1,\dots,m, \ j=1,\dots,n_i)
$$
このモデルにおいて$\tau^2 = 0.5^2, \sigma^2 = 1^2$を既知として，$(\beta_0, \beta_1) = (1,2)$をEMアルゴリズムで推定してみます．以下がRのコード例です．

library(Matrix)

RR <- 100  # 実験の繰り返し数
m  <- 30
ni <- 5
n  <- m * ni

# 真のパラメータの値
tau2 <- 0.5^2
sigma2 <- 1^2
beta <- c(1,2)

# 行列XとZの作成
X <- cbind(1, rbinom(n, 1, 0.3))
i <- 1:n
j <- rep(1:m, each = ni)
Z <- sparseMatrix(i, j, x = 1, dims = c(n,m))

# Sigmaとその逆行列の計算
Sigma <- tau2 * Z %*% t(Z) + sigma2 * diag(n)
Sigma_inv <- solve(Sigma)

res <- data.frame(beta0 = rep(0,RR), beta1 = rep(0,RR), iter = rep(0,RR))
# 実験の繰り返し
for(rr in 1:RR){
  # 乱数の発生
  b <- rnorm(m, sd = sqrt(tau2))
  y <- X %*% beta + Z %*% b + rnorm(n, sd = sqrt(sigma2))
  
  d <- 1
  iter <- 0
  beta_new <- c(0,0)  # 初期値
  # EMの繰り返し
  while(d > 0.001){
    iter <- iter + 1
    beta_old <- beta_new  # 1期前のM-stepで求めたbeta^kの値
    
    # E-step
    b_tilde <- tau2 * t(Z) %*% Sigma_inv %*% (y - X %*% beta_old)
    
    # M-step
    beta_new <- solve(t(X) %*% X) %*% t(X) %*% (y - Z %*% b_tilde)
    
    # 収束判定
    d <- sqrt( sum( ( beta_new - beta_old )^2 ) )
  }
  # 結果の格納
  res[rr,1:2] <- beta_new
  res[rr,3] <- iter
}

実験を100回繰り返し，その結果をresというデータフレームに格納しています．resの1列目と2列目は$(\beta_0,\beta_1)$の推定値，3列目はEMアルゴリズムが何回で収束したかを保存しています．収束判定は，

d = \| \beta^{(k+1)} - \beta^{(k)} \|

が0.001より小さくなれば収束したと判定しています．結果は以下です．

> mean(res$beta0)
[1] 1.00784
> mean(res$beta1)
[1] 2.025926
> var(res$beta0)
[1] 0.0196979
> var(res$beta1)
[1] 0.03556165
>
> table(res$iter)

12 13 
26 74 

真のパラメータ値が$(\beta_0, \beta_1) = (1,2)$であることから，妥当な結果が得られていると言えそうです．EMアルゴリズムは，12回か13回の繰り返しで収束しています．

おわりに

株式会社Nospareには，統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospare までお問い合わせください．