R语言实现:混合正态分布EM最大期望估计法
时间:2022-07-22
本文章向大家介绍R语言实现:混合正态分布EM最大期望估计法,主要内容包括其使用实例、应用技巧、基本知识点总结和需要注意事项,具有一定的参考价值,需要的朋友可以参考一下。
原文链接: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)
然后绘制相应的变量相关图:
从图上我们可以大概估计出初始的平均点为(-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,我们可以筛选出适合的聚类数参数值。(欢迎转载,请注明出处。)
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法
- Flutter延时任务、Flutter通过Future与Timer实现延时任务
- DDIA 笔记
- 工作流和状态机
- CentOS 6.x 搭建:Headless Chrome + ChromeDriver + Selenium的爬虫环境系统
- 聊聊dubbo-go的registryAwareCluster
- 同样是空值,null和undefined有什么异同?
- 强大到没朋友的mysql-shell及插件
- android JavaPoet记录
- JavaScript里的分号,你加还是不加?
- 技术干货 | Docker 容器逃逸案例汇集
- 一张千万级别数据的表想做分页,如何优化?
- 一文学会爬虫技巧
- 为什么机器学习应用交易那么难(中)
- 消息队列的消费幂等性如何保证
- js中数组Array.reduce方法介绍及使用场景