论文R语言复现

高斯混合概率在众多领域都有重要应用，依据已知观测数据估计高斯模型中未知参数就显得尤为重要，由于观测值具体来自于高斯分布的哪个分模型是未知的，那么利用传统的极大似然( MLE) 方法进行参数估计就变得十分困难。引入 EM 算法，该方法通过构造分布已知的潜变量对模型进行参数估计，经过多次迭代优化可以使估计值逐渐逼近真实值。本文主要复现该篇文章的实证部分~

一、复现内容

二、复现代码

alpha_1<-0.5
alpha_2<-0.5
mu_1<-0.36
mu_2<-0.25
sigma_1<-1.4
sigma_2<-1.3

fi<-function(x,mu,sigma){
  y<-1/(sqrt(2*pi)*sigma)*exp(-(x-mu)**2/(2*sigma**2))
  return(y)
}

x<-c(0.900,-2.618,1.235,-0.382,0.850,0.443,0.444,0.185,1.502,1.071)


gauss_em<-function(x,alpha,mu,sigma){
  z<-alpha*fi(x,mu,sqrt(sigma))/(alpha_1*fi(x,mu_1,sqrt(sigma_1))+alpha_2*fi(x,mu_2,sqrt(sigma_2)))
  mu<-sum(z*x)/sum(z)
  sigma<-sum(z*((x-mu)**2))/sum(z)
  alpha<-mean(z)  
  parms<-list(alpha,mu,sigma)
  return(parms)
}

gauss_em(x,alpha_1,mu_1,sigma_1)  
gauss_em(x,alpha_2,mu_2,sigma_2)

alpha1_vec<-c()
alpha2_vec<-c()
mu1_vec<-c()
mu2_vec<-c()
sigma1_vec<-c()
sigma2_vec<-c()

for(i in 1:20){
  result_1<-gauss_em(x,alpha_1,mu_1,sigma_1)
  result_2<-gauss_em(x,alpha_2,mu_2,sigma_2)
  alpha_1<-result_1[[1]]
  mu_1<-result_1[[2]]
  sigma_1<-result_1[[3]]
  
  alpha_2<-result_2[[1]]
  mu_2<-result_2[[2]]
  sigma_2<-result_2[[3]]
  
  alpha1_vec<-c(alpha1_vec,alpha_1)
  alpha2_vec<-c(alpha2_vec,alpha_2)
  mu1_vec<-c(mu1_vec,mu_1)
  mu2_vec<-c(mu2_vec,mu_2)
  sigma1_vec<-c(sigma1_vec,sigma_1)
  sigma2_vec<-c(sigma2_vec,sigma_2)
  
}



em_rst<-data.frame(alpha1_vec,alpha2_vec,
                   mu1_vec,mu2_vec,sigma1_vec,sigma2_vec)
em_rst


plot_rst<-function(x,xlab,title){
  plot(1:20,
       x,
       type='o',
       xlab=xlab,
       main=title,
       xaxt='n',
       col='red')
  axis(1,1:20)
}

par(mfrow=c(2,3))
params_names<-c("alpha1","alpha2","mu1","mu2","sigma1","sigma2")
params_list<-list(alpha1_vec,alpha2_vec,mu1_vec,mu2_vec,sigma1_vec,sigma2_vec)

for(i in 1:length(params_list)){
  plot_rst(params_list[[i]],
           params_names[i],
           paste0('EM算法——',params_names[i],'迭代趋势'))
}

注意：论文表格中的 σ \sigma σ 和 μ \mu μ 标记反了！

参考文章：
[1]梁盛楠.基于EM算法的高斯混合模型参数估计[J].黔南民族师范学院学报,2020,40(04):5-8.