EM算法对鸢尾花无监督分类学习

引子

  EM算法，也叫最大期望算法，是一种聚类算法，常用于无监督学习。我认为EM算法是一种很神奇的算法，比如说，你做了五轮抛硬币的实验，每轮抛10次，实验结果如下

图a

  如果告诉你每轮实验都是用A硬币或者B硬币抛的，但是每一轮实验到底是用A抛的，还是B抛的就不知道了。请你计算A硬币抛出正面的概率。

  事实上，这个问题并不好解答。但是如果实验结果是这样呢

图b

  补全了表之后，可知用A硬币抛了20次，出现9次正面，那么A硬币抛出正面的概率就是$9/20$。

  通过这个例子我们发现，要解决刚开始提出的那个问题，只需要知道每轮投掷的硬币即可。怎么才能知道呢？

  我们似乎忽略了一个已知条件，每一次实验不管选择的是哪一枚硬币，出现正面的次数都服从二项分布，而且N=10。不同的是：如果使用A硬币，出现正面次数服从$B(10,p_{A正})$，如果使用B硬币，出现正面次数服从$B(10,p_{B正})$，其中$p_{A正}、p_{B正}$分别代表用硬币A、B抛出正面的概率，即模型的未知参数。

  我们不妨先初始化模型的未知参数$p_{A正}=0.4，p_{B正}=0.9$，这样一来，A/B硬币的正面次数分布已知，计算产生第一轮实验结果的概率就很容易了：

$$\begin{align*} &用硬币A抛：p(抛出5次正面)=C_{10}^{5}\cdot0.4^5\cdot0.6^5=0.20066 \\ &用硬币B抛：p(抛出5次正面)=C_{10}^{5}\cdot0.9^5\cdot0.1^5=0.001488 \end{align*} \tag{1.1}$$

  容易看的出来，用硬币A产生第一轮实验结果的概率远远大于使用硬币B。我们就认为第一轮实验更有可能使用的是硬币A。同样的我们可以计算出剩余四轮实验更有可能使用的硬币。这样每一轮实验投掷硬币的期望(是硬币A还是硬币B)我们就求出来了，这也是EM算法中的E步。

  当然，这还未结束。别忘了我们之所以能求出这个期望，初始化的参数功不可没。如今我们有了每一轮的类别期望后，就把图a的问题转换成图b的问题了，通过重新计算$p_{A正}$和$p_{B正}$来更新模型的参数，这一步就是M步。之后根据新的参数再来计算每一轮实验的期望，循环往复，直到参数不再发生改变或者到达一定的迭代次数。这便是EM算法的思想。可以总结为下图

  真实的EM算法要比我讲的更复杂一点，主要复杂在E步。我们原先认为每一轮投掷硬币的期望不是A就是B，这种非此即彼的想法还不能准确描述真正的期望。事实上，通过梳理1.1式我们还能得到更有趣的结果

$$\begin{align*} &有\frac{0.20066}{0.20066+0.001488} = 0.992639的概率选择硬币A \\ &有\frac{0.20066}{0.20066+0.001488} = 0.007361的概率选择硬币B \tag{1.2} \end{align*}$$

  用和为1的两个概率来作为两个类别的权重，使得期望更准确，随后更新的参数自然也更快收敛。

  EM算法的思想大概就是这样，至于为什么EM算法通过迭代会逼近最优解，即算法收敛，将在理论部分展开推导。

实战

  本次实战仅使用鸢尾花数据集的特征部分，而标签部分则在最后评价模型时用作参考。

  同样的，我们先从sklearn的datasets库中导入数据集，并将特征数组转换为DataFrame类型。

from sklearn import datasets
import pandas as pd
ds = datasets.load_iris()
df_ds = pd.DataFrame(ds.data, columns=ds.feature_names)

  EM算法是聚类算法，常用的模型是GMM高斯混合模型。sklearn中的高斯混合模型创建时有以下参数

n_components : 高斯混合模型的个数，也就是我们要聚类的个数，默认值为1
convariance_type : 协方差类型。一个高斯混合模型的分布是由均值向量和协方差矩阵决定的，所以协方差的类型也代表了不同的高斯混合模型的特征。
   full类型，每个分量有各自不同的完全协方差矩阵，即元素都不为0
   tied类型，所有分量有相同的完全协方差矩阵
   diag类型，每个分量有各自不同的对角协方差矩阵，即除对角线外所有元素为0
   spherical类型，每个分量有各自不同的球面协方差矩阵，即对角线元素相同，其余元素为0
max_iter : 代表最大迭代次数，EM算法是由E步和M步迭代求得最终的模型参数，默认值为100

  在这里，我们将鸢尾花分成三类。调用fit()函数训练，调用predict()函数预测。并查看预测的分类和实际的标签分类是不是一致(EM算法聚类后打上的标签序号可能和实际标签序号不一致，可多尝试一下)

from sklearn.mixture import GaussianMixture
from sklearn.metrics import calinski_harabasz_score
gmm = GaussianMixture(n_components=3) #分为三类
gmm.fit(df_ds)
label=gmm.predict(df_ds)
print("准确率:{:.2%}".format(len(label[label==ds.target])/len(label)))
print("得分:{:}".format(calinski_harabasz_score(df_ds, label)))

  多次训练，找到标签序号一致的预测结果，输出如下

准确率:96.67%
得分:481.78070899745234

  将预测的分类标签，以及真实标签画成散点图

import matplotlib.pyplot as plt
plt.scatter([i for i in range(len(ds.target))], ds.target, color='black', marker='v')
plt.scatter([i for i in range(len(label))], label, color='red', marker='.')
plt.yticks([0, 1, 2])
plt.xlabel('samples')
plt.ylabel('category')
plt.legend(['actual', 'predict'])
plt.savefig('predict.jpg', dpi=800)

  输出如下

  大约150个数据点中只有5个数据点是分类错误的，正确率高达96%以上。效果还是比较好的。

理论