廿TT

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

EMアルゴリズム:ランダムな欠測があるデータから多変量正規分布のパラメータ推定

今日の川柳

http://ebsa.ism.ac.jp/ebooks/sites/default/files/ebook/1881/pdf/vol3_ch9.pdf
アルゴリズムはこれに出ています。

R のコードだけ貼ります。

library(mvtnorm)
EM_mGauss <- function(y,muini,sigmaini) {
for(j in 1:1000){
  z <- y
  m <-is.na(y)
  N<-nrow(y)
  K<-ncol(y)
  Sigma_new <- matrix(0,K,K)
  for(i in 1:N){
    if(any(m[i,])){
      sigma00 <- sigmaini[m[i,],m[i,]]
      sigma01 <- sigmaini[m[i,],!m[i,]]
      sigma10 <- sigmaini[!m[i,],m[i,]]
      sigma11 <- sigmaini[!m[i,],!m[i,]]
      z[i,m[i,]] <- muini[m[i,]] + sigma01%*%solve(sigma11)%*%(y[i,!m[i,]]-muini[!m[i,]]) 
      Sigma_new[m[i,],m[i,]] <- Sigma_new[m[i,],m[i,]] + z[i,m[i,]]%*%t(z[i,m[i,]]) + sigma00 -sigma01%*%solve(sigma11)%*%sigma10
      Sigma_new[m[i,],!m[i,]] <- Sigma_new[m[i,],!m[i,]] + z[i,m[i,]]%*%t(z[i,!m[i,]])
      Sigma_new[!m[i,],m[i,]] <- Sigma_new[!m[i,],m[i,]] + z[i,!m[i,]]%*%t(z[i,m[i,]])
      Sigma_new[!m[i,],!m[i,]] <- Sigma_new[!m[i,],!m[i,]] + z[i,!m[i,]]%*%t(z[i,!m[i,]])
    }else{
      Sigma_new <- Sigma_new + z[i,]%*%t(z[i,])
    }
  }
  mu_new <- colMeans(z)
  Sigma_new <- Sigma_new/N - mu_new%*%t(mu_new)
  
  if(all(abs(muini-mu_new)<1e-8)&all(abs(sigmaini-Sigma_new)<1e-8)){
    break
  }
  muini <- mu_new
  sigmaini <- Sigma_new  
}
return(list(mu=mu_new,sigma=Sigma_new,iter=j))
}
sigma <- matrix(c(1,0.1,0.6,
                  0.1,2,0.5,
                  0.6,0.5,3),byrow=TRUE,ncol=3)
y <- rmvnorm(n=100, mean=c(1,2,3), sigma=sigma)
y[sample.int(100*3,30)] <- NA
y <- y[apply(y, 1, function(x)!all(is.na(x))),]
muini <- c(0,0,0)
sigmaini <- diag(1,3)
EM_mGauss(y,muini,sigmaini)

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

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