（EM算法）The EM Algorithm 和 混合高斯模型Demo(转载)
2017-03-04 00:43
459 查看
转载来源：http://www.cnblogs.com/jerrylead/archive/2011/04/06/2006936.html
这篇讨论使用期望最大化算法（Expectation-Maximization）来进行密度估计（density estimation）。
与k-means一样，给定的训练样本是
，我们将隐含类别标签用
表示。与k-means的硬指定不同，我们首先认为
是满足一定的概率分布的，这里我们认为满足多项式分布，
，其中
，
有k个值{1,…,k}可以选取。而且我们认为在给定
后，
满足多值高斯分布，即
。由此可以得到联合分布
。
整个模型简单描述为对于每个样例
，我们先从k个类别中按多项式分布抽取一个
，然后根据
所对应的k个多值高斯分布中的一个生成样例
，。整个过程称作混合高斯模型。注意的是这里的
仍然是隐含随机变量。模型中还有三个变量
和
。最大似然估计为
。对数化后如下：
这个式子的最大值是不能通过前面使用的求导数为0的方法解决的，因为求的结果不是close form。但是假设我们知道了每个样例的
，那么上式可以简化为：
这时候我们再来对
和
进行求导得到：
就是样本类别中
的比率。
是类别为j的样本特征均值，
是类别为j的样例的特征的协方差矩阵。
实际上，当知道
后，最大似然估计就近似于高斯判别分析模型（Gaussian
discriminant analysis model）了。所不同的是GDA中类别y是伯努利分布，而这里的z是多项式分布，还有这里的每个样例都有不同的协方差矩阵，而GDA中认为只有一个。
之前我们是假设给定了
，实际上
是不知道的。那么怎么办呢？考虑之前提到的EM的思想，第一步是猜测隐含类别变量z，第二步是更新其他参数，以获得最大的最大似然估计。用到这里就是：
在E步中，我们将其他参数
看作常量，计算
的后验概率，也就是估计隐含类别变量。估计好后，利用上面的公式重新计算其他参数，计算好后发现最大化最大似然估计时，
值又不对了，需要重新计算，周而复始，直至收敛。
的具体计算公式如下：
这个式子利用了贝叶斯公式。
这里我们使用
代替了前面的
，由简单的0/1值变成了概率值。
对比K-means可以发现，这里使用了“软”指定，为每个样例分配的类别
是有一定的概率的，同时计算量也变大了，每个样例i都要计算属于每一个类别j的概率。与K-means相同的是，结果仍然是局部最优解。对其他参数取不同的初始值进行多次计算不失为一种好方法。
下面主要介绍EM的整个推导过程。
1. Jensen不等式
回顾优化理论中的一些概念。设f是定义域为实数的函数，如果对于所有的实数x，
，那么f是凸函数。当x是向量时，如果其hessian矩阵H是半正定的（
），那么f是凸函数。如果
或者
，那么称f是严格凸函数。
Jensen不等式表述如下：
如果f是凸函数，X是随机变量，那么
特别地，如果f是严格凸函数，那么
当且仅当
，也就是说X是常量。
这里我们将
简写为
。
如果用图表示会很清晰：
图中，实线f是凸函数，X是随机变量，有0.5的概率是a，有0.5的概率是b。（就像掷硬币一样）。X的期望值就是a和b的中值了，图中可以看到
成立。
当f是（严格）凹函数当且仅当-f是（严格）凸函数。
Jensen不等式应用于凹函数时，不等号方向反向，也就是
。
2. EM算法
给定的训练样本是
，样例间独立，我们想找到每个样例隐含的类别z，能使得p(x,z)最大。p(x,z)的最大似然估计如下：
第一步是对极大似然取对数，第二步是对每个样例的每个可能类别z求联合分布概率和。但是直接求
一般比较困难，因为有隐藏变量z存在，但是一般确定了z后，求解就容易了。
EM是一种解决存在隐含变量优化问题的有效方法。竟然不能直接最大化
，我们可以不断地建立
的下界（E步），然后优化下界（M步）。这句话比较抽象，看下面的。
对于每一个样例i，让
表示该样例隐含变量z的某种分布，
满足的条件是
。（如果z是连续性的，那么
是概率密度函数，需要将求和符号换做积分符号）。比如要将班上学生聚类，假设隐藏变量z是身高，那么就是连续的高斯分布。如果按照隐藏变量是男女，那么就是伯努利分布了。
可以由前面阐述的内容得到下面的公式：
（1）到（2）比较直接，就是分子分母同乘以一个相等的函数。（2）到（3）利用了Jensen不等式，考虑到
是凹函数（二阶导数小于0），而且
就是
的期望（回想期望公式中的Lazy
Statistician规则）
对应于上述问题，Y是
，X是
，
是
，g是
到
的映射。这样解释了式子（2）中的期望，再根据凹函数时的Jensen不等式：
可以得到（3）。
这个过程可以看作是对
求了下界。对于
的选择，有多种可能，那种更好的？假设
已经给定，那么
的值就决定于
和
了。我们可以通过调整这两个概率使下界不断上升，以逼近
的真实值，那么什么时候算是调整好了呢？当不等式变成等式时，说明我们调整后的概率能够等价于
了。按照这个思路，我们要找到等式成立的条件。根据Jensen不等式，要想让等式成立，需要让随机变量变成常数值，这里得到：
c为常数，不依赖于
。对此式子做进一步推导，我们知道
，那么也就有
，（多个等式分子分母相加不变，这个认为每个样例的两个概率比值都是c），那么有下式：
至此，我们推出了在固定其他参数
后，
的计算公式就是后验概率，解决了
如何选择的问题。这一步就是E步，建立
的下界。接下来的M步，就是在给定
后，调整
，去极大化
的下界（在固定
后，下界还可以调整的更大）。那么一般的EM算法的步骤如下：
那么究竟怎么确保EM收敛？假定
和
是EM第t次和t+1次迭代后的结果。如果我们证明了
，也就是说极大似然估计单调增加，那么最终我们会到达最大似然估计的最大值。下面来证明，选定
后，我们得到E步
这一步保证了在给定
时，Jensen不等式中的等式成立，也就是
然后进行M步，固定
，并将
视作变量，对上面的
求导后，得到
，这样经过一些推导会有以下式子成立：
解释第（4）步，得到
时，只是最大化
，也就是
的下界，而没有使等式成立，等式成立只有是在固定
，并按E步得到
时才能成立。
况且根据我们前面得到的下式，对于所有的
和
都成立
第（5）步利用了M步的定义，M步就是将
调整到
，使得下界最大化。因此（5）成立，（6）是之前的等式结果。
这样就证明了
会单调增加。一种收敛方法是
不再变化，还有一种就是变化幅度很小。
再次解释一下（4）、（5）、（6）。首先（4）对所有的参数都满足，而其等式成立条件只是在固定
，并调整好Q时成立，而第（4）步只是固定Q，调整
，不能保证等式一定成立。（4）到（5）就是M步的定义，（5）到（6）是前面E步所保证等式成立条件。也就是说E步会将下界拉到与
一个特定值（这里
）一样的高度，而此时发现下界仍然可以上升，因此经过M步后，下界又被拉升，但达不到与
另外一个特定值一样的高度，之后E步又将下界拉到与这个特定值一样的高度，重复下去，直到最大值。
如果我们定义
从前面的推导中我们知道
，EM可以看作是J的坐标上升法，E步固定
，优化
，M步固定
优化
。
3. 重新审视混合高斯模型
我们已经知道了EM的精髓和推导过程，再次审视一下混合高斯模型。之前提到的混合高斯模型的参数
和
计算公式都是根据很多假定得出的，有些没有说明来由。为了简单，这里在M步只给出
和
的推导方法。
E步很简单，按照一般EM公式得到：
简单解释就是每个样例i的隐含类别
为j的概率可以通过后验概率计算得到。
在M步中，我们需要在固定
后最大化最大似然估计，也就是
这是将
的k种情况展开后的样子，未知参数
和
。
固定
和
，对
求导得
等于0时，得到
这就是我们之前模型中的
的更新公式。
然后推导
的更新公式。看之前得到的
在
和
确定后，分子上面的一串都是常数了，实际上需要优化的公式是：
需要知道的是，
还需要满足一定的约束条件就是
。
这个优化问题我们很熟悉了，直接构造拉格朗日乘子。
还有一点就是
，但这一点会在得到的公式里自动满足。
求导得，
等于0，得到
也就是说
再次使用
，得到
这样就神奇地得到了
。
那么就顺势得到M步中
的更新公式：
的推导也类似，不过稍微复杂一些，毕竟是矩阵。结果在之前的混合高斯模型中已经给出。
4. 总结
如果将样本看作观察值，潜在类别看作是隐藏变量，那么聚类问题也就是参数估计问题，只不过聚类问题中参数分为隐含类别变量和其他参数，这犹如在x-y坐标系中找一个曲线的极值，然而曲线函数不能直接求导，因此什么梯度下降方法就不适用了。但固定一个变量后，另外一个可以通过求导得到，因此可以使用坐标上升法，一次固定一个变量，对另外的求极值，最后逐步逼近极值。对应到EM上，E步估计隐含变量，M步估计其他参数，交替将极值推向最大。EM中还有“硬”指定和“软”指定的概念，“软”指定看似更为合理，但计算量要大，“硬”指定在某些场合如K-means中更为实用（要是保持一个样本点到其他所有中心的概率，就会很麻烦）。
另外，EM的收敛性证明方法确实很牛，能够利用log的凹函数性质，还能够想到利用创造下界，拉平函数下界，优化下界的方法来逐步逼近极大值。而且每一步迭代都能保证是单调的。最重要的是证明的数学公式非常精妙，硬是分子分母都乘以z的概率变成期望来套上Jensen不等式，前人都是怎么想到的。
在Mitchell的Machine Learning书中也举了一个EM应用的例子，明白地说就是将班上学生的身高都放在一起，要求聚成两个类。这些身高可以看作是男生身高的高斯分布和女生身高的高斯分布组成。因此变成了如何估计每个样例是男生还是女生，然后在确定男女生情况下，如何估计均值和方差，里面也给出了公式，有兴趣可以参考。
% HSW
% 2017-3-3
% 混合高斯模型 Mixed Gaussian Model
%
clc;
clear all;
close all;
% 生成学习样本
feature_num = 3;
mu = [2 3 4];
SIGMA = [1,0,0;0,1,0;0,0,1];
training_x_positive = mvnrnd(mu,SIGMA,1000); % 正例
mu = [4 5 6];
SIGMA = [1,0,0;0,1,0;0,0,1];
training_x_negative = mvnrnd(mu,SIGMA,1000); % 反例
mu = [2 3 4];
SIGMA = [1,0,0;0,1,0;0,0,1];
testing_x_positive = mvnrnd(mu,SIGMA,1000); % 正例
mu = [5 6 7];
SIGMA = [1,0,0;0,1,0;0,0,1];
testing_x_negative = mvnrnd(mu,SIGMA,1000); % 反例
training_y_positive = ones(1,size(training_x_positive,1)); % 正例标签
training_y_negative = 2*ones(1,size(training_x_negative,1));% 反例标签
testing_y_positive = ones(1,size(training_x_positive,1)); % 正例标签
testing_y_negative = 2*ones(1,size(training_x_negative,1));% 反例标签
figure;
plot(training_x_positive(:,1),training_x_positive(:,2),'r+');
hold on;
plot(training_x_negative(:,1),training_x_negative(:,2),'b*');
title('正例(+)和反例(*)');
%
training_x = cat(2,training_x_positive',training_x_negative');
testing_x = cat(2,testing_x_positive',testing_x_negative');
training_y = cat(2,training_y_positive,training_y_negative);
testing_y = cat(2,testing_y_positive,testing_y_negative);
% 样本个数
Sample_num = size(training_x,2);
% 打乱样本
re_rank = randperm(Sample_num);
training_x = training_x(:,re_rank);
training_y = training_y(:,re_rank);
% 混合高斯模型的个数
Gaussian_num = 2;
% % 初始化，采用随机指定每个样本的类别，但是结果不是很好
% for iter_Gaussian = 1:Gaussian_num
% mu(:,iter_Gaussian) = rand(feature_num,1);
% fi(1,iter_Gaussian) = 1/Gaussian_num;
% sigma{1,iter_Gaussian} = eye(feature_num);
% end% for iter_Gaussian
% 初始化采用k-mean聚类
Cluster_sample = training_x;
Cluster_sample_num = Sample_num;
ClusterNum = Gaussian_num;
Max_iter_Num = 1000;
Cluster_index = zeros(Cluster_sample_num,1);% 样本的类别
mu = ones(feature_num, ClusterNum);% 聚类中心
Loss = zeros(Max_iter_Num,1);
Distortion_Loss_Old = 1e10;
Distortion_Loss = 0;
Count = 0;
epsion = 1e-6;
for iter = 1:Max_iter_Num
Count = Count + 1;
if iter == 1
% 随机初始化聚类中心
[cluster_index] = randperm(Cluster_sample_num,2);
mu(:,1) = Cluster_sample(:,cluster_index(1));
mu(:,2) = Cluster_sample(:,cluster_index(2));
cluster_mu_1 = sum((Cluster_sample - repmat(mu(:,1),[1,Cluster_sample_num])).^2,1);
cluster_mu_2 = sum((Cluster_sample - repmat(mu(:,2),[1,Cluster_sample_num])).^2,1);
cluster_one_index = find( (cluster_mu_1 - cluster_mu_2) < 0);
cluster_two_index = find( (cluster_mu_1 - cluster_mu_2) >= 0);
Cluster_index(cluster_one_index) = 1;
Cluster_index(cluster_two_index) = 2;
Distortion_Loss = sum(sum((Cluster_sample(:,cluster_one_index) - repmat(mu(:,1),[1,length(cluster_one_index)])).^2,1)) + sum(sum((Cluster_sample(:,cluster_two_index) - repmat(mu(:,2),[1,length(cluster_two_index)])).^2,1)); % 初始的
mu(:,1) = mean(Cluster_sample(:,cluster_one_index),2); % 更新聚类中心
mu(:,2) = mean(Cluster_sample(:,cluster_two_index),2);
else
cluster_mu_1 = sum((Cluster_sample - repmat(mu(:,1),[1,Cluster_sample_num])).^2,1);
cluster_mu_2 = sum((Cluster_sample - repmat(mu(:,2),[1,Cluster_sample_num])).^2,1);
cluster_one_index = find( (cluster_mu_1 - cluster_mu_2) < 0);
cluster_two_index = find( (cluster_mu_1 - cluster_mu_2) >= 0);
Cluster_index(cluster_one_index) = 1;
Cluster_index(cluster_two_index) = 2;
Distortion_Loss = sum(sum((Cluster_sample(:,cluster_one_index) - repmat(mu(:,1),[1,length(cluster_one_index)])).^2,1)) + sum(sum((Cluster_sample(:,cluster_two_index) - repmat(mu(:,2),[1,length(cluster_two_index)])).^2,1));
mu(:,1) = mean(Cluster_sample(:,cluster_one_index),2); % 更新聚类中心
mu(:,2) = mean(Cluster_sample(:,cluster_two_index),2);
end
Loss(Count) = Distortion_Loss;
if abs(Distortion_Loss_Old - Distortion_Loss) < epsion
break;
end
Distortion_Loss_Old = Distortion_Loss;
end
%显示聚类结果
figure;
for iter = 1:Cluster_sample_num
hold on;
if Cluster_index(iter) == 1
plot3(Cluster_sample(1,iter),Cluster_sample(2,iter),Cluster_sample(3,iter),'g+');
else
plot3(Cluster_sample(1,iter),Cluster_sample(2,iter),Cluster_sample(3,iter),'b*');
end
end
% 显示聚类中心
hold on
plot3(mu(1,1),mu(2,1),mu(3,1),'rx');
hold on
plot3(mu(1,2),mu(2,2),mu(3,2),'rx');
figure;
plot(1:Count,Loss(1:Count,1),'r*');
title('畸变函数收敛曲线');
sigma = cell(1,Gaussian_num);
mu_old = 100000*mu;
epsion = 1e-8;
fi = zeros(1,Gaussian_num);
for iter_Gaussian = 1:Gaussian_num
Cluster_index_k = find(Cluster_index == iter_Gaussian);
fi(1,iter_Gaussian) = length(Cluster_index_k)/Sample_num;
sigma{1,iter_Gaussian} = zeros(feature_num);
for iter_Sample = 1:Sample_num
if Cluster_index(iter_Sample) == iter_Gaussian
sigma{1,iter_Gaussian} = sigma{1,iter_Gaussian} + (training_x(:,iter_Sample) - mu(:,iter_Gaussian))*(training_x(:,iter_Sample) - mu(:,iter_Gaussian))';
end
end
sigma{1,iter_Gaussian} = sigma{1,iter_Gaussian}/length(Cluster_index_k);
end
disp('初始高斯混合模型1参数');
fi(1,1)
mu(:,1)
sigma{1,1}
disp('初始高斯混合模型2参数');
fi(1,2)
mu(:,2)
sigma{1,2}
W = zeros(Gaussian_num,Sample_num);
% 进行迭代更新
iter_Gaussian_Max_num = 1000;
for iter = 1:iter_Gaussian_Max_num
% E-M算法学习混合高斯模型
% E-Step: 计算W_i_j = p(x^i | z^i = j; mu, sigma)*p(z^i = j,fi)/sum_l={1,,,k}(p(x^i | z^i = l; mu, sigma)*p(z^i = l, fi))
for iter_Sample = 1:Sample_num
for iter_Gaussian = 1:Gaussian_num
W(iter_Gaussian,iter_Sample) = 1/((2*pi)^(feature_num/2)*sqrt(det(sigma{1,iter_Gaussian})))*exp(-1/2*(training_x(:,iter_Sample)- mu(:,iter_Gaussian))'/sigma{1,iter_Gaussian}*(training_x(:,iter_Sample) - mu(:,iter_Gaussian)))*fi(1,iter_Gaussian);
end
end
W = W./repmat(sum(W,1),[Gaussian_num,1]);
% M-Step:
for iter_Gaussian =1:Gaussian_num
fi(1,iter_Gaussian) = sum(W(iter_Gaussian,:))/Sample_num;
mu(:,iter_Gaussian) = sum(repmat(W(iter_Gaussian,:),[feature_num,1]).*training_x,2);
mu(:,iter_Gaussian) = mu(:,iter_Gaussian)./sum(W(iter_Gaussian,:));
sigma{1,iter_Gaussian} = zeros(feature_num);
for iter_Sample = 1:Sample_num
sigma{1,iter_Gaussian} = sigma{1,iter_Gaussian} + W(iter_Gaussian,iter_Sample)*(training_x(:,iter_Sample) - mu(:,iter_Gaussian))*(training_x(:,iter_Sample) - mu(:,iter_Gaussian))';
end% for iter_Sample
sigma{1,iter_Gaussian} = sigma{1,iter_Gaussian}./sum(W(iter_Gaussian,:));
end% for iter_Gaussian
if mod(iter,100) == 0
disp(['迭代第',num2str(iter),'次']);
disp('高斯混合模型1参数');
fi(1,1)
mu(:,1)
sigma{1,1}
disp('高斯混合模型2参数');
fi(1,2)
mu(:,2)
sigma{1,2}
end
if sum(sum(abs(mu_old - mu))) < epsion
break;
end
mu_old = mu;
end
disp('学习获得的高斯混合模型1参数');
fi(1,1)
mu(:,1)
sigma{1,1}
disp('学习获得的高斯混合模型2参数');
fi(1,2)
mu(:,2)
sigma{1,2}
结果如下：
初始高斯混合模型1参数
ans =
0.5125
ans =
3.9788
4.9553
5.9916
ans =
0.9897 -0.0364 -0.0609
-0.0364 0.9018 -0.0459
-0.0609 -0.0459 0.9617
初始高斯混合模型2参数
ans =
0.4875
ans =
1.9684
3.0232
3.9805
ans =
1.0290 0.0047 -0.0530
0.0047 0.8828 -0.0519
-0.0530 -0.0519 0.9517
迭代第100次
高斯混合模型1参数
ans =
0.4930
ans =
3.9830
4.9807
6.0065
ans =
1.0221 -0.0224 -0.0337
-0.0224 0.8995 -0.0430
-0.0337 -0.0430 0.9804
高斯混合模型2参数
ans =
0.5070
ans =
2.0418
3.0729
4.0436
ans =
1.1297 0.0778 0.0340
0.0778 0.9309 0.0141
0.0340 0.0141 1.0275
学习获得的高斯混合模型1参数
ans =
0.4930
ans =
3.9830
4.9807
6.0065
ans =
1.0221 -0.0224 -0.0337
-0.0224 0.8995 -0.0430
-0.0337 -0.0430 0.9804
学习获得的高斯混合模型2参数
ans =
0.5070
ans =
2.0418
3.0729
4.0436
ans =
1.1297 0.0778 0.0340
0.0778 0.9309 0.0141
0.0340 0.0141 1.0275
这篇讨论使用期望最大化算法（Expectation-Maximization）来进行密度估计（density estimation）。
与k-means一样，给定的训练样本是
，我们将隐含类别标签用
表示。与k-means的硬指定不同，我们首先认为
是满足一定的概率分布的，这里我们认为满足多项式分布，
，其中
，
有k个值{1,…,k}可以选取。而且我们认为在给定
后，
满足多值高斯分布，即
。由此可以得到联合分布
。
整个模型简单描述为对于每个样例
，我们先从k个类别中按多项式分布抽取一个
，然后根据
所对应的k个多值高斯分布中的一个生成样例
，。整个过程称作混合高斯模型。注意的是这里的
仍然是隐含随机变量。模型中还有三个变量
和
。最大似然估计为
。对数化后如下：
这个式子的最大值是不能通过前面使用的求导数为0的方法解决的，因为求的结果不是close form。但是假设我们知道了每个样例的
，那么上式可以简化为：
这时候我们再来对
和
进行求导得到：
就是样本类别中
的比率。
是类别为j的样本特征均值，
是类别为j的样例的特征的协方差矩阵。
实际上，当知道
后，最大似然估计就近似于高斯判别分析模型（Gaussian
discriminant analysis model）了。所不同的是GDA中类别y是伯努利分布，而这里的z是多项式分布，还有这里的每个样例都有不同的协方差矩阵，而GDA中认为只有一个。
之前我们是假设给定了
，实际上
是不知道的。那么怎么办呢？考虑之前提到的EM的思想，第一步是猜测隐含类别变量z，第二步是更新其他参数，以获得最大的最大似然估计。用到这里就是：
| 循环下面步骤，直到收敛： { （E步）对于每一个i和j，计算  （M步），更新参数：  } |
看作常量，计算
的后验概率，也就是估计隐含类别变量。估计好后，利用上面的公式重新计算其他参数，计算好后发现最大化最大似然估计时，
值又不对了，需要重新计算，周而复始，直至收敛。
的具体计算公式如下：
这个式子利用了贝叶斯公式。
这里我们使用
代替了前面的
，由简单的0/1值变成了概率值。
对比K-means可以发现，这里使用了“软”指定，为每个样例分配的类别
是有一定的概率的，同时计算量也变大了，每个样例i都要计算属于每一个类别j的概率。与K-means相同的是，结果仍然是局部最优解。对其他参数取不同的初始值进行多次计算不失为一种好方法。
下面主要介绍EM的整个推导过程。
1. Jensen不等式
回顾优化理论中的一些概念。设f是定义域为实数的函数，如果对于所有的实数x，
，那么f是凸函数。当x是向量时，如果其hessian矩阵H是半正定的（
），那么f是凸函数。如果
或者
，那么称f是严格凸函数。
Jensen不等式表述如下：
如果f是凸函数，X是随机变量，那么
特别地，如果f是严格凸函数，那么
当且仅当
，也就是说X是常量。
这里我们将
简写为
。
如果用图表示会很清晰：
图中，实线f是凸函数，X是随机变量，有0.5的概率是a，有0.5的概率是b。（就像掷硬币一样）。X的期望值就是a和b的中值了，图中可以看到
成立。
当f是（严格）凹函数当且仅当-f是（严格）凸函数。
Jensen不等式应用于凹函数时，不等号方向反向，也就是
。
2. EM算法
给定的训练样本是
，样例间独立，我们想找到每个样例隐含的类别z，能使得p(x,z)最大。p(x,z)的最大似然估计如下：
第一步是对极大似然取对数，第二步是对每个样例的每个可能类别z求联合分布概率和。但是直接求
一般比较困难，因为有隐藏变量z存在，但是一般确定了z后，求解就容易了。
EM是一种解决存在隐含变量优化问题的有效方法。竟然不能直接最大化
，我们可以不断地建立
的下界（E步），然后优化下界（M步）。这句话比较抽象，看下面的。
对于每一个样例i，让
表示该样例隐含变量z的某种分布，
满足的条件是
。（如果z是连续性的，那么
是概率密度函数，需要将求和符号换做积分符号）。比如要将班上学生聚类，假设隐藏变量z是身高，那么就是连续的高斯分布。如果按照隐藏变量是男女，那么就是伯努利分布了。
可以由前面阐述的内容得到下面的公式：
（1）到（2）比较直接，就是分子分母同乘以一个相等的函数。（2）到（3）利用了Jensen不等式，考虑到
是凹函数（二阶导数小于0），而且
就是
的期望（回想期望公式中的Lazy
Statistician规则）
设Y是随机变量X的函数 （g是连续函数），那么 （1） X是离散型随机变量，它的分布律为  ，k=1,2,…。若  绝对收敛，则有  （2） X是连续型随机变量，它的概率密度为  ，若  绝对收敛，则有  |
，X是
，
是
，g是
到
的映射。这样解释了式子（2）中的期望，再根据凹函数时的Jensen不等式：
可以得到（3）。
这个过程可以看作是对
求了下界。对于
的选择，有多种可能，那种更好的？假设
已经给定，那么
的值就决定于
和
了。我们可以通过调整这两个概率使下界不断上升，以逼近
的真实值，那么什么时候算是调整好了呢？当不等式变成等式时，说明我们调整后的概率能够等价于
了。按照这个思路，我们要找到等式成立的条件。根据Jensen不等式，要想让等式成立，需要让随机变量变成常数值，这里得到：
c为常数，不依赖于
。对此式子做进一步推导，我们知道
，那么也就有
，（多个等式分子分母相加不变，这个认为每个样例的两个概率比值都是c），那么有下式：
至此，我们推出了在固定其他参数
后，
的计算公式就是后验概率，解决了
如何选择的问题。这一步就是E步，建立
的下界。接下来的M步，就是在给定
后，调整
，去极大化
的下界（在固定
后，下界还可以调整的更大）。那么一般的EM算法的步骤如下：
| 循环重复直到收敛 { （E步）对于每一个i，计算  （M步）计算  |
和
是EM第t次和t+1次迭代后的结果。如果我们证明了
，也就是说极大似然估计单调增加，那么最终我们会到达最大似然估计的最大值。下面来证明，选定
后，我们得到E步
这一步保证了在给定
时，Jensen不等式中的等式成立，也就是
然后进行M步，固定
，并将
视作变量，对上面的
求导后，得到
，这样经过一些推导会有以下式子成立：
解释第（4）步，得到
时，只是最大化
，也就是
的下界，而没有使等式成立，等式成立只有是在固定
，并按E步得到
时才能成立。
况且根据我们前面得到的下式，对于所有的
和
都成立
第（5）步利用了M步的定义，M步就是将
调整到
，使得下界最大化。因此（5）成立，（6）是之前的等式结果。
这样就证明了
会单调增加。一种收敛方法是
不再变化，还有一种就是变化幅度很小。
再次解释一下（4）、（5）、（6）。首先（4）对所有的参数都满足，而其等式成立条件只是在固定
，并调整好Q时成立，而第（4）步只是固定Q，调整
，不能保证等式一定成立。（4）到（5）就是M步的定义，（5）到（6）是前面E步所保证等式成立条件。也就是说E步会将下界拉到与
一个特定值（这里
）一样的高度，而此时发现下界仍然可以上升，因此经过M步后，下界又被拉升，但达不到与
另外一个特定值一样的高度，之后E步又将下界拉到与这个特定值一样的高度，重复下去，直到最大值。
如果我们定义
从前面的推导中我们知道
，EM可以看作是J的坐标上升法，E步固定
，优化
，M步固定
优化
。
3. 重新审视混合高斯模型
我们已经知道了EM的精髓和推导过程，再次审视一下混合高斯模型。之前提到的混合高斯模型的参数
和
计算公式都是根据很多假定得出的，有些没有说明来由。为了简单，这里在M步只给出
和
的推导方法。
E步很简单，按照一般EM公式得到：
简单解释就是每个样例i的隐含类别
为j的概率可以通过后验概率计算得到。
在M步中，我们需要在固定
后最大化最大似然估计，也就是
这是将
的k种情况展开后的样子，未知参数
和
。
固定
和
，对
求导得
等于0时，得到
这就是我们之前模型中的
的更新公式。
然后推导
的更新公式。看之前得到的
在
和
确定后，分子上面的一串都是常数了，实际上需要优化的公式是：
需要知道的是，
还需要满足一定的约束条件就是
。
这个优化问题我们很熟悉了，直接构造拉格朗日乘子。
还有一点就是
，但这一点会在得到的公式里自动满足。
求导得，
等于0，得到
也就是说
再次使用
，得到
这样就神奇地得到了
。
那么就顺势得到M步中
的更新公式：
的推导也类似，不过稍微复杂一些，毕竟是矩阵。结果在之前的混合高斯模型中已经给出。
4. 总结
如果将样本看作观察值，潜在类别看作是隐藏变量，那么聚类问题也就是参数估计问题，只不过聚类问题中参数分为隐含类别变量和其他参数，这犹如在x-y坐标系中找一个曲线的极值，然而曲线函数不能直接求导，因此什么梯度下降方法就不适用了。但固定一个变量后，另外一个可以通过求导得到，因此可以使用坐标上升法，一次固定一个变量，对另外的求极值，最后逐步逼近极值。对应到EM上，E步估计隐含变量，M步估计其他参数，交替将极值推向最大。EM中还有“硬”指定和“软”指定的概念，“软”指定看似更为合理，但计算量要大，“硬”指定在某些场合如K-means中更为实用（要是保持一个样本点到其他所有中心的概率，就会很麻烦）。
另外，EM的收敛性证明方法确实很牛，能够利用log的凹函数性质，还能够想到利用创造下界，拉平函数下界，优化下界的方法来逐步逼近极大值。而且每一步迭代都能保证是单调的。最重要的是证明的数学公式非常精妙，硬是分子分母都乘以z的概率变成期望来套上Jensen不等式，前人都是怎么想到的。
在Mitchell的Machine Learning书中也举了一个EM应用的例子，明白地说就是将班上学生的身高都放在一起，要求聚成两个类。这些身高可以看作是男生身高的高斯分布和女生身高的高斯分布组成。因此变成了如何估计每个样例是男生还是女生，然后在确定男女生情况下，如何估计均值和方差，里面也给出了公式，有兴趣可以参考。
% HSW
% 2017-3-3
% 混合高斯模型 Mixed Gaussian Model
%
clc;
clear all;
close all;
% 生成学习样本
feature_num = 3;
mu = [2 3 4];
SIGMA = [1,0,0;0,1,0;0,0,1];
training_x_positive = mvnrnd(mu,SIGMA,1000); % 正例
mu = [4 5 6];
SIGMA = [1,0,0;0,1,0;0,0,1];
training_x_negative = mvnrnd(mu,SIGMA,1000); % 反例
mu = [2 3 4];
SIGMA = [1,0,0;0,1,0;0,0,1];
testing_x_positive = mvnrnd(mu,SIGMA,1000); % 正例
mu = [5 6 7];
SIGMA = [1,0,0;0,1,0;0,0,1];
testing_x_negative = mvnrnd(mu,SIGMA,1000); % 反例
training_y_positive = ones(1,size(training_x_positive,1)); % 正例标签
training_y_negative = 2*ones(1,size(training_x_negative,1));% 反例标签
testing_y_positive = ones(1,size(training_x_positive,1)); % 正例标签
testing_y_negative = 2*ones(1,size(training_x_negative,1));% 反例标签
figure;
plot(training_x_positive(:,1),training_x_positive(:,2),'r+');
hold on;
plot(training_x_negative(:,1),training_x_negative(:,2),'b*');
title('正例(+)和反例(*)');
%
training_x = cat(2,training_x_positive',training_x_negative');
testing_x = cat(2,testing_x_positive',testing_x_negative');
training_y = cat(2,training_y_positive,training_y_negative);
testing_y = cat(2,testing_y_positive,testing_y_negative);
% 样本个数
Sample_num = size(training_x,2);
% 打乱样本
re_rank = randperm(Sample_num);
training_x = training_x(:,re_rank);
training_y = training_y(:,re_rank);
% 混合高斯模型的个数
Gaussian_num = 2;
% % 初始化，采用随机指定每个样本的类别，但是结果不是很好
% for iter_Gaussian = 1:Gaussian_num
% mu(:,iter_Gaussian) = rand(feature_num,1);
% fi(1,iter_Gaussian) = 1/Gaussian_num;
% sigma{1,iter_Gaussian} = eye(feature_num);
% end% for iter_Gaussian
% 初始化采用k-mean聚类
Cluster_sample = training_x;
Cluster_sample_num = Sample_num;
ClusterNum = Gaussian_num;
Max_iter_Num = 1000;
Cluster_index = zeros(Cluster_sample_num,1);% 样本的类别
mu = ones(feature_num, ClusterNum);% 聚类中心
Loss = zeros(Max_iter_Num,1);
Distortion_Loss_Old = 1e10;
Distortion_Loss = 0;
Count = 0;
epsion = 1e-6;
for iter = 1:Max_iter_Num
Count = Count + 1;
if iter == 1
% 随机初始化聚类中心
[cluster_index] = randperm(Cluster_sample_num,2);
mu(:,1) = Cluster_sample(:,cluster_index(1));
mu(:,2) = Cluster_sample(:,cluster_index(2));
cluster_mu_1 = sum((Cluster_sample - repmat(mu(:,1),[1,Cluster_sample_num])).^2,1);
cluster_mu_2 = sum((Cluster_sample - repmat(mu(:,2),[1,Cluster_sample_num])).^2,1);
cluster_one_index = find( (cluster_mu_1 - cluster_mu_2) < 0);
cluster_two_index = find( (cluster_mu_1 - cluster_mu_2) >= 0);
Cluster_index(cluster_one_index) = 1;
Cluster_index(cluster_two_index) = 2;
Distortion_Loss = sum(sum((Cluster_sample(:,cluster_one_index) - repmat(mu(:,1),[1,length(cluster_one_index)])).^2,1)) + sum(sum((Cluster_sample(:,cluster_two_index) - repmat(mu(:,2),[1,length(cluster_two_index)])).^2,1)); % 初始的
mu(:,1) = mean(Cluster_sample(:,cluster_one_index),2); % 更新聚类中心
mu(:,2) = mean(Cluster_sample(:,cluster_two_index),2);
else
cluster_mu_1 = sum((Cluster_sample - repmat(mu(:,1),[1,Cluster_sample_num])).^2,1);
cluster_mu_2 = sum((Cluster_sample - repmat(mu(:,2),[1,Cluster_sample_num])).^2,1);
cluster_one_index = find( (cluster_mu_1 - cluster_mu_2) < 0);
cluster_two_index = find( (cluster_mu_1 - cluster_mu_2) >= 0);
Cluster_index(cluster_one_index) = 1;
Cluster_index(cluster_two_index) = 2;
Distortion_Loss = sum(sum((Cluster_sample(:,cluster_one_index) - repmat(mu(:,1),[1,length(cluster_one_index)])).^2,1)) + sum(sum((Cluster_sample(:,cluster_two_index) - repmat(mu(:,2),[1,length(cluster_two_index)])).^2,1));
mu(:,1) = mean(Cluster_sample(:,cluster_one_index),2); % 更新聚类中心
mu(:,2) = mean(Cluster_sample(:,cluster_two_index),2);
end
Loss(Count) = Distortion_Loss;
if abs(Distortion_Loss_Old - Distortion_Loss) < epsion
break;
end
Distortion_Loss_Old = Distortion_Loss;
end
%显示聚类结果
figure;
for iter = 1:Cluster_sample_num
hold on;
if Cluster_index(iter) == 1
plot3(Cluster_sample(1,iter),Cluster_sample(2,iter),Cluster_sample(3,iter),'g+');
else
plot3(Cluster_sample(1,iter),Cluster_sample(2,iter),Cluster_sample(3,iter),'b*');
end
end
% 显示聚类中心
hold on
plot3(mu(1,1),mu(2,1),mu(3,1),'rx');
hold on
plot3(mu(1,2),mu(2,2),mu(3,2),'rx');
figure;
plot(1:Count,Loss(1:Count,1),'r*');
title('畸变函数收敛曲线');
sigma = cell(1,Gaussian_num);
mu_old = 100000*mu;
epsion = 1e-8;
fi = zeros(1,Gaussian_num);
for iter_Gaussian = 1:Gaussian_num
Cluster_index_k = find(Cluster_index == iter_Gaussian);
fi(1,iter_Gaussian) = length(Cluster_index_k)/Sample_num;
sigma{1,iter_Gaussian} = zeros(feature_num);
for iter_Sample = 1:Sample_num
if Cluster_index(iter_Sample) == iter_Gaussian
sigma{1,iter_Gaussian} = sigma{1,iter_Gaussian} + (training_x(:,iter_Sample) - mu(:,iter_Gaussian))*(training_x(:,iter_Sample) - mu(:,iter_Gaussian))';
end
end
sigma{1,iter_Gaussian} = sigma{1,iter_Gaussian}/length(Cluster_index_k);
end
disp('初始高斯混合模型1参数');
fi(1,1)
mu(:,1)
sigma{1,1}
disp('初始高斯混合模型2参数');
fi(1,2)
mu(:,2)
sigma{1,2}
W = zeros(Gaussian_num,Sample_num);
% 进行迭代更新
iter_Gaussian_Max_num = 1000;
for iter = 1:iter_Gaussian_Max_num
% E-M算法学习混合高斯模型
% E-Step: 计算W_i_j = p(x^i | z^i = j; mu, sigma)*p(z^i = j,fi)/sum_l={1,,,k}(p(x^i | z^i = l; mu, sigma)*p(z^i = l, fi))
for iter_Sample = 1:Sample_num
for iter_Gaussian = 1:Gaussian_num
W(iter_Gaussian,iter_Sample) = 1/((2*pi)^(feature_num/2)*sqrt(det(sigma{1,iter_Gaussian})))*exp(-1/2*(training_x(:,iter_Sample)- mu(:,iter_Gaussian))'/sigma{1,iter_Gaussian}*(training_x(:,iter_Sample) - mu(:,iter_Gaussian)))*fi(1,iter_Gaussian);
end
end
W = W./repmat(sum(W,1),[Gaussian_num,1]);
% M-Step:
for iter_Gaussian =1:Gaussian_num
fi(1,iter_Gaussian) = sum(W(iter_Gaussian,:))/Sample_num;
mu(:,iter_Gaussian) = sum(repmat(W(iter_Gaussian,:),[feature_num,1]).*training_x,2);
mu(:,iter_Gaussian) = mu(:,iter_Gaussian)./sum(W(iter_Gaussian,:));
sigma{1,iter_Gaussian} = zeros(feature_num);
for iter_Sample = 1:Sample_num
sigma{1,iter_Gaussian} = sigma{1,iter_Gaussian} + W(iter_Gaussian,iter_Sample)*(training_x(:,iter_Sample) - mu(:,iter_Gaussian))*(training_x(:,iter_Sample) - mu(:,iter_Gaussian))';
end% for iter_Sample
sigma{1,iter_Gaussian} = sigma{1,iter_Gaussian}./sum(W(iter_Gaussian,:));
end% for iter_Gaussian
if mod(iter,100) == 0
disp(['迭代第',num2str(iter),'次']);
disp('高斯混合模型1参数');
fi(1,1)
mu(:,1)
sigma{1,1}
disp('高斯混合模型2参数');
fi(1,2)
mu(:,2)
sigma{1,2}
end
if sum(sum(abs(mu_old - mu))) < epsion
break;
end
mu_old = mu;
end
disp('学习获得的高斯混合模型1参数');
fi(1,1)
mu(:,1)
sigma{1,1}
disp('学习获得的高斯混合模型2参数');
fi(1,2)
mu(:,2)
sigma{1,2}
结果如下：
初始高斯混合模型1参数
ans =
0.5125
ans =
3.9788
4.9553
5.9916
ans =
0.9897 -0.0364 -0.0609
-0.0364 0.9018 -0.0459
-0.0609 -0.0459 0.9617
初始高斯混合模型2参数
ans =
0.4875
ans =
1.9684
3.0232
3.9805
ans =
1.0290 0.0047 -0.0530
0.0047 0.8828 -0.0519
-0.0530 -0.0519 0.9517
迭代第100次
高斯混合模型1参数
ans =
0.4930
ans =
3.9830
4.9807
6.0065
ans =
1.0221 -0.0224 -0.0337
-0.0224 0.8995 -0.0430
-0.0337 -0.0430 0.9804
高斯混合模型2参数
ans =
0.5070
ans =
2.0418
3.0729
4.0436
ans =
1.1297 0.0778 0.0340
0.0778 0.9309 0.0141
0.0340 0.0141 1.0275
学习获得的高斯混合模型1参数
ans =
0.4930
ans =
3.9830
4.9807
6.0065
ans =
1.0221 -0.0224 -0.0337
-0.0224 0.8995 -0.0430
-0.0337 -0.0430 0.9804
学习获得的高斯混合模型2参数
ans =
0.5070
ans =
2.0418
3.0729
4.0436
ans =
1.1297 0.0778 0.0340
0.0778 0.9309 0.0141
0.0340 0.0141 1.0275
内容来自用户分享和网络整理，不保证内容的准确性，如有侵权内容，可联系管理员处理 
相关文章推荐
- （EM算法）The EM Algorithm
- （EM算法）The EM Algorithm
- （EM算法）The EM Algorithm
- （EM算法）The EM Algorithm
- （EM算法）The EM Algorithm
- The EM Algorithm（EM算法）
- Weka:call for the EM algorithm to achieve clustering.（EM算法）
- （EM算法）The EM Algorithm
- （EM算法）The EM Algorithm
- （EM算法）The EM Algorithm
- （EM算法）The EM Algorithm
- （EM算法）The EM Algorithm
- （EM算法）The EM Algorithm
- （EM算法）The EM Algorithm
- EM算法[The EM Algorithm]
- （EM算法）The EM Algorithm
- 【机器学习-斯坦福】学习笔记12 （EM算法）The EM Algorithm
- （EM算法）The EM Algorithm
- （EM算法）The EM Algorithm
- （EM算法）The EM Algorithm