混合模型初探
1. 混合模型简介
如果我们定义观测变量和潜在变量的一个联合概率分布，那么对应的观测变量本身的概率分布可以通过求边缘概率的方法得到。这使得观测变量上的复杂的边缘概率分布可以通过观测与潜在变量组成的扩展空间上的更加便于计算的联合概率分布来表示。
因此，潜在变量的引入使得复杂的概率分布可以由简单的分量组成。
2. 混合模型的一个特例 - KMean聚类
除了提供一个构建更复杂的概率分布的框架之外，混合模型也可以用于数据聚类，或者说聚类是混合模型一种特定形式。
我们从这小结开始，尝试从“数据聚类问题”切入对混合模型的讨论，我们来讨论数据点集合中的聚类的问题。
0x1：K均值聚类
我们首先使用一个非概率的方法解决这个问题，即K均值算法。
在我们后面谈到混合概率分布的潜在变量的时候，其中离散潜在变量可以被看作将数据点分配到了混合概率分布具体成分当中。K均值算法对应于用于高斯混合模型的EM算法的一个特定的非概率极限
1. 问题场景
我们考虑寻找多维空间中数据点的分组或聚类的问题。假设我们有一个数据集
，它由 D 维欧几里德空间中的随机变量 x 的 N 次观测（N次独立分布采样）组成。我们的目标是将数据集划分为 K 个类别。现阶段我们假定 K 的值是给定的。
2. 形式化描述K均值
直观上讲，我们会认为由一组数据点构成的一个聚类中，聚类内部点之间的距离应该小于数据点与聚类外部的点之间的距离。
为了形式化地说明这个概念，引入一组 D 维向量
，其中 k = 1，2，....，K，且
是与第 k 个聚类关联的一个代表，我们可以认为
表示了聚类的中心。我们的目标是找到数据点分别属于的哪一个
类，使得每个数据点和与它最近的向量
之间的误差评估函数（例如平方和）最小。我们定义一些记号来描述数据点的聚类情况，对于每个数据点
，我们引入一组对应的二值指示变量
，其中k = 1，2，....，K 表示数据点
属于 K 个聚类中的哪一个，从而如果数据点
被分配到类别 k，那么
=1，且对于 j != k，有
=0。这其实就是一种”“one-hot”编码方法。
之后我们可以定义一个目标函数，也即损失函数（loss function）：
它表示每个数据点与它被分配的向量
之间的距离的平方和。我们的目标是找到
和
的值，使得 J 达到最小值。
3. 求解目标函数极值的过程 - EM过程
我们可以用一种迭代的方法完成K均值的目标函数，其中每次迭代涉及到两个连续的步骤，分别对应于
和
的最优化。1）首先，我们为
选择一些初始值。即随机选择一些初始的质心。2）然后，在第一阶段，我们关于
最小化 J，保持
固定。对
求偏导我们就会发现，最优化公式可以转化为
。很显然，即在现有质心的前提下，将数据点归类为距离最近的聚类点，即使 J 最小。3）在第二阶段，我们关于
最小化 J，保持
固定。目标函数 J 是
的一个二次函数，直接求导求极值即可，
，可以很容易地解出
，即
。这个表达式有一个简单质朴的含义，即在现有聚类概率分布的前提下，重新计算各类中新的质心（取均值），这样即使得 J 最小。
4）不断重复这个二阶段优化直到收敛（直到聚类的分配不改变，或直到迭代次数超过了某个最大值）。
我们看到，更新
和更新
的两个阶段分别对应于 EM 算法中的 E（期望）步骤和 M（最大化）步骤。
由于每个阶段都减小了目标函数 J 的值，因此算法的收敛性得到了保证。
0x2：K均值改进 - K中心点算法（K-medoids algorithm）
K均值算法的基础是将平方欧几里德距离作为数据点与代表向量之间不相似度的度量，这不仅限制了能够处理的数据变量的类型（例如离散标签情况），而且使得聚类中心的确定对于异常点不具有鲁棒性。
我们可以这样推广K均值算法：引入连个向量 x 和 x′ 之间的一个更加一般的不相似度的度量
，然后最小化如下的失真度量（损失函数）：
与普通的K均值一样
1）对于给定的聚类代表
，E步骤涉及到为每个数据点分配聚类，使得与对应的聚类集中数据点间的不相似程序最小。这一步的计算代价为
，与标准的K均值算法的情形相同。2）对于不相似程度度量的一般选择，M步骤通常比K均值的情形更加复杂，因此通常会将聚类原型限制为等于某个分配到那个聚类的数据向量，因为这使得算法可以适用于任何不相似程度的度量
，只要它能够被计算。因此，对于每个聚类 k，M步骤涉及到在分配到那个聚类的
个点上的离散搜索，这需要
次对于
的计算。
0x3：硬划分与概率分布软划分
K均值算法的一个值得注意的特征是，在每次迭代中，每个数据点都被分配到一个唯一的聚类中，这友类类似经典决策树每层叶子的硬划分。虽然某些数据点与其聚类中心
的距离远小于其他中心的聚类，但是也存在其他的数据点，位于两个聚类中心的大概中间的位置。在后一种情形中，强行将数据点分配到最近的聚类不是最合适的。
类似于XGboost对决策树的改进（采用权值的方式代替决策树的硬分类），通过对K均值采用概率的方法，我们可以得到对数据点聚类的“软分配”，它反映了在最合适聚类分配上的不确定性，这个概率形式带来了一些数值计算上的优势。
3. 混合高斯
我们知道，高斯混合模型可以看成高斯分量的简单线性叠加，目标是提供一类比单独的高斯分布更强大的概率模型（多峰分布概率分布）。更进一步探求本质，我们可以使用离散潜在变量来描述高斯混合模型，这会让我们更加深刻地认识混合高斯模型，也会让我们了解期望最大化算法。
0x1：通过一个简单的“1-of-K”例子理解高斯混合分布的联合概率形式
让我们引入一个 K 为二值随机变量 z，这个变量采用了“1-of-K”表示方法，其中一个特定的元素
等于1，其余所有的元素都等于0。我们看到
的值满足概率分布的要求，即
，我们根据边缘概率分布
和条件概率分布
定义联合概率分布
。
z 的边缘概率分布根据混合系数
进行赋值，即：
，其中
以及
因此，给定 z 的一个特定的值（即选定一个特定的基模型），x 的条件概率分布是一个高斯分布：
联合概率分布为
，从而 x 的边缘概率分布可以通过将联合概率分布对所有可能的 z 求和的方式得到，即：
对于每个观测数据点
，存在一个对应的潜在变量
这非常关键，因为这意味着潜在变量不在是一个未知的变量，我们可以将其转化为对 x 的计算优化过程中（例如求导）。
得到这个联合概率分布的好处体现在几个方面：
1. 我们找到了高斯混合分布的一个等价的公式，可以将潜在变量 z 显示地写出
这么做的意义在于，我们现在能够对联合概率分布
操作。而不是对边缘概率分布
操作，这会产生极大的计算上的简化，后面我们将要讨论的期望最大化（EM）算法，将会看到这点。
2. 给定 x 的条件下，潜在变量 z 的条件概率可以用贝叶斯定义求出
我们可以用
表示
，它的值可以使用贝叶斯定理求出：
我们将
看成
=1的先验概率将
看成观测到 x 之后，对应的后验概率，也可以理解为分量 k 对于“解释”观测值 x 的“责任”（responsibility）
0x2：最大似然公式 - 高斯混合模型求解参数的基本策略
假设我们有一个观测的数据集
，我们希望使用混合高斯模型来对数据进行建模。我们可以将这个数据集表示为一个
的矩阵
，其中第 n 行为
。类似地，对应的隐含变量会被表示一个
的矩阵
，其中第 n 行表示为
。
如果我们假定数据点独立地从概率分布中抽取，那么我们可以使用下图所示的图模型来表示这个独立同分布数据集的高斯混合模型。
它对应的对数似然函数为：
最大化高斯混合模型的对数似然函数，比单一的高斯分布的情形更加复杂。困难的来源在于公式中，对 k 的求和出现在对数计算内部，从而对数函数不再直接作用于高斯分布。如果我们令对数似然函数的导数等于零，那么我们不会得到一个解析解，这给求导计算带来很大的困难。
一种方法是使用基于梯度的优化方法。虽然梯度的方法是可行的，并且梯度方法在混合密度网络的计算中起到了重要作用。但是我们这里讨论另一种主流的算法，被称为EM算法。它具有广泛的适用性，同时也是变分推断的基础。
0x3：用于高斯混合模型的EM
一种优雅并且强大的，用于寻找带有潜在变量的模型的最大似然解方法，被称为期望最大化算法（expectation-maximization algorithm），或者EM算法 。实际上，EM算法是可以推广到更一般化的变分推断框架中的，但是我们在这小结不着急，我们从高斯混合模型作为切入点，给出EM算法的一种非形式化的描述。
首先，让我们写下似然函数的最大值必须满足的条件，令公式
中
关于高斯分量的均值
的均值等于零，我们有：
值得注意的是，后验概率（或者说是“责任”）
很自然地出现在了等式右侧分母中。两侧同乘
（假设矩阵是非奇异的），整理，可得：
，其中
我们可以将
看作分配到聚类 k 的数据点的有效数量。仔细观察这个公式会发现，第 k 个高斯分量的均值
通过对数据集里所有的数据点求加权平均的方式得到，而这里所谓的权因子又是由后验概率
给出，而
表示分量 k 对生成
的责任（也可以理解为数据集分布的概率分布）。如果我们对
关于
的导数为零，然后用一个类似的推理过程，使用单一高斯分别协方差矩阵的最大似然结果，我们有：
这与一元高斯分布的对应的结果具有相同的函数形式，但是所不同的是，每个数据点都有一个自己的权值（不像K均值那样是硬划分），权值等于对应的后验概率，分母为与对应分量相关联的数据点的有效数量（分类集中的样本数量）。
最后，我们关于混合系数
最大化
，把限制条件考虑进去，即混合系数的累加和等于1。使用拉格朗日乘数法，最大化下面的量：
求导可得：
在右侧分母中，我们再次看到了“责任”
这一项。将两侧乘以
，对 k 求和，会发现
。使用这个结果消去
，整理可得：
。
一个非常质朴的表征含义：第 k 个分量的混合系数为那个分量对于解释数据点的“责任”的平均值。
值得注意的是，上式并没有给出混合模型参数的一个解析解，因为“责任”
在公式
中是以一种复杂的方式依赖这些参数。
然而，这些结果给出了一个简单的迭代方法来寻找问题的最大似然解。整个迭代过程是EM算法应用于高斯混合模型的一个实例。
1）我们首先为均值
2）然后，我们交替进行两个更新，被称为
2.1）E步骤（使用当前参数值计算责任）：在E步骤中，我们使用参数的当前计算公式 
2.2）M步骤（使用当前的责任重新估计参数）：然后在M步骤中，我们将E步骤计算出的概率用于最大化步骤，重新对模型参数（均值、协方差、混合系数）进行估值。 
。然后使用新的模型参数进入下一轮E步骤
3）计算对数似然函数： 
每次通过E步骤和接下来的M步骤对参数的更新确保了对数似然函数的增大，在实际应用中，当对数似然函数的变化量或者参数的变化量地域某个阈值时，我们就认为算法收敛。
上图给出了将两个⾼斯分布组成的混合概率分布的EM算法应⽤于⽼忠实间歇喷泉数据集的情形。这⾥，我们使⽤了两个⾼斯分布的混合（K=2）。
分布中⼼的初始值与图中的K均值算法使⽤了相同的初始值，精度矩阵被初始化为正⽐于单位矩阵。
图(a)⽤绿⾊标记出了数据点，以及初始的混合模型的配置，其中两个⾼斯分量的⼀个标准差位置的轮廓线分别⽤红⾊圆圈和蓝⾊圆圈标记。
图(b)给出了初始E步骤的结果，其中每个数据点的颜⾊中，蓝⾊所占的⽐重等于由蓝⾊分量⽣成对应数据点的后验概率，红⾊所占的⽐重等于由红⾊分量⽣成对应数据点的后验概率。因此，对于属于两个聚类的后验概率都较⼤的数据点来说，颜⾊看起来是紫⾊的。
图(c)给出了第⼀个M步骤之后的结果，其中蓝⾊⾼斯分布的均值被移⾄数据点的均值，同时根据属于蓝⾊类别的每个数据点的概率进⾏加权。换句话说，它被移到了蓝⾊标记数据点的质⼼。类似地，蓝⾊⾼斯分布的协⽅差被设置为蓝⾊标记数据点的协⽅差。红⾊分量的情形与此类似。图(d),(e)和(f)分别给出了2次、5次、20次完整的EM循环之后的结果。在图(f)中，算法接近收敛。
注意：
1. 与K均值算法相比，EM算法在达到（近似）收敛之前，经历了更多次的迭代，每个迭代需要更多的计算量。因此，通常运行K均值算法找到高斯混合模型的一个合适的初始化值，输入EM，接下来使用EM算法进行微调节。
2. 协方差矩阵可以初始化为通过K均值算法找到的聚类的样本协方差。
3. 混合系数可以被设置为分配到对应类别中的数据点所占的比例。
4. EM的另一种观点 - 泛化讨论
我们继续讨论EM算法的另⼀种观点，其中潜在变量起着重要的作。我们⾸先使⽤⼀种抽象的⽅式讨论这种⽅法，然后我们再次考虑⾼斯混合模型的例⼦，来具体说明这个模型。
EM算法的目标是找到具有潜在变量的模型的最大似然解。
我们将所有观测数据的集合记作
，其中第 n 行表示
，所有潜在变量的集合记作
，对应的行为
。所有模型参数的集合被记作
。因此对数似然函数为：
，该式同样适用于连续潜在变量的情形，只需要对
的求和替换为积分即可。一个关键的现象是，对于潜在变量的求和位于对数的内部。即使联合概率分布
属于指数族分布，但是由于这个求和式的存在，求和得到的结果边缘概率分布
通常也不是指数族分布。求和式的出现阻止了对数运算直接作用于联合概率分布，使得最大似然接的形式更加复杂。现在假设对于
中的每个观测，我们都有潜在变量
的对应值。我们将
称为完整（completely）数据集，并且我们称实际的观测数据集
是不完整的（incomplete）。完整数据集的对数似然函数的形式为
，并且我们假定对这个完整数据的对数似然函数进行最大化是很容易的。然而遗憾的是，在实际应用中，我们没有完整数据集
。我们关于潜在变量
的取值的知识仅仅来源于后验概率分布
。1）由于我们不能使用完整数据的对数似然函数，因此我们反过来考虑使用当前的参数值
，在潜在变量的后验概率分布
下，使用这个后验概率分布计算完整数据对数似然函数对于一般的参数值
的期望。这个期望被记作
，这对应于EM算法中的E步骤：
2）在接下来的M步骤中，我们最大化该式：
来确定修正后的参数估计
。如果当前对于参数的估计为
，那么一次连续的E步骤和M步骤会产生一个修正的估计
。注意，在
的定义中，对数操作直接作用于联合概率分布
，因此根据假设，对应的M步骤的最大化是可以计算的。
每个EM循环都会增大不完整数据的对数似然函数（除非已经达到局部极大值），EM的迭代过程是由数据集本身的概率分布驱动的，如果数据集本身存在局部最优区域，EM算法就有可能在EM过程中落入该区域。
0x1：重新考察高斯混合模型
我们继续考虑将EM算法的潜在变量的观点应用于高斯混合模型例子。我们知道，在高斯混合模型的求解中，我们的目标是最大化对数似然函数：
它是使用观测数据集
进行计算的，我们看到这个计算比单一高斯分布的情形更加困难，因为对 k 的求和出现在对数运算内部。现在考虑对完整数据
进行最大化。根据公式
和公式
，似然函数的形式为：
，其中
表示
的第 k 个分量。取对数，我们有：
（2）
（2）公式和（1）公式对比，我们看到在 k 上的求和与对数运算的顺序交换了。对数运算现在直接作用于高斯分布上，而高斯分布本身是指数族分布的一个成员，这个方法产生了最大似然问题的一个简单得多的解。
因此我们看到，完整数据的对数似然函数可以用一种简单的方法求出最大值的解析解。然而，在实际应用中，我们并没有潜在变量的值，因此，与之前讨论的一样，我们考虑完整数据对数似然函数关于潜在变量后验概率分布的期望，形式为：
。因此后验概率分布可以在 n 上进行分解，从而
是独立的，指示值
的期望为：
它就是 k 分量对于数据点
的“责任”，于是，完整数据的对数似然函数的期望值为：
，接下来可以使用EM步骤进行迭代求解。在后面的讨论中，我们会继续深入讨论完整数据的对数似然函数的期望的作用。
0x2：与K均值的关系
高斯模型的EM算法和K均值算法，在本质上有很强的相似性。
1. K均值算法对数据点的聚类进行了“硬”分配，即每个数据点只属于唯一的聚类。
2. 而EM算法基于后验概率分布，进行了一个“软”分配
实际上，我们可以将K均值算法看成高斯混合模型的EM算法的一个特殊的极限情况。下面我们来讨论这种情况：
考虑一个稍微特殊的高斯混合模型，其中混合分量的协方差矩阵为
，
是一个被所有分量共享的方差参数，
是单位矩阵，从而：
。我们现在考虑K个这种形式的高斯分布组成的混合模型的EM算法，其中我们将
看作一个固定的常数，而不是一个需要重新评估的参数。对一个特定的数据点
，后验概率（“责任”）为：
。如果我们考虑极限情况
（所有的数据点都有且只属于一个分类，不存在概率分布，则方差为0），那么我们看到，在分母中，
最小的项将会慢慢地趋近于零，因此对于数据点
，只有项 j 的“责任”
趋近于1，其他的项的责任都趋近于0.因此，在这种极限情况下，我们得到对数据点聚类的一个硬分类，与K均值算法相同，从而
，因此，每个数据点都被分配为距离最近的均值的聚类。这样，EM重估计就简化为了K均值的结果。最后，在极限
的情况下，公式
给出的完整数据的对数似然函数变成了：
。
因此在极限的情况下，最大化完整对数似然函数的期望等价于最小化K均值算法的失真度量J。
0x3：与伯努利分布的混合
我们之前的讨论集中在由混合高斯模型描述的连续变量的概率分布上，为了更加泛化的讨论混合模型（混合模型不依赖具体的概率分布先验假设），我们现在讨论由伯努利分布描述的离散二值变量的混合，这个模型也被称为潜在分布（latent class analysis）。这个模型不仅具有实际应用的价值，还是我们讨论离散变量上的马尔科夫模型的基础。
考虑D个二值变量
组成的集合，其中
，每个变量都由一个参数为
的伯努利分布控制，即：
，其中，
，且
。我们看到，在给定
的条件下，各个变量
是独立的。因此，在独立同分布情况下，这个分布的均值和方差可以很容易求得：
。我们先来尝试讨论这种分布的有限混合，即：
，其中，
，且
。这个混合分布的均值和方差为：
，其中
。
由于协方差矩阵 
如果我们有一个数据集
，那么这个模型的对数似然函数为：
。
和在高斯混合模型中一开始的情况一样，我们看到求和运算位于对数运算内部，从而最大似然解没有解析解。
为了解决这个问题，我们需要引入混合伯努利分布的最大化似然函数的EM算法。
1. 混合伯努利分布的最大化似然函数的EM算法推导
我们首先引入一个潜在变量 z，它与 x 的每个实例相关联。与高斯混合模型的情形相同，
是一个二值 K 维变量，其中只有一个元素等于1，其余元素等于0。这样，给定潜在变量，我们可以写出 x 的条件概率分布，形式为：
。而潜在变量的先验概率分布与高斯混合模型的形式相同，即：
。如果我们将
和
相乘，然后对 z 求和，我们同样得到该公式：
为了推导EM算法，我们首先写出完整数据的对数似然函数，形式为：
，其中
。
接下来我们取完整数据对数似然函数关于潜在变量后验概率分布的期望，得：
，其中
是给定数据点
的条件下，分量 k 的后验概率分布，或者“责任”。
在E步骤中，这些后验概率使用贝叶斯定理计算，形式为：
在上面关于潜在变量后验概率分布的期望公式中，我们可以看到，如果我们对 n 求和，“责任”只出现在两项中，这两项可以写成：
，其中
是与分量 k 的均值组成的集合等于数据的加权平均值，权系数为分量 k 对于数据点的“责任”。对于关于
的最大化，我们需要引入一个拉格朗日乘数来满足限制条件
。采用与高斯混合模型中类似的步骤，我们有：
，即分量 k 的混合系数等于数据集里那个分量的数据点所占的比例，
伯努利分布参数的共轭先验是Beta分布。我们已经看到一个Beta先验分布等价于引入 x 的额外的有效观测。类似地，我们可以引入伯努利混合模型的先验分布，然后使用EM算法最大化后验概率分布。
Relevant Link:
http://sklearn.apachecn.org/cn/stable/modules/mixture.html
5. 一般形式的EM算法
期望最大化算法，或者EM算法，是寻找具有潜在变量的概率模型的最大似然解的一种通用的方法。
这里，我们考虑一个概率模型，其中我们将所有的观测变量联合起来记作
，将所有的隐含变量记作
。联合概率分布
由一组参数控制，记作
。我们的目标是最大化似然函数：
这里，我们假设
是离散的，但是当
是连续变量或者离散变量与连续变量的组合时，方法是完全相同的，只需要把求和换成适当的积分即可。我们假设直接优化
比较困难，但是最优化完整数据似然函数
就容易得多。接下来，我们引入一个定义在潜在变量上的分布
。我们观察到，对于任意的
，下面的分解成立：
注意，
是概率分布
的一个泛函，并且是参数
的一个函数。根据上式，我们看到
是
和后验概率分布
之间的Kullback-Leibler散度。我们知道，Kullback-Leibler散度满足
，当且仅当
时等号成立。因此，根据公式
我们得
。换句话说，
是
的一个下界，如下图：
0x1：从KL散度角度看EM一般化过程
EM算法是一个两阶段的迭代优化算法，用于寻找最大似然解。我们可以使用公式
来定义EM算法，证明它确实最大化了对数似然函数。假设参数向量的当前值为
在E步骤中，下界
关于
被最大化，而
保持固定。最大化问题的解很容易看出来。我们注意到
不依赖于
，因此
的最大值出现在Kullback-Leibler散度等于0的时候，换句话说，最大值出现在
与后验概率分布
相等的时候。此时，下界等于对数似然函数，如下图所示：
在接下来的M步骤中，分布
保持固定，下界
关于
进行最大化，得到了某个新值
。这会使得下界
增大（除非已经达到了极大值），这会使得对应的对数似然喊增大。由于概率分布 q 由旧的参数值确定，并且在M步骤中保持固定，因为它不会等于新的后验概率分布
，从而KL散度非零。于是，对数似然函数的增大量大于下界的增加量，如下图所示：
此时，下界的形式为：
其中，常数就是分布 q 的熵，因此与
无关。从而在M步骤中，最大化的量是完整数据对数似然函数的期望，正如我们在之前的混合高斯模型的情形中看到的一样。
0x2：从参数空间角度看EM一般化过程
EM算法的计算也可以被看作参数空间中的运算，如下图所示：
红色曲线表示（不完整数据）对数似然函数，它的最大值是我们想要得到的。我们首先选择某个初始的参数值
，然后在第一个E步骤中，我们计算潜在变量上的后验概率分布，得到
的一个更小的下界，它的值等于在
处的对数似然函数值，用蓝色曲线表示。注意，下界与对数似然函数在
处以切线的方式连接，因此两条曲线的梯度相同。这个界是一个凹函数，对于指数族分布的混合分布来说，有唯一的最大值。在M步骤中，下界被最大化，得到新的值
，这个值给出了比
处更大的对数似然函数只。接下来的E步骤构建了一个新的下界，它在
处与对数似然函数切线连接，用绿色线白哦是。如此循环往复。