聚类之EM算法

定义：在统计计算中，最大期望（EM）算法是在概率模型中（E步）寻找参数最大似然估计或者最大后验估计（M步）的算法，其中概率模型依赖于无法观测的隐藏变量（LatentVariable），聚类应用中就是隐藏的类别。

期望最大算法是一种从不完全数据或有数据丢失的数据集（存在隐含变量）中求解概率模型参数的最大似然估计方法。

该算法既简单又复杂，简单的是其仅包含了两个步骤（E步和M步）就能完成强大的功能，复杂的是数学推理涉及到比较繁杂的概率公式。下面先举一个简单的例子也许就能大致明白这里的E步和M步了。

如下图所示，就是一个完整数据的参数估计：

假设某研究人员做了5次试验，每次选择A硬币或B硬币进行10次投掷，最后记录5次试验中硬币正面朝上的次数。那么问题来了，A硬币正面朝上的概率与B硬币正面朝上的概率一样吗？常识告诉我们，两种硬币正面朝上的概率均为0.5，事实如何呢？其实，根据试验我们是可以计算出A硬币和B硬币正面朝上的概率的，最简单粗暴的计算方法是：

P(A)=硬币A正面朝上的次数/硬币A投掷的次数=(5+7+5)/30=0.57

P(B)=硬币B正面朝上的次数/硬币B投掷的次数=(3+8)/20=0.55

还可以使用最大似然法计算出P(A)和P(B)。具体步骤如下：

1）由于只有两种硬币，可以假设选择硬币的分布服从二项分布，二项分布的概率公式为：

2）根据二项分布，构造5次试验的似然函数：

3）对含有参数P(A)和P(B)的似然函数求偏导，并令偏导函数为0，可进一步计算出P(A)和P(B)的估计值。

如果还是这个试验，但是收集到的数据只是正面朝上的次数，并不知道是投掷的哪一种硬币，该如何求解这5个试验分别使用哪种硬币呢？正如下图所示就是一种不完整数据的参数估计：

对于这种问题其实就是一种聚类，即5次试验该聚为哪两类？可以通过EM算法将5次试验进行聚类，其聚类的思路需要解决两个问题：

1）每次试验使用A硬币或B硬币的概率是多少？

2）根据计算出来的概率，判别A或B，由此得到的新的参数P(A)和P(B)值又是多少？

而这两个问题又是互相依赖的，即知道其中一个就可以算的另一个。为了更清楚的理解EM算法，不妨做一次EM算法的步骤：

1）假定初始的P(A)和P(B)值，不妨P(A)=0.6，P(B)=0.5。

2）对于第一次试验，A硬币或B硬币出现5次正面朝上的概率，该步为E步

从而断定，第一次可能投掷B硬币的概率比较大，依次可以计算后四次用某种硬币投掷的可能性，不妨为如下结果：

3）由第一次的E步可以实现完整数据的参数估计，接下来就要计算似然函数和参数估计，该步为M步

进一步可以根据新的参数估计P(A)和P(B)再次计算E步，这样不断的迭代优化，最终实现收敛完成聚类。

上面的过程只是易于理解EM算法的一个简单例子，有关理论部分的知识课参考后文的链接。

应用

R语言中提供了EM算法的软件包，通过使用mclust包中的Mclust()函数可以实现EM聚类，首先看一下该函数的参数解释：

Mclust(data, G = NULL, modelNames = NULL,

prior = NULL,

control = emControl(),

initialization = NULL,

warn = mclust.options("warn"), ...)

data提供用于分析的数据集，切记数据集中不可有分类变量；

G指定聚类的数目，默认为1:9类；

modelNames指定EM聚类过程中的拟合模型；

prior指定先验值，默认不指定；

control指定EM算法的控制参数，如收敛阈值、最大迭代次数等

initialization为算法指定初始值；

warn逻辑值，是否反馈某些警告信息。

这里使用R重自带的iris数据集进行EM算法的聚类：

#加载mclust包

library(mclust)

#创建训练数据集和测试数据集

index <- sample(1:2, size = nrow(iris), replace = TRUE, prob = c(0.7,0.3))

train <- iris[index == 1,]

test <- iris[index == 2,]

#EM算法进行聚类

EM <- Mclust(data = train[,-5], G=3)

summary(EM)

#判断样本内的预测准确率

Freq1 <- table(EM$classification, train[,5])

accuracy1 <- sum(diag(Freq1))/sum(Freq1)

accuracy1

#判断样本外的预测准确率

pred <- predict(object = EM, newdata = test[,-5])

Freq2 <- table(pred$classification, test[,5])

accuracy2 <- sum(diag(Freq2))/sum(Freq2)

accuracy2

算法最大的缺陷是数据量越大时，算法的效率会越低。

参考资料

http://blog.csdn.net/crazyhacking/article/details/18982587

http://wenku.baidu.com/link?url=y80Pc7JQcw_3uTJBr6VCB1m9uLUADrT-Xg3vQD80OxZNjB9RzkSWcrKyIiGSXL2uVZN1MnMZg1yaKMMFGUqTc0w_6qThP9b4ml7VbnK--SC

http://www.cnblogs.com/zhangchaoyang/articles/2623364.html

http://www.tuicool.com/articles/Av6NVzy

数据挖掘：R语言实战