作者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