提供 EM 演算法 (R code)

[aflilfesy (...)]

這是EM演算法 適用於
Model: Yi ~ mixed Normal(Σw_i*N(mu_i,sigma^2_i)) ,i=1,2,3,...,g , Σw_i=1
可估計 w_i,mu_i,sigma_i

--

ECM.normix=function(y, g, tol)
{
begin = proc.time()[1]
n=length(y)
class=kmeans(y, g)$cluster
w=mu=sigma=rep(NA, g)
 for(i in 1:g)
 {
 w[i]=length(y[class==i])/n;
 mu[i]=mean(y[class==i])
 sigma[i]=sqrt(var(y[class==i]))
 }
yg=matrix(rep(y, g), ncol=g)
den.M= dnorm(t(yg), mean=mu, sd=sigma)
w.den=t(w*den.M)
log.like.old = sum( log(rowSums(w.den)) )
cat("log.like=", log.like.old, "\n")
iter=0

repeat{
iter=iter+1
# E-step
Z=w.den/rowSums(w.den)
# M-step
n.z=colSums(Z)
w=n.z/n
mu=colSums(Z*y)/n.z
yg.cent = t(t(yg)-mu)
sigma=sqrt(colSums(Z*yg.cent^2)/n.z)
den.M= dnorm(t(yg), mean=mu, sd=sigma)
w.den=t(w*den.M)
log.like.new=sum(log(rowSums(w.den)))
if(abs(log.like.new-log.like.old)/abs(log.like.old)<tol) break
log.like.old=log.like.new
cat("Iter=", iter, "\n")
cat("w = ", round(w, 4), "\n")
cat("mu = ", round(mu, 4), "\n")
cat("sigma = ", round(sigma, 4), "\n")
cat("log.like = ", round(log.like.new, 12), "\n")
cat(paste(rep("-", 30), collapse = ""), "\n")
 }

end = proc.time()[1]
cat("The spending time is ", round( end-begin, 6 ) , "sec", "\n")
cat(paste(rep("-", 50), collapse = ""), "\n")
den.M= dnorm(t(yg), mean=mu, sd=sigma)
w.den=t(w*den.M)
Z = w.den /rowSums(w.den)
list(w=w, mu=mu, sigma=sigma, Z=Z)
}

# Compute the MLEs
mle=ECM.normix(y, g, tol=1e-10)
w=mle$w; mu=mle$mu; sigma=mle$sigma; Z=mle$Z

--
以上就是在R環境下的指令
其中原文的
w就是我們作業的p,
g就是我們作業的q,
y就是我們作業的error,

tol 是指tolerance

希望各位受用了...

這是我們班同學提供的喔,因為這次作業她剛好有這方面的code,
但是她害羞,不好意思公佈姓名..

大家加油吧~^^
--

一個小例子
E = rnorm(1000)
mle=ECM.normix(y=c(E),g=2,tol=1e-10)   #這邊要特別加c 用成向量表示
w=mle$w; mu=mle$mu; sigma=mle$sigma


--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 140.113.90.16
推 josephw:謝囉!!!真厲害呀!!哈
推 uranuss:謝啦  真厲害 @@
→ ym7226:太強
推 Y0SHIKI:高手@@a
→ aflilfesy:ㄜ~~~我只是用別人的code  高手是別人@@
推 peggyant:我用這個跑出來都不準ㄟ~是error生成的問題嗎?
推 aflilfesy:恩~~我們還頗準的.. 有時不准 可能跟你輸入的參數值有關
推 Oyin:謝謝^0^
推 peggyant:謝謝 ^+++++++^
推 mrliang:感謝
※ 編輯: aflilfesy       來自: 140.113.90.16        (08/18 11:44)