EM算法
leengsmile
2016年9月24日
EM 算法
本文档介绍如何在R语言中，通过EM算法，估计高斯混合模型的参数。首先通过简单的例子，用简单的程序描述EM算法估计高斯混合模型参数的过程，再介绍如何使用第三方包实现相应的估计。
为保证数据结果的可重复性，设置随机数种子
set.seed(123)
首先需要高斯混合模型的数据
n <- 1000
mean_s <- c(1, 7)
y <- sample(c("head", "tail"), size = n, replace = TRUE, prob = c(0.25, 0.75))
x <- rnorm(n = 1000, mean = mean_s[1])
tails <- y %in% c("tail")
x[tails] <- rnorm(sum(tails), mean = mean_s[2])
上述产生的混合模型是均值分别为1和7，标准差均为1的混合模型，且混合的概率为(0.25, 0.75)。 也就是说，混合模型中的观测值有0.25的概率来自于均值为1的高斯分布，有0.75的概率来自于均值为7的高斯分布。
其概率概率密度函数为
require(lattice)
densityplot(~x, par.settings = list(plot.symbol = list(col = factor(y))))
数据分布明显，呈现很好的可分性。下面需要估计对应的正太分布的均值，以及混合概率。这里假定方差恒定且相等，均为1。
probs <- c(0.5, 0.5)
mu_s <- c(0, 1)
sigma_s <- c(1, 1)
for(i in seq(10))
{
ps <- matrix(0, ncol = 2, nrow = n)
for(j in seq(2))
{
ps[, j] <- probs[j] * dnorm(x, mean = mu_s[j], sd = sqrt(sigma_s[j]))
}
ps <- ps / rowSums(ps)
for(j in seq(2))
{
sigma_s[j] <- sum( ps[, j] * (x - mu_s[j])^2) / sum(ps[, j])
mu_s[j] <- sum(x * ps[, j]) / sum(ps[, j])
probs[j] <- mean(ps[, j])
}
}
cat(
"mean:", mean_s, "\n",
"sigma:", sqrt(sigma_s), "\n",
"prob:", probs, "\n",
sep = " "
)
## mean: 1 7
## sigma: 0.9415133 0.9913065
## prob: 0.2469451 0.7530549
估计的均值分别为0.9748991, 6.9999522，混合概率为0.2469451, 0.7530549。经过10次迭代，估计值已经很接近精确值。
可以将上述求解的过程封装成一个函数
gmm <- function(x, mean, sd = NULL)
{
num <- length(mean)
if(is.null(sd))
{
sd <- rep(1, num)
}
epsilon <- 1e-4
probs <- rep(1/num, num)
mu_s <- mean
sigma_s <- sd ^ 2
n <- length(x)
while(TRUE)
{
ps <- matrix(0, ncol = num, nrow = n)
for(j in seq(num))
{
ps[, j] <- probs[j] * dnorm(x, mean = mu_s[j], sd = sqrt(sigma_s[j]))
}
ps <- ps / rowSums(ps)
sigma_s_p <- sigma_s
for(j in seq(num))
{
sigma_s[j] <- sum( ps[, j] * (x - mu_s[j])^2) / sum(ps[, j])
mu_s[j] <- sum(x * ps[, j]) / sum(ps[, j])
probs[j] <- mean(ps[, j])
}
if (max(abs(sigma_s_p - sigma_s)) < epsilon)
{
break
}
}
return (list(mu = mu_s, sd = sqrt(sigma_s), prob = probs))
}
上述封装的函数gmm用以估计高斯混合模型的参数，包括各个混合成分的均值mu，标准差sd，混合成分的概率prob。
用gmm估计前面提到的数据x
gmm(x, mean = c(0, 1), sd = c(1, 1))
## $mu
## [1] 0.9749062 6.9999548
##
## $sd
## [1] 0.9415230 0.9913026
##
## $prob
## [1] 0.2469457 0.7530543
在R语言中，可以通过mixtools包实现上述的EM算法估计过程。
首先载入mixtools
require(mixtools)
mixtools的normalmixEM可以实现高斯混合模型的参数估计。
em <- normalmixEM(x, mu = c(0, 1), sigma = c(1, 1), sd.constr = c(1, 1))
## number of iterations= 6
估计的结果中，lambda含有混合比例，mu是混合成分的均值。
print(em$lambda)
## [1] 0.2471721 0.7528279
print(em$mu)
## [1] 0.9777237 7.0008419
从上面的结果可知，normalmixEM的估计结果与前面编写的程序估计出的参数一致。
plot(em, whichplots = 2)
同时，返回的em变量，含有许多有用的信息
str(em)
## List of 9
## $ x : num [1:1000] 6.179 0.0063 6.6927 1.7511 -0.5092 ...
## $ lambda : num [1:2] 0.247 0.753
## $ mu : num [1:2] 0.978 7.001
## $ sigma : num [1:2] 1 1
## $ loglik : num -1955
## $ posterior : num [1:1000, 1:2] 6.14e-07 1.00 2.78e-08 1.00 1.00 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:2] "comp.1" "comp.2"
## $ all.loglik: num [1:7] -15574 -2260 -1955 -1955 -1955 ...
## $ restarts : num 0
## $ ft : chr "normalmixEM"
## - attr(*, "class")= chr "mixEM"
后验概率是一个$n \times k$的矩阵，是每个观测值由各个混合成分产生的概率，可以据此得到每个观测值的可能类别。
label <- c("head", "tail")[apply(em$posterior, 1, which.max)]
数据真正的标签在y中，可以得到混淆矩阵
xtabs( ~ y + label)
## label
## y head tail
## head 247 1
## tail 0 752
248个head，有1被标记为tail，标记错误。
可以进一步查看该观测值
x[y != label]
## [1] 4.390371
该值更靠近tail，所以高斯混合模型的结果判定为tail也不足为奇。
参考文献
[1]. http://rstudio-pubs-static.s3.amazonaws.com/1001_3177e85f5e4840be840c84452780db52.html
[2]. https://en.wikibooks.org/wiki/Data_Mining_Algorithms_In_R/Clustering/Expectation_Maximization_(EM)
