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

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

因为近期在分析数据时用到了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)

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

从图上我们可以大概估计出初始的平均点为(-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语言提供了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,我们可以筛选出适合的聚类数参数值。(欢迎转载,请注明出处。)

相关文章
|
14天前
|
资源调度 前端开发 数据可视化
R语言参数自抽样法Bootstrap:估计MSE、经验功效、杰克刀Jackknife、非参数自抽样法可视化自测题
R语言参数自抽样法Bootstrap:估计MSE、经验功效、杰克刀Jackknife、非参数自抽样法可视化自测题
|
14天前
|
资源调度 前端开发 数据可视化
R语言GARCH模型对股市sp500收益率bootstrap、滚动估计预测VaR、拟合诊断和蒙特卡罗模拟可视化
R语言GARCH模型对股市sp500收益率bootstrap、滚动估计预测VaR、拟合诊断和蒙特卡罗模拟可视化
|
14天前
|
vr&ar Python
R语言风险价值:ARIMA,GARCH,Delta-normal法滚动估计VaR(Value at Risk)和回测分析股票数据
R语言风险价值:ARIMA,GARCH,Delta-normal法滚动估计VaR(Value at Risk)和回测分析股票数据
|
14天前
|
vr&ar
R语言风险价值:ARIMA,GARCH,Delta-normal法滚动估计VaR(Value at Risk)和回测分析股票数据-1
R语言风险价值:ARIMA,GARCH,Delta-normal法滚动估计VaR(Value at Risk)和回测分析股票数据
R语言风险价值:ARIMA,GARCH,Delta-normal法滚动估计VaR(Value at Risk)和回测分析股票数据-1
|
14天前
|
算法 Python
R语言随机波动模型SV:马尔可夫蒙特卡罗法MCMC、正则化广义矩估计和准最大似然估计上证指数收益时间序列
R语言随机波动模型SV:马尔可夫蒙特卡罗法MCMC、正则化广义矩估计和准最大似然估计上证指数收益时间序列
|
14天前
|
算法 vr&ar Python
【视频】什么是梯度下降?用线性回归解释和R语言估计GARCH实例
【视频】什么是梯度下降?用线性回归解释和R语言估计GARCH实例
|
14天前
|
算法
R语言非参数PDF和CDF估计、非参数分位数回归分析间歇泉、GDP增长数据
R语言非参数PDF和CDF估计、非参数分位数回归分析间歇泉、GDP增长数据
|
14天前
|
前端开发 数据库
R语言基于Bootstrap的线性回归预测置信区间估计方法分析汽车制动距离|数据分享
R语言基于Bootstrap的线性回归预测置信区间估计方法分析汽车制动距离|数据分享
|
14天前
|
算法 vr&ar Windows
R语言广义矩量法GMM和广义经验似然GEL估计ARMA、CAPM模型分析股票收益时间序列
R语言广义矩量法GMM和广义经验似然GEL估计ARMA、CAPM模型分析股票收益时间序列
|
14天前
|
机器学习/深度学习 算法 数据可视化
R语言:EM算法和高斯混合模型聚类的实现
R语言:EM算法和高斯混合模型聚类的实现

热门文章

最新文章