2012年5月19日土曜日

『計算統計学の方法 』読書ノート2;EMアルゴリズムの練習問題(生存時間分析)

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



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


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


「6.4 中途打ち切りデータと単回帰」より.
前回は一応,本を読んでなくても理解できるように文章を書いたけど,なんか帯に短したすきに長し的な感じになってしまったので,今回は計算メモとコードだけを羅列することにした

計算メモ

  • p.98 (6.26) 式 $z^{(k)}_{i} = E_{\boldsymbol{\theta}^{(k)}}[Z_i | \boldsymbol{Y}=\boldsymbol{y}]$ の計算

いま,計算したいのは $i=m+1,\ldots ,n$ の場合なので, \begin{align} z^{(k)}_{i} &= E_{\boldsymbol{\theta}^{(k)}}[Z_i | \boldsymbol{Y}=\boldsymbol{y}] \\ &= E_{\boldsymbol{\theta}^{(k)}}[X_i | X_i > c_i] \end{align} 連続型の分布の時,条件付き期待値$E[X | X > c]$は,
\begin{align} E[X | X > c] &= \frac{\int^{\infty}_{c}x \cdot f(x)dx}{P(X>c)}\\ &= \frac{\int^{\infty}_{c}x \cdot f(x)dx}{1-F(c)} \end{align} なので,p.98 の \begin{align} \int^{\infty}_{c} x \phi \left(\frac{x-\mu}{\sigma} \right) \, dx= \mu\left\{1-\Phi \left(\frac{c-\mu}{\sigma} \right) \right\} +\sigma \phi \left(\frac{c-\mu}{\sigma} \right) \end{align} より \begin{align} z^{(k)}_{i} &= \frac{\int^{\infty}_{c_i} x_i \phi \left(\frac{x_i-\mu^{(k)}_i}{\sigma^{(k)}} \right)\, dx_i}{1-\Phi \left(\frac{c_i-\mu^{(k)}_{i}} {\sigma^{(k)}} \right)} \\ &=\frac{ \mu^{(k)}_{i}\left\{1-\Phi \left(\frac{c_i-\mu^{(k)}_{i}}{\sigma^{(k)}} \right) \right\} +\sigma^{(k)} \phi \left(\frac{c_i-\mu^{(k)}_{i}}{\sigma^{(k)}} \right)}{1-\Phi \left(\frac{c_i-\mu^{(k)}_{i}} {\sigma^{(k)}} \right)} \\ &= \mu^{(k)}_{i}+\sigma^{(k)} \frac{\phi \left(\frac{c_i-\mu^{(k)}_{i}} {\sigma^{(k)}}\right) }{1-\Phi \left(\frac{c_i-\mu^{(k)}_{i}} {\sigma^{(k)}}\right)}. \end{align}

  •  p.99-100 Mステップの計算
 \begin{align} Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(k)})=& - \frac{n}{2} \log2\pi - \frac{n}{2} \log \sigma^2 -\frac{1}{ \sigma^2}\sum^{n}_{i=1}(\mu_i)^2 \\ &- \frac{1}{ \sigma^2} \left\{ \sum^{m}_{i=1} x^2_1 -2 \sum^{m}_{i=1} \mu_i x_i + \sum^{n}_{i=m+1} (z^2_i)^{(k)} -2 \sum^{n}_{i=m+1} \mu_i z^{(k)}_{i}\right\} \end{align} ベクトルはベクトルのまま微分する方法があったような気がするが、よく分からないので一つずつ計算する. \begin{align} \mu_i &= \boldsymbol{v}^{T}\boldsymbol{\beta}\\ &=\beta_0+\beta_1v_{1i}+\beta_2v_{2i} \end{align} なので, \begin{align} 0 &=\frac{\partial}{\partial \beta _0} Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(k)}) \\ &=- \frac{1}{2 \sigma^2} \sum^{n}_{i=1}2\mu _i - \frac{1}{ 2\sigma^2}\left\{ -2 \sum^{m}_{i=1} x_i + -2 \sum^{n}_{i=m+1} z^{(k)}_{i} \right\} \\ &= - \frac{1}{\sigma^2}\left\{ \sum^{n}_{i=1}(\mu_i -x^{\ast}_{i})\right\}, \end{align} \begin{align} 0 &=\frac{\partial}{\partial \beta _1} Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(k)}) \\ &=- \frac{1}{2 \sigma^2} \sum^{n}_{i=1}2\mu _i v_{1i}- \frac{1}{ 2\sigma^2}\left\{ -2 \sum^{m}_{i=1} v_{1i} x_i + -2 \sum^{n}_{i=m+1} v_{1i} z^{(k)}_{i} \right\} \\ &= - \frac{1}{\sigma^2}\left\{ \sum^{n}_{i=1}v_{1i} (\mu_i -x^{\ast}_{i})\right\}, \end{align} 同様, \begin{align} 0 &=\frac{\partial}{\partial \beta _2} Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(k)}) \\ &= - \frac{1}{\sigma^2}\left\{ \sum^{n}_{i=1}v_{2i} (\mu_i -x^{\ast}_{i})\right\}, \end{align} ここで, \begin{align} \boldsymbol{x}^{\ast} =(x^{\ast}_{1},\cdots , x^{\ast}_{n})^T = (c_{1}, \cdots ,c_m, z^{(k)}_{m+1}, \cdots, z^{(k)}_{n})^T . \end{align} まとめると, \begin{align} \boldsymbol{0} &=\frac{\partial}{\partial \boldsymbol{\beta}} Q(\boldsymbol{\theta},\boldsymbol{\theta}^{(k)}) \\ &= - \frac{1}{\sigma^2}\left\{ \sum^{n}_{i=1} \boldsymbol{v}_i (\mu_i -x^{\ast}_{i})\right\}. \end{align} これを解いて, \begin{align} \boldsymbol{0}&= - \frac{1}{\sigma^2}\left\{ \sum^{n}_{i=1} \boldsymbol{v}_i (\mu_i -x^{\ast}_{i})\right\} \\ &= \sum^{n}_{i=1} \boldsymbol{v}_i (\boldsymbol{v}^T_i \boldsymbol{\beta} -x^{\ast}_{i}) \\ &= \sum^{n}_{i=1} \boldsymbol{v}_i \boldsymbol{v}^T_i \boldsymbol{\beta} - \sum^{n}_{i=1} \boldsymbol{v}_i x^{\ast}_{i} \\ \boldsymbol{\beta}^{(k+1)} &= \left( \sum^{n}_{i=1} \boldsymbol{v}_i \boldsymbol{v}^T_i \right)^{-1} \left( \sum^{n}_{i=1} \boldsymbol{v}_i x^{\ast}_{i} \right) . \end{align} その他の計算は本にある通り.

R コード
aml <-data.frame(
  d =c( 1, 1, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0),
  t =c(18, 9,28,31,39,19,45, 6, 8,15,23,28, 7,12, 9, 8, 2,26,10, 4, 3, 4,18, 8, 3,14, 3,13,13,35),
  v1=c( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
  v2=c( 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 1, 0))
aml <-aml[order(aml$d,decreasing =T),]
aml <-transform(aml,x=log(t))
attach(aml)
v <- t(cbind(1, v1, v2))
#初期値
beta <-c(beta_0 =2.43, beta_1=-0.59, beta_2=-0.11)
sigma2 =0.57
theta<-matrix(0, nrow=19, ncol=4)
c <-x[d==1] #打ち切りでない部分
z <-x[d==0] #中途打ち切りデータの部分
n <-length(d)
for(k in 1:19){
  mu <-t(v) %*% beta
  h= dnorm((z-mu[d==0])/sqrt(sigma2))/(1-pnorm((z-mu[d==0])/sqrt(sigma2)))
  zk=mu[d==0]+sqrt(sigma2)*h
  z2k=mu[d==0]^2+sigma2+sqrt(sigma2)*(z+mu[d==0])*h
  xast <-c(c,zk)
  beta2 <- solve(v %*% t(v)) %*% (v %*% xast)
  mu2 =t(v) %*% beta2
  sigma22 <- (sum(c^2-2*mu2[d==1]*c+mu2[d==1]^2)+sum(z2k-2*mu2[d==0]*zk+mu2[d==0]^2))/n
  mu =mu2
  beta =beta2
  sigma2 =sigma22
  theta[k,] =c(t(beta), sigma2)
}
detach(aml)
#library(xtable)
#print(xtable(theta,digits=7),type="html")

実行結果

収束の様子はこんな感じ



$k$ $\beta^{(k)}_{0}$ $\beta^{(k)}_{1}$ $\beta^{(k)}_{2}$ $(\sigma^2)^{(k)}$
1 3.1489857 -0.9346340 -0.2141950 0.6825710
2 3.2484719 -1.0066516 -0.2489479 0.7693717
3 3.2811932 -1.0289521 -0.2613609 0.8030753
4 3.2930759 -1.0369361 -0.2659463 0.8158525
5 3.2974980 -1.0398975 -0.2676593 0.8206739
6 3.2991564 -1.0410071 -0.2683023 0.8224909
7 3.2997799 -1.0414242 -0.2685441 0.8231754
8 3.3000147 -1.0415812 -0.2686352 0.8234332
9 3.3001030 -1.0416404 -0.2686695 0.8235303
10 3.3001363 -1.0416626 -0.2686824 0.8235669
11 3.3001489 -1.0416710 -0.2686872 0.8235806
12 3.3001536 -1.0416742 -0.2686891 0.8235858
13 3.3001554 -1.0416753 -0.2686898 0.8235878
14 3.3001560 -1.0416758 -0.2686900 0.8235885
15 3.3001563 -1.0416760 -0.2686901 0.8235888
16 3.3001564 -1.0416760 -0.2686902 0.8235889
17 3.3001564 -1.0416761 -0.2686902 0.8235889
18 3.3001564 -1.0416761 -0.2686902 0.8235890
19 3.3001564 -1.0416761 -0.2686902 0.8235890


答え合わせ
> sr <- survreg(Surv(x,d)~v1+v2,dist="gaussian", data=aml)
> sr
Call:
survreg(formula = Surv(x, d) ~ v1 + v2, data = aml, dist = "gaussian")

Coefficients:
(Intercept)          v1          v2 
  3.3001564  -1.0416761  -0.2686902 

Scale= 0.907518 

Loglik(model)= -36.8   Loglik(intercept only)= -40.7
 Chisq= 7.71 on 2 degrees of freedom, p= 0.021 
n= 30 
#survreg は分散でなく標準偏差で表示しているので,
> (sr$scale)^2
[1] 0.823589
> sr$it
[1] 4
survreg は 4回のイテレーションで収束している.優秀.

参考文献(孫引き)
$\int^{\infty}_{c} x \phi \left(\frac{x-\mu}{\sigma} \right) \, dx$ と $\int^{\infty}_{c} x^2 \phi \left(\frac{x-\mu}{\sigma} \right) \, dx$ の計算はこの本,解析に使用したAMLデータはこの本に出ているらしい.

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\}\] とは別物である.