Expectation-maximization： 高斯混合模型的期望最大化算法

高斯混合模型的期望最大化算法

本文大部分内容翻译自Answers.com的Expectation-maximization algorithm 词条。

对于一个不完整的样本集，假设y为观测（已知）数据，z为丢失数据，z和y构成了完整的数据。z可以为丢失的观测量，也可以是隐含的变量，如果这些变量已知的话将会简化问题的求解。例如，在混合模型中，如果产生样本点的各个成分的分布为已知，则可以大大简化似然函数的计算。

假设p为完整数据的联合概率密度函数，该函数具有参数θ：，它也可以看成是完整数据的似然函数（即样本集(y,z)具有θ分布函数的相似性），一个关于θ的函数。进一步，根据贝叶斯公式和全概率公式，丢失数据对于观测数据的条件分布为：

这个式子仅要求计算观测数据对于丢失数据的条件概率（似然函数）和丢失数据的分布概率。（分母仅为归一化需要）

EM算法通过迭代方法，构造更好的 θ₁,θ₂,... ，逐步优化初始的 θ₀， θ的递推优化公式通过下式得到：

Q(θ)是对数似然函数ln的期望。换句话说，我们不知道所有数据的值，所以我们不知道似然函数的准确值，但是对于给定的观测数据y，我们可以得到未知的z的概率分布，对于每一个丢失的观测（z中的一个元素），都有一个关于θ的似然值，这样我们就可以计算似然函数的数学期望。这个期望显然受到θ的值的影响。

Q(θ)通过下式得到：

或者，更普遍地，可以写成条件期望的形式

应该清楚， 这里关于的条件期望中用到的z的条件分布（前一个式子中的第一项）是在固定θ_n的情况下取得的。用对数似然代替直接的似然函数，可以在保证极值位置不变的情况下简化计算。也就是说，θ_{n + 1}是通过在θ_n 的基础上，最大化所有数据关于已知数据的条件期望E而取得的。

当样本集合的最大似然估计比较容易获得时，EM算法更为有用。如果似然估计具有给定的形式，则M步就非常简单。一个经典的例子是用于有限的高斯混合模型的最大似然估计，如果混合分布已知，则可以很容易得到各个成分的分布。

高斯混合模型

假设n个高斯分布，从中产生m个样本,每一个 y_j 所从属的分布用 z_j 表示，z_j 的取值范围为1到n。对于任意的y，其来自第i个高斯分布的概率为

我们的任务是估计未知的参数，即每一个高斯分布的均值、方差和先验概率。

首先我们必须设定初始值。

对于无监督的聚类，我们首先必须为每一个样本指定一个类别，这一步可以通过其他的聚类方法实现，如K-means方法，求出各个类别的中心和每一个样本的类别。然后求出各个类别中样本的协方差阵，可以用每个类别中样本的个数来表示该类别的权重（先验概率）。

E-step

使用上一次M步中得到的参数θ，计算未知的z（即每一个分类）对于观测数据的条件分布。即，在当前参数θ_n 下，对每一个样本，计算各个类别的后验概率。 下式中第三个等号后面的分子的第一项，为每一个样本在类别i中的分布概率，分子中的第二项为每一个类别的先验概率。
分母为归一化，即对每一个样本，其在各个类别中的分布概率的和为1.

for ti = 1:Ncluster
PyatCi = mvnpdf(featurespace,Mu(ti,:),Sigma(:,:,ti));
E(:,ti) = PyatCi*Weight(ti);
end
Esum = sum(E,2);
for ti = 1:Ncluster
E(:,ti) = E(:,ti)./Esum;
end

M-step

为了最大化联合分布的对数似然函数的期望

上式中最后一项的联合分布可以写成条件分布和边缘分布的乘积的形式，得到

上式需要满足先验概率归一化的要求：

为了求得Q(θ)的极大值的位置，引入拉格朗日算子 ，然后带入概率密度函数，

为了求得新的 θ_{t + 1},我们需要找出上式取得极值时对应的θ，即其一阶导数。

下面推导满足上式的均值，

for ti = 1:Ncluster
muti = E(:,ti)'*featurespace;
muti = muti/sum(E(:,ti));
Mu(ti,:) = muti;
end

产生新的协方差矩阵

新的各个分类的权重

进行归一化处理

Inserting λ into our estimate:

将新的参数，均值、方差、分类的先验概率，带入E步，重新计算，直到样本集合对各个分类的似然函数不再有明显的变化为止。

下面的程序提供了一种快速的似然函数计算方法，摘自Patrick Tsui的EM_GM样例程序：
http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=8636

function L = Likelihood(X,k,W,M,V)
% Compute L based on K. V. Mardia, "Multivariate Analysis", Academic Press, 1979, PP. 96-97
% to enchance computational speed
% a subfunction from Patrick P. C. Tsui's EM_GM code.
[n,d] = size(X);
U = mean(X)';
S = cov(X);
L = 0;
for i=1:k,
iV = inv(V(:,:,i));
L = L + W(i)*(-0.5*n*log(det(2*pi*V(:,:,i))) ...
-0.5*(n-1)*(trace(iV*S)+(U-M(:,i))'*iV*(U-M(:,i))));
end