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

廿TT

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

EMアルゴリズムの練習:右打ち切りデータからの指数分布のパラメータ推定

アルゴリズム

独立に平均 \mu の指数分布に従う大きさ n の標本 x_i があり, 内 n_1 は完全な観測が得られ, n_2=n-n_1 は右打ち切りされているとする.

完全な観測が得られたとして, このときの対数尤度は,

\log L (\mu) = -n \log \mu - \frac{1}{\mu} \left\{\sum_{i=1}^{n_1}x_i + \sum_{j=1}^{n_2} z_j \right\}

である. ただし, z_i は観測されない打ち切られた線分の長さを拡大して表したものである.

右打ち切りされた観測値 x^c_j が得られた下での, z_j の期待値は  \mu + x^c_j である( 左切断指数分布の期待値 - 廿TT ).

z_j をこの期待値で置き換え(E ステップ), 完全な対数尤度を最大化する  \mu^{(k+1)} は,

 \mu^{(k+1)}=(\sum x_i +n_2\mu^{(k)}+\sum x^c_j)/n

である(M ステップ).

適当な初期値  \mu^{(0)} からはじめて, 繰り返し計算を回すことで, 推定値  \hat \mu が得られる.

ところで, この場合の指数分布のパラメータの最尤推定量は,

 \left( \sum x_i +\sum x^c_j \right)/ n_1

である.

続く R によるシミュレーションで, この最尤推定と EM アルゴリズムによる推定を比較する.

シミュレーション

シミュレーションの詳細は以下の R コードを参照.

type2censor <- function(n,mu,w){
  z <-rexp(n,1/mu)
  d <- z<w
  x <- ifelse(z>w,w,z)
  data.frame(x,d)
}
estExp  <- function(dat,mu0,w){
  x <-dat$x[dat$d]
  y <-dat$x[!dat$d]
  n1 <-sum(dat$d)
  n2 <-sum(!dat$d)
  for(i in 1:100){
    muhat <- (sum(x)+n2*mu0+sum(y))/(n1+n2)
    if(abs(muhat-mu0)<1e-6)break
    mu0 <-muhat
  }
  muhat
}
N <- 10000
est1 <- estEM <- numeric(N)
for(i in 1:N){
  dat <-type2censor(20,3,3)
  est1[i] <- estExp(dat,1,3)
  estEM[i] <- sum(dat$x)/sum(dat$d)
}
plot(est1,estEM)

縦軸に EM アルゴリズムによる推定値, 横軸に最尤推定による推定値を取りプロットすると, 両者が一致していることがわかる.

f:id:abrahamcow:20161101222315p:plain