OpenCV4.1.1版本的EM算法实现

看了《模式识别与机器学习》，希望能够将其中的一些概念性较强的内容加以实现，来加强和巩固自己的理解

EM算法——《模式识别与机器学习 第九章》

1. 潜在变量的引入：使得复杂的概率分布可以有简单的分量组成。观测变量上的复杂的边缘概率分布可以通过观察变量与与潜在变量组成的扩展空间上的更加便于计算的联合概率分布来表示
2. EM算法——期望最大化：寻找带有潜在变量的模型的最⼤似然估计的⼀个⼀般的⽅法是期望最⼤化（EM）算法。其中，在许多应⽤中，参数由最⼤似然⽅法确定，通常会使⽤EM算法
   解释来说就是：E，求期望（Expectation），利用概率模型参数的现有估计值，计算潜在变量的期望；M，求极大（Maximization），利用E 步上求得的潜在变量的期望，对参数模型进行最大似然估计。所得参数估计值用于下个E步的计算，重复至收敛。
3. K均值聚类中的EM思想：其中，uk表示聚类的中心，xn表示D维欧几里得空间的随机变量x的N次观测结果。rnk表示二值指示变量，使用1-of-K表示方式，rnk表示数据点xn被分配到类别k。
   J为目标函数。最终的目标是找到{rnk}和{uk}的值，使得J达到最小值。
   在这个过程中，由于涉及到两个待优化变量，所以，分阶段对其进行最优化，同时在对一个变量进行优化时，保持另一个值固定不变。其中，更新rnk和更新uk的两个阶段分别对应于EM算法中的E步骤和M步骤。
   
   对于rnk的最优化：固定uj，即将聚类中心固定住，将每个观测值赋给最近的聚类中心
   
   uk的分母相当于聚类k中的数据点数目，此时uk即为类别k的所有数据点的均值
4. K均值算法：在实际应⽤中，⼀个更好的初始化步骤是将聚类中心选择为由K个随机数据点组成的子集。还有⼀点值得注意的地⽅，K均值算法本身经常被用于在EM算法之前初始化高斯混合模型的参数。
5. K均值算法的online random算法：使用顺序更新算法，对于每个数据点xn，使用下式更新最近的代表向量uk
6. k中心点算法：K均值算法中使用欧几里得距离作为数据点与代表向量之间的不相似程度的度量，此时限制了能够处理的数据变量的类型，也会使得中心点的确定对异常点不具有鲁棒性。此时，引入更加一般的不相似程度的度量
   此时，E步骤涉及到为每个数据点分配聚类，使得与对应的聚类代表的不相似程度最小。
   但是，M步骤通常比K均值的情形更加复杂，因此通常会将聚类原型 限制为等于 某个分配到那个聚类的数据向量，从而使得算法可以适⽤于任何不相似程度的度量V(.,.)。即，M步骤涉及到在分配到那个聚类的Nk个点上的离散搜索
7. 用于高斯混合模型的EM：
   高斯混合分布函数：其中pai k为混合系数
   
   那么对于每个观测数据点，都存在一个对应的潜在变量zn。该潜在变量为一个K为二值随机变量，采用1-of-K表示方法。根据哪个元素非零，向量z有K个可能的状态。
   
   将pai k看做zk=1的先验概率，那么上式即为观测到x之后，对应的后验概率。
8. 可以使用EM方法来解决由于奇异性的存在造成的应用于高斯混合模型的最大似然框架中的一个大问题
9. 高斯混合模型的对数似然函数
   令其关于高斯分量的均值uk的均值等于0， 则有
   整理得其中，Nk看作分配到聚类k的数据点的有效数量。 此时后验概率表示分量k对生成xn的责任
   同理，关于混合系数pai k最大化lnp，使用拉格朗日乘数法，
   得到即第k个分量的混合系数为那个分量对于解释数据点的责任的平均值
   
   
10. EM算法的另一种观点：在实际应用中，没有完整数据集{X, Z}，只有不完整的数据X。
    此时，适用于未观测的变量对应于数据集里的缺失值的清晰。观测值的概率分布可以通过对所有变量的联合概率分布关于缺失变量求和或积分的方式得到
11. 一般形式的EM算法：目标是最大化似然函数
    引入一个定义在潜在变量上的分布q(Z)，对于任意的q(Z)有
    
    其中，KL散度满足KL(q||p)>=0，当且仅当q(Z)=p(Z|X,θ)时等号成立
    在E步骤中，下界L(q, θ旧)关于q(Z)被最大化，而θ旧保持不变
    在M步骤中，分布q(Z)保持古典，下界L(q, θ)关于θ进行最大化，得到θ新

EM算法的OpenCV实现

1. EM算法属于非监督学习算法，可以从给定的样本集中，计算出高斯混合系数的最大似然估计，也能得到每个样本对应的标注值。
2. 关键是使用下述语句
   Ptr em = EM::create();
   em->setClustersNumber(3);
   em->setCovarianceMatrixType(EM::COV_MAT_SPHERICAL);
   em->setTermCriteria(TermCriteria(TermCriteria::COUNT+TermCriteria::EPS, 500, 0.1));
   em->trainEM(data, noArray(), modifiedImage, noArray());