EM与高斯混合分布

EM与高斯混合分布
03 May 2013

从问题入手，现在有一样本，这样产生：


> set.seed(123)
> data <- NULL
> ind <- NULL
> library(mvtnorm)
> mean1 <- c(-2, -2)
> mean2 <- c(2, 2)
> sigma <- matrix( c(3.5, 0, 0,3),2,2)
> for (i in 1:300) {
+   if (runif(1) < 0.6) {
+     data <- rbind(data, rmvnorm(1, mean1, sigma))
+     ind <- c(ind, 0)
+   } else {
+     data <- rbind(data, rmvnorm(1, mean2, sigma))
+     ind <- c(ind, 1)
+   }
+ }
> plot(data, xlab = "", ylab = "")

样本是这个样子的如下图，我们可以掌握的信息是，这是两个二元正态分布f1,f2混合后的样本,其中混合参数为p.即总体密度函数为：

\[f(x1,x2)=(1-p)\times f1(x1, x2) + p \times f2(x1, x2)\]

acma

我们不仅需要估计f1,f2中的$\mu, \sigma$参数，还要估计p值。传统的极大似然估计对此束手无策，可以采用EM算法，EM算法有两步，分别是E步和M步，E步是固定$\mu, \sigma$估计p,M步是固定p估计$\mu, \sigma$，E步和M步不断循环，最终会逼近真实结果。


pi <- rep(0.4, 300)
u1 <- rbind(NULL, c(6, 6)) #储存均值
u2 <- rbind(NULL, c(2, 3)) #储存均值
sigma <-  list()
sigma[[1]] <- matrix( c(2, 0, 0,1),2,2) #储存协方差
pi.matrix <- NULL
pi.matrix <- rbind(pi.matrix, pi) #储存标记
for (i in 1:300) {
  pi <- pi.matrix[i, ]
  
  pi <- pi * dmvnorm(data, u2[i,], sigma[[i]]) / ((1-pi) * dmvnorm(data, u1[i,], sigma[[i]]) +
    pi * dmvnorm(data, u2[i,], sigma[[i]]))#E步
  print(pi)
  pi.matrix <- rbind(pi.matrix, pi)

  #M步
  u1.temp <- colSums((1-pi) * data )/(n- sum(pi))
  print(u1.temp)
  u2.temp <- colSums(pi * data) / sum(pi)
  print(u2.temp)
  #sigma.temp <- ((colSums((1-pi)*(data-u1.temp)^2)+colSums(pi*(data-u2.temp)^2))/n)^0.5
  sigma.temp <- matrix(apply(sapply(1:n, function(i) (1-pi[i])*((data[i, ]-u1.temp)%*%t(data[i, ]-u1.temp)) +
                                         pi[i]*((data[i, ]-u2.temp)%*%t(data[i, ]-u2.temp)) ), 1, mean), nrow = 2)
  u1 <- rbind(u1, u1.temp)
  u2 <- rbind(u2, u2.temp)
  print(1)
  sigma[[i+1]] <- sigma.temp
}

均值的收敛轨迹，红色点和绿色点的路径即为均值的收敛，它们逐渐逼近样本的中心： acma

参数估计：


> u1[300,]
[1] 2.417023 2.127362
> u2[300,]
[1] -1.956428 -1.955138
> sigma[[300]]
           [,1]       [,2]
[1,]  3.1124007 -0.1008249
[2,] -0.1008249  3.1173438

结果与生成数据时的参数接近。而且还能准确估计出分布参数p