ITPub博客

首页 > 大数据 > 数据挖掘 > R语言实现:混合正态分布EM最大期望估计法

R语言实现:混合正态分布EM最大期望估计法

原创 数据挖掘 作者:tecdat 时间:2020-10-26 15:07:05 0 删除 编辑

原文链接: http://tecdat.cn/?p=4815

 

因为近期在分析数据时用到了EM最大期望估计法这个算法,在参数估计中也用到的比较多。然而,发现国内在R软件上实现高斯混合分布的EM的实例并不多,大多数是关于1到2个高斯混合分布的实现,不易于推广,因此这里分享一下自己编写的k个高斯混合分布的EM算法实现请大神们多多指教。并结合EMCluster包对结果进行验算。

      本文使用的密度函数为下面格式:


   对应的函数原型为 em.norm(x,means,covariances,mix.prop)

x为原数据,means为初始均值,covariances为数据的协方差矩阵,mix.prop为混合参数初始值。

使用的数据为MASS包里面的synth.te数据的前两列

x <- synth.te[,-3]

首先安装需要的包,并读取原数据。

install.packages("MASS")library(MASS)install.packages("EMCluster")library(EMCluster)install.packages("ggplot2")library(ggplot2)Y=synth.te[,c(1:2)]qplot(x=xs, y=ys, data=Y)

然后绘制相应的变量相关图:

R语言实现:混合正态分布EM最大期望估计法

从图上我们可以大概估计出初始的平均点为(-0.7,0.4) (-0.3,0.8)(0.5,0.6)

当然 为了试验的严谨性,我可以从两个初始均值点的情况开始估计

首先输入初始参数:

mustart = rbind(c(-0.5,0.3),c(0.4,0.5))    covstart = list(cov(Y), cov(Y))probs = c(.01, .99)

 

然后编写em.norm函数,注意其中的clusters值需要根据不同的初始参数进行修改,

em.norm = function(X,mustart,covstart,probs){  params = list(mu=mustart, var=covstart, probs=probs)     clusters = 2   tol=.00001  maxits=100  showits=T  require(mvtnorm)  N = nrow(X)  mu = params$mu  var = params$var  probs = params$probs  ri = matrix(0, ncol=clusters, nrow=N)           ll = 0                                          it = 0                                           converged = FALSE                              if (showits)                                     cat(paste("Iterations of EM:", "\n"))  while (!converged & it < maxits) {     probsOld = probs        llOld = ll    riOld = ri           # Compute responsibilities    for (k in 1:clusters){      ri[,k] = probs[k] * dmvnorm(X, mu[k,], sigma = var[[k]], log=F)    }    ri = ri/rowSums(ri)          rk = colSums(ri)                                 probs = rk/N    for (k in 1:clusters){      varmat = matrix(0, ncol=ncol(X), nrow=ncol(X))               for (i in 1:N){        varmat = varmat + ri[i,k] * X[i,]%*%t(X[i,])      }      mu[k,] = (t(X) %*% ri[,k]) / rk[k]      var[[k]] =  varmat/rk[k] - mu[k,]%*%t(mu[k,])      ll[k] = -.5 * sum( ri[,k] * dmvnorm(X, mu[k,], sigma = var[[k]], log=T) )    }    ll = sum(ll)             parmlistold =  c(llOld, probsOld)                parmlistcurrent = c(ll, probs)                 it = it + 1    if (showits & it == 1 | it%%5 == 0)               cat(paste(format(it), "...", "\n", sep = ""))    converged = min(abs(parmlistold - parmlistcurrent)) <= tol  }    clust = which(round(ri)==1, arr.ind=T)         clust = clust[order(clust[,1]), 2]             out = list(probs=probs, mu=mu, var=var, resp=ri, cluster=clust, ll=ll)}

结果,可以用图像化来表示:

qplot(x=xs, y=ys, data=Y) ggplot(aes(x=xs, y=ys), data=Y) +   geom_point(aes(color=factor(test$cluster)))

R语言实现:混合正态分布EM最大期望估计法

R语言实现:混合正态分布EM最大期望估计法

 

 

 类似的其他情况这里不呈现了,另外r语言提供了EMCluster包可以比较方便的实现EM进行参数估计和结果的误差分析。

ret <- init.EM(Y, nclass = 2)em.aic(x=Y,emobj=list(pi = ret$pi, Mu = ret$Mu, LTSigma = ret$LTSigma))#计算结果的AIC

通过比较不同情况的AIC,我们可以筛选出适合的聚类数参数值。(欢迎转载,请注明出处。 )

 

 


来自 “ ITPUB博客 ” ,链接:http://blog.itpub.net/69982319/viewspace-2729821/,如需转载,请注明出处,否则将追究法律责任。

请登录后发表评论 登录
全部评论
拓端tecdat分享最in的大数据资讯,提供“一站式”的数据分析学习和咨询体验,让我们一起做有态度的数据人。关注公众号:拓端数据部落、拓端数据。

注册时间:2020-08-21

  • 博文量
    87
  • 访问量
    35747