统计学习方法—EM算法及其推广

一、相关概念

1.先验概率

在贝叶斯统计学中,某个未知参数p的概率先验分布是在不考虑观察到的数据的情况下,所表示对p不确定性程度的概率分布。这种分布主要用于描述该变量本身的不确定性而非其随机性质,既可以依据历史观察到的数据进行统计分析,也可以依据领域内的背景知识确定,还可以由个人主观判断提供。通常都是用来表示单次事件的概率,例如P(x)、P(y)等。

2.后验概率

在贝叶斯统计中，在考虑到并提供了相关证据或数据之后所获得的是条件概率 。同样地，在试验与调查的基础上得出的概率分布即为后验概率分布；或者说是根据先验概率所获得的逆向条件概率，并与条件概率的形式相一致。

    P(y|x) = ( P(x|y) * P(y) ) / P(x)

其中：

    P(y) 是先验概率，一般都是人主观给出的。贝叶斯中的先验概率一般特指它。
    
    P(x|y) 是条件概率，又叫似然概率，一般是通过历史数据统计得到。一般不把它叫做先验概率，但从定义上也符合先验定义。
    
    P(y|x) 是后验概率，一般是我们求解的目标。（已知x的前提下，然后求y发生的概率）
    
    P(x) 其实也是先验概率，只是在贝叶斯的很多应用中不重要（因为只要最大后验不求绝对值），需要时往往用全概率公式计算得到。

3. 似然概率

通常所说的概率是指在给定参数的情况下预测即将发生的事件的可能性。拿抛掷均匀硬币的例子来说,已知其正反面出现的概率均为0.5,则想知道两次抛掷结果均为正面的概率是多少

    H代表Head，表示头朝上
    
    p(HH | pH = 0.5) = 0.5*0.5 = 0.25.

这种表达方式可能带来误导性信息。其中这个p实际上是作为参数存在的，在其值上相当于一个未知的参数（其值可取如p=0.4/0.5/0.6等），它并不是一个随机变量，在传统条件下不能将其表示为条件概率的形式，在传统的条件概率框架中都是基于随机变量的概念。因此更为稳妥的做法应该是写作 p(HH;p=0.5)。

(2) 似然概率

这两个事件的概率被称为似然函数；也就是基于观察到的结果推测某个参数（p）的可能性；这相当于将其视为给定条件下计算概率的方式。

    **L**(pH=0.5|HH) =P(pH=0.5|HH) = P(pH=0.5,HH)/P(HH)；既然HH是已经发生的事件，理所当然P(HH) = 1

    所以：**P(pH=0.5|HH)  = P(pH=0.5,HH) = P(HH;pH=0.5)**

所以，我们可以safely得到:

    L(pH=0.5|HH) = P(HH|pH=0.5) = 0.5*0.5=0.25.

这个0.25的意思是，在已知抛出两个正面的情况下，pH = 0.5的概率等于0.25。

再算一下

    L(pH=0.6|HH) = P(HH|pH=0.6) = 0.6*0.6=0.36.

…

把pH从0~1的取值所得到的似然函数的曲线画出来得到这样一张图：

注意到当pH = 1时似然函数（概率）达到最大值。
使似然函数达到最大值的参数是p=1。
这样就容易理解了。
最后，在观察到两次抛掷结果均为正面（HH）的情况下，在这种情况下，默认假设硬币在未抛掷前正面朝上的概率为1。
显然这个结果可能不符合真实情况（数据量过少所致）。

4. 最大似然估计

具体而言，在实际应用中我们通常会遇到这样的情形：当我们面对一组随机采样的观测数据时，并不清楚其背后的分布规律或相关参数的真实取值范围；此时就需要借助统计推断的方法来建立合理的概率模型并求解相应的未知参数。”

基于这种基本原理，在最大似然估计理论框架下我们假设存在某个特定的参数θ能够使得观测数据X出现的概率密度函数P(X|θ)达到全局最大值；基于这一核心假说，在缺乏其他有效信息的情况下我们有理由相信θ即为待估参数的真实值。”

抽取到样本的概率 ，也就是样本集X中各个样本的联合概率，用下式表示：

这个概率表征了，在给定参数θ的情况下，观察到样本集X的概率。由于我们已经观测到样本集X，并且参数θ尚未确定，则该公式仅包含待估计的参数θ。因此它成为了一个关于θ的函数。这便定义了参数θ相对于样本集X的似然函数（likelihood function），通常表示为L(θ)。

通常情况下，我们能够观察到L(θ)是由多个因子相乘构成的。为了简化计算过程，在分析过程中我们可能会定义对数似然函数，并将其转换为累加的形式。

现在我们明确了，在统计学中为了估计参数θ的最大似然估计量θMLE的过程中，我们需要最大化对数似然函数L(θ)。这其实就是寻求该函数的最大值问题。那么问题来了：求取一个函数的最大值或最小值的方法是什么？通常的做法是计算导数，并令其导数等于零。这一步骤会得到一系列方程组的结果，在连续可微的情况下这些方程组就可以被解出来从而得到对应的参数估计值θMLE。当参数向量维度为n时我们需要计算所有元素关于各个变量的一阶偏导数组成的梯度向量∇_θ L(θ)，进而建立n个联立方程组并求解这些方程组就可以找到使得对数似然函数取得极大值或者极小值的位置点（当然这里我们关注的是极大情况）。因此在实际操作中无论是单参数还是多参数的情形都需要通过这种方法来进行推断和计算得出最终的结果

求最大似然函数估计值的一般步骤：

（1）写出似然函数；

（2）对似然函数取对数，并整理；

（3）求导数，令导数为0，得到似然方程；

（4）解似然方程，得到的参数即为所求；

5.如何求解边缘概率分布

考虑二维离散型随机向量(X, Y)，其中X与Y为两个离散型随机变量，并令{Xi}及{Yi}分别代表X与Y的所有可能取值集合，则其联合概率分布P(X=xi, Y=yi)可用以下函数表达式或表格形式来描述它们的概率分布情况：

则只包含部分变量的概率分布称为边缘分布为：

二、EM算法

1.引入实例

（9.1）式公式推导：

由以上三个等式可以得到：

该公式基于单次实验设计，在这种情况下

解释：为何第三个等式的左侧中没有出现变量Z？这是因为对变量Z采用了求和运算（sum），即等价于对所有可能的Z值进行了求和操作，并因此无需单独写出变量Z项。这种处理方式表明，在统计建模中可以通过观测数据y以及隐藏随机变量z之间的关系来推导出相应的似然函数表达式。

好了, 现在我们清楚了, 估计θ的方法是通过最大化θ的似然函数L(θ), 然后使得其最大值对应的参数即为我们所需的估计

解释1：
该式子就是似然函数，现在要求解参数

，即即可通过最极大似然估计，求最优参数

，使得以上条件概率最大（或者对数条件概率最大）。即：

解释2：
此处通过取似然函数对数的形式来最大化。

该问题不具备解析解（闭式解），仅能借助迭代方法进行求解。EM算法即该算法的一种迭代方法以解决该问题。

2.EM算法推导

EM算法通过迭代求

的极大似然估计，在每次迭代的过程中会被划分为两个步骤：首先进行期望步骤（Expectation step），计算条件期望；接着进行最大化步骤（Maximization step），求取参数的最大值以优化目标函数）。随后将详细阐述EM算法的推导过程

假设我们有一个样本集{x(1),…,x(m)}，包含m个独立的样本。但每个样本i对应的类别z(i)是未知的（相当于聚类），也即隐含变量。故我们需要估计概率模型p(x,z)的参数θ，但是由于里面包含隐含变量z，所以很难用最大似然求解，但如果z知道了，那我们就很容易求解了。
对于参数估计，我们本质上还是想获得一个使得似然函数最大化的那个参数θ，现在与最大化似然不同的只是似然函数式中多了一个未知的变量z。也就是说我们的目标是：找到合适的θ和z让L(θ)最大 。那我们也许会想，你就是多了一个未知的变量而已啊，我们也可以分别对未知的变量θ和z分别求偏导，然后再令其为0，求解出来不也一样吗？
现将一个样本(里面含有很多子样)的似然函数

扩展至样本集m上，并将其中包含的m个样本的似然函数进行求和运算，请参考公式（1）。

解释：

对（1）式而言，在离散情况下掌握联合分布情况的基础上推导边缘分布函数（这一操作中变量z同样扮演随机变量的角色）。具体来说，在计算过程中我们需针对每个样本i的所有可能类别z计算右边的联合分布之和；这个结果即为我们定义的似然函数值（虽然在实际操作中会遇到将"和"与"对数"结合使用的情况）。值得注意的是，在这一过程中我们遇到了一些复杂性——特别是在处理包含多个项相加后再取对数的操作时（想象一下log(f1(x)+ f2(x)+ f3(x)+...)这种复合函数的形式），这使得直接求解变得困难重重。因此我们需要另辟蹊径寻找解决方案。
观察到（2）式的构造方式后我们会发现其实它并未带来本质性的改进——仍然存在"和的对数"这一关键问题；因此这种方法也难以奏效。或许我们应该考虑其他途径。
转而研究（3）式我们会发现它实际上实现了另一种转换——将原本复杂的"对数运算"转换为了更为简单的"对数相加"形式；这样一来整个运算过程就变得简单明了许多——这正是Jensen不等式的魅力所在之处。

jensen不等式

Jensen 不等式的定义可表述为：
（1）若考虑 f 为凸函数且 X 作为随机变量，则有 E[f(X)] ≥ f(E[X])。
（2）特别地 ，若 f 为严格凸函数，则仅在 X 为常量时方程两边相等。
此外，
对于凹函数的情形而言 ，Jensen 不等式的应用会导致其方向与通常情况相反。
具体而言，
这表明期望值的变换与输入值的变换之间存在逆向关系。

解释：

图标说明：
曲线f表示一个凸函数,其中随机变量X以相等的概率取值a或b,因此其期望值\mathbb{E}[X]必然位于a与b之间。相应地,\mathbb{E}[f(X)]作为f(X)的整体期望值,必然位于f(a)与f(b)之间。

回到公式（2），因为f(x)=log x为严格凹函数，（其二次导数为-1/x^2<0）。

上式子正好是

的期望（考虑到

，得到：

，又因为：

，所以就可以得到公式（3）的不等式了，即通过凹函数f(E[X]) >= E[f(X)]得出。

（Q是随机变量z(i)的概率密度函数 ）
现在式（3）就容易地求导了，但是式（2）和式（3）之间是不等号啊，要求解式子（2）的最大值，该如何利用式（3）来求解式（2）的最大值呢？现在我们就需要一点想象力了，上面的式（2）和式（3）不等式可以写成：似然函数L(θ)>=J(z,Q)，那么我们可以通过不断的最大化这个下界J，来使得L(θ)不断提高，最终达到它的最大值。即不断最大化（3）式来不断使式子（2）逼近它的最大值。
方法解释图如下：

解释：

对于不等式 L(\theta) \geq J(z, Q)
我们首先固定\theta_1
随后寻找使得 L(\theta_1) = J(z, Q)成立的最优Q
接着逐步优化J
选取优化后的新参数值\theta_2
同时保证 L(\theta_2)\geq J(z, Q)仍然成立
并令\theta_1 = \theta_2
重复上述步骤以获得新的参数值\theta_3
同时保证 L(\theta_3)\geq J(z, Q)
持续下去，则L将逐步增大直至收敛到其最大可能值

这里有两个问题：什么时候下界J(z,Q)与L(θ)在此点θ处相等？为什么一定会收敛？

（1）固定θ值并使其在该点θ处与L(θ)相等

首先对分母分子关于z求和，因为分母求和为，

，所以分子求和为：

所以可得：

即得到：

转化上述公式：

所以，在固定参数θ之后，在满足上述方程的前提下实现了对Q(z)的选择问题的解决，并使得即使在提升下界的情况下计算得到的Q(z)仍然等于其后验概率分布。这一步骤即为E步，在此过程中我们建立了L(θ)的一个下界表达式

计算公式：

在给定Q(z)的情况下进行M步迭代，在这一步中我们需要通过调整参数θ来最大化L(θ)的下界J（即使固定了Q(z)，我们仍然有机会进一步提高J）。通过最大化似然函数来更新模型参数。

该优化过程通过逐步迭代的方式实现对似然函数L(θ)的最大化求解，在此过程中我们深入探讨其收敛性以及参数值θ的变化规律。
（2）为了深入探讨其收敛性，请问L(θ)最大化的过程中是否会趋于稳定？换句话说，在迭代过程中参数值θ是否会趋向于一个稳定的极限点？
直观上来看，在每一次迭代中下界都会持续提升这一事实表明极大似然估计具有严格的单调递增特性，并由此可推出该估计量必定存在一个有限的极限点作为全局最大值的位置。

理性分析的话，就会得到下面的东西：

请提供详细说明证明过程，并参考推导过程

3.EM算法举例说明

（1）利用坐标上升法理解EM算法

图中的线性迭代优化路径显示，在每一次迭代过程中系统都将朝着最优解的方向前进。值得注意的是该移动轨迹与坐标轴保持平行这是因为每次仅优化单一变量（θ或Q）。

这类似于在一个二维空间中寻找一条曲线的最值点，在此过程中我们发现该函数无法直接计算其导数从而使得基于梯度的信息无法直接应用到此问题中。然而一旦固定其中一个变量则另一个变量可通过求导确定其最优解从而可以采用坐标上升法逐步逼近全局最优解的具体位置这一过程被广泛应用于解决复杂的优化问题例如在统计学习中的期望最大化算法框架下我们可以将这一过程具体化为E步：固定θ参数优化隐含变量Q；M步：固定Q重新估计参数θ；如此往复交替迭代直至收敛达到最优状态

（2）例子

在这个案例中使用抛硬币来说明问题。其中H代表正面朝上的结果而T代表反面朝上的结果。参数θ被定义为正面出现的概率。我们有两枚硬币标记为A和B这两者都是有偏的的。在本次实验中总共设置了5组每组随机选取一枚硬币并连续进行10次抛掷操作。如果能够确定每次抛掷的具体硬币那么计算参数θ的过程就显得相对简单明了如下图所示

如果无法识别每次抛掷使用的硬币类型，则需要采用EM算法。该方法的主要包含以下内容：1. 给定参数θA和θB设定初始值；2. （E-step）计算每组实验更可能使用硬币A的概率（相应地，使用硬币B的概率则为1减去前者）。基于此概率分别计算两者的期望值；3. （M-step）基于上述第三步求得的期望值重新评估参数θA和θB；4. 当达到预设迭代次数或收敛到指定精度时终止运算流程。

稍微解释一下上图的计算过程。初始值θA=0.6,θB=0.5。

图中标注的数值0.45源自何处？基于两个硬币初始设定值分别为0.6和0.5的情况下，在投掷过程中出现五次正面和五次反面的概率为

其中, 0\.{ }{ }{ }{ }{ }{ }{ }{ }{ }{ }{ }{ }{}{}{}{}{}{}{}{}{}{}{}{}{}, 也就是来源于 .}{}{{}}{}{{}}{}{{}}{}{{}}{}{{}}{}{{}}{}{{}}{}{{}}{}{{}}{}{{}}{}{{}}{}{{}}{}{{}}}{}}}}}}}}}} 的一个近似值, 并用来表示在第一组实验中选择硬币 A 的概率为 .}{}{}{}{}{}{}{}{}{}{}{}{} 。图中标注的数据点（如 . H\ 和 . T$\）是如何计算得到的呢？

在第一组实验中采用A硬币并获得正面朝上的次数期望值为2.2。其余各项数值依次类推

解释：

E步：通过初始化参数（初始值θA，θB），分别计算隐藏变量（A或B，未知）的期望值。

M步：基于已知的结果进行分析，在具体实施时例如针对第一组数据而言，则意味着在10次投掷中观察到有2.2H + 2.2枚A型硬币以及2.8H + 2.8H枚B型硬币。继而利用这些结果进行最大似然估计来更新原有的参数（\theta_A, \theta_B），从而实现对初始值的更新。

往复计算，从而收敛到最优参数值。

至此，EM算法总结完毕！

三.高斯混合模型

1.一个例子

Gaussian mixture model represents a class of probabilistic models that combine multiple probability distributions to represent complex data. Theoretically, GMM can approximate a wide range of distributions, making it particularly effective for modeling datasets that exhibit mixed characteristics. This includes scenarios where data conforms to the same underlying distribution type but with differing parameters, or cases involving fundamentally distinct distribution types such as normal and Bernoulli distributions.

从图1可以看出，在我们眼中这些点明显分为两个聚类群。这两个聚类群中的点分别由两个不同的正态分布生成而来。然而，在没有混合高斯模型的情况下（GMM），我们只能用一个二维高斯分布在图1中描述这些数据特征。其等高线代表了两倍标准差范围的区域边界线。这显然显得不合理，因为直观观察就能看出应该将它们分为两类

此时最适合应用GMM的情况是当数据呈现明显的聚类特征时

2.高斯混合模型（GMM）

设有随机变量X ，则混合高斯模型可以用下式表示：

其中N(x | μ_k, Σ_k) 被称为混合模型中的第k个成分（component）。例如，在图2所示的例子中，请注意有两个聚类情况：我们可以用两个二维高斯分布来表示，则分量数目K=2；其中π_k 代表混合系数（mixture coefficient），并且满足：

可以看到πk 相当于每个分量N(x∣μk,Σk) 的权重。

3.GMM的应用

4.GMM参数估计过程

（1）GMM的贝叶斯理解

在介绍GMM参数估计之前，我们先改写GMM的形式，改写之后的GMM模型可以方便地使用EM估计参数。GMM的原始形式如下（1）式：

前面我们已经表示π_k对应于第k个类别的被选中概率。为了方便后续讨论，在此我们定义一个K维随机变量z=(z₁,z₂,…,z_K)，其中每个z_k仅取0或1两个值；当zk=1时，则表示第k个类别被选中（即p(z_k=1)=π_k）；反之，则表示该类别未被选中（即p(z_k=0)=1−π_k）。进一步地，在数学上我们可以定义这些约束条件：

如图2所示的例子中存在两类情况，则变量z的维度为2个维度。若从第一类中选取任一点，则二维变量z=(1,0)，而若从第二类中选取任一点，则变量z=(0,1)= (0,1)。由此可得：zk=1的概率即为πk；假定zk之间相互独立且服从同一分布（independently and identically distributed,iid），则根据公式(2)，我们能够推导出变量Z的联合概率密度函数表达式

由于zk仅为二进制变量，在向量空间中仅有一个分量取值为1而其余分量均为0的状态下所述等式得以满足因此上式成立

按照图2所示的数据分布特征我们可以将这些数据划分为两个类别并且每一类别的数据都呈现出明显的正态分布特征这种分类关系可以通过条件概率模型来进行描述

这一论述可以用条件概率模型来描述：

即第k kk类中的数据服从正态分布。进而上式有可以写成如下形式（3）式：

上述内容分别陈述了变量z及其条件概率p(x∣z)的表达式，在运用条件概率的基本原理后，则能够推导并获得样本边缘概率分布函数f_x(t)，其数学表达如（4）所示：

（注 ：上式第二个等号由于对某个k值而言只要i≠k，则有zi=0 ，因此zk=0 的项为1而无需考虑其影响），最终得到第三个等号）
可以看出GMM模型的第一式与第四式具有相同的结构，并且第四式中引入了一个新的隐含变量z ，通常称为潜在变量（latent variable） 。对于图2所示的数据集，“潜在”的含义是：我们知道数据可以分为两类，并且从外部随机抽取一个数据点时我们无法立即判断该数据点属于哪一类类别成员身份不可见 ，因此通过引入潜在变量z 来描述这一现象。

基于贝叶斯理论的视角来看，在模型构建过程中我们通常假设数据遵循某个特定的概率分布形式。具体而言，在贝叶斯框架下：

• p(z) 代表先验概率分布；
• p(x|z) 被视为似然函数；
  因此，在贝叶斯框架下推导后验分布 p(z|x|y)；
  其中，

\gamma(z_k)=\int p(y|z_k)p(z_k|x,y)\mathrm{d}z_k

定义为第 k 个分量对应的后验分布的概率密度函数；
将其代入贝叶斯公式中，
最终得到方程(5)，即：

\hat{f}_k(x)=\sum_{j=1}^n \gamma_j^{(k)}(x)\phi_j^{(k)}(x)

上述内容对GMM的形式进行了重新表述，并在已知x的情况下引入了隐含变量z及其后验概率γ(zk)，这样做的目的是为了便于应用EM算法来估计GMM的参数。

（2）EM算法估计GMM参数

该算法分为两个阶段，在第一个阶段中预估需要估计的参数的大致数值，在第二个阶段中则通过最大化似然函数来优化这些参数。因此必须首先计算混合高斯模型（GMM）的似然函数才能进行后续优化步骤。

设输入样本集为X={x₁,x₂,…,x_N}。在图2中,X中的每个样本点具有二维坐标（以粗体形式表示），因此变量如x₁,x₂等均采用此方式标记。其概率密度函数由下式给出：
p(x|π,μ,Σ)=∑_{k=1}^K π_k N(x|μ_k,Σ_k)
其中参数估计问题涉及三个关键量：混合比例π、均值向量μ以及协方差矩阵Σ. 通过重新排列公式项可得（6）式：

为了估算这三个参数的值 需要依次计算出每个参数的最大似然函数 然后计算\mu_k的最大似然函数 接着将(6)式取自然对数后并对\mu_k进行微分运算并将导数值设为零从而得出最大似然估计值

（部分求导过程可以参考正态分布的求导 ）

参考文献

[1] 对似然函数的认识：https://www.cnblogs.com/zhsuiy/p/4822020.html
[2] 先验概率、后验概率及似然概率的概念分别来自不同的统计学理论。<>
[3] 统计学习笔记（九）——对EM算法的深入分析：<>
[4] 联合分布的概率模型在机器学习中具有重要应用。https://baike.baidu.com/item/联合概率分布/12611723?fr=aladdin
[5] 高斯混合模型（GMM）及其EM算法的掌握：<>