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

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

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

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

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)


mustart = rbind(c(-0.5,0.3),c(0.4,0.5))

covstart = list(cov(Y), cov(Y))

probs = c(.01, .99)


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)))


ret <- init.EM(Y, nclass = 2)

em.aic(x=Y,emobj=list(pi = ret$pi, Mu = ret$Mu, LTSigma = ret\$LTSigma))#计算结果的AIC

|
14天前
|

R语言参数自抽样法Bootstrap：估计MSE、经验功效、杰克刀Jackknife、非参数自抽样法可视化自测题
R语言参数自抽样法Bootstrap：估计MSE、经验功效、杰克刀Jackknife、非参数自抽样法可视化自测题
18 0
|
14天前
|

R语言GARCH模型对股市sp500收益率bootstrap、滚动估计预测VaR、拟合诊断和蒙特卡罗模拟可视化
R语言GARCH模型对股市sp500收益率bootstrap、滚动估计预测VaR、拟合诊断和蒙特卡罗模拟可视化
36 2
|
14天前
|
vr&ar Python
R语言风险价值：ARIMA，GARCH，Delta-normal法滚动估计VaR（Value at Risk）和回测分析股票数据
R语言风险价值：ARIMA，GARCH，Delta-normal法滚动估计VaR（Value at Risk）和回测分析股票数据
17 0
|
14天前
|
vr&ar
R语言风险价值：ARIMA，GARCH，Delta-normal法滚动估计VaR（Value at Risk）和回测分析股票数据-1
R语言风险价值：ARIMA，GARCH，Delta-normal法滚动估计VaR（Value at Risk）和回测分析股票数据
26 1
|
14天前
|

R语言随机波动模型SV：马尔可夫蒙特卡罗法MCMC、正则化广义矩估计和准最大似然估计上证指数收益时间序列
R语言随机波动模型SV：马尔可夫蒙特卡罗法MCMC、正则化广义矩估计和准最大似然估计上证指数收益时间序列
27 7
|
14天前
|

【视频】什么是梯度下降？用线性回归解释和R语言估计GARCH实例
【视频】什么是梯度下降？用线性回归解释和R语言估计GARCH实例
20 0
|
14天前
|

R语言非参数PDF和CDF估计、非参数分位数回归分析间歇泉、GDP增长数据
R语言非参数PDF和CDF估计、非参数分位数回归分析间歇泉、GDP增长数据
23 0
|
14天前
|

R语言基于Bootstrap的线性回归预测置信区间估计方法分析汽车制动距离|数据分享
R语言基于Bootstrap的线性回归预测置信区间估计方法分析汽车制动距离|数据分享
26 0
|
14天前
|

R语言广义矩量法GMM和广义经验似然GEL估计ARMA、CAPM模型分析股票收益时间序列
R语言广义矩量法GMM和广义经验似然GEL估计ARMA、CAPM模型分析股票收益时间序列
15 0
|
14天前
|

R语言：EM算法和高斯混合模型聚类的实现
R语言：EM算法和高斯混合模型聚类的实现
7 1