作者aflilfesy (...)
站內NCTU-STAT95G
標題提供 EM 演算法 (R code)
時間Thu Aug 17 22:08:58 2006
這是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
1F：推 josephw:謝囉!!!真厲害呀!!哈 08/17 22:17
2F：推 uranuss:謝啦 真厲害 @@ 08/17 22:24
3F：→ ym7226:太強 08/17 23:03
4F：推 Y0SHIKI:高手@@a 08/17 23:12
5F：→ aflilfesy:ㄜ~~~我只是用別人的code 高手是別人@@ 08/18 00:33
6F：推 peggyant:我用這個跑出來都不準ㄟ~是error生成的問題嗎? 08/18 00:53
7F：推 aflilfesy:恩~~我們還頗準的.. 有時不准 可能跟你輸入的參數值有關 08/18 00:58
8F：推 Oyin:謝謝^0^ 08/18 01:32
9F：推 peggyant:謝謝 ^+++++++^ 08/18 09:06
10F：推 mrliang:感謝 08/18 11:34
※ 編輯: aflilfesy 來自: 140.113.90.16 (08/18 11:44)
11F：推 mangogogo:哈哈 謝啦^^ 08/18 22:08