程序博客网 > 设计一款人工智能产品

GMM的EM算法实现

来源：互联网 发布：设计一款人工智能产品 编辑：程序博客网 时间：2024/05/18 01:31

在 聚类算法K-Means, K-Medoids, GMM, Spectral clustering，Ncut一文中我们给出了GMM算法的基本模型与似然函数，在EM算法原理中对EM算法的实现与收敛性证明进行了详细说明。本文主要针对如何用EM算法在混合高斯模型下进行聚类进行代码上的分析说明。


1. GMM模型：

每个 GMM 由 K 个 Gaussian 分布组成，每个 Gaussian 称为一个“Component”，这些 Component 线性加成在一起就组成了 GMM 的概率密度函数：


根据上面的式子，如果我们要从 GMM 的分布中随机地取一个点的话，实际上可以分为两步：首先随机地在这 K个Gaussian Component 之中选一个，每个 Component 被选中的概率实际上就是它的系数 pi(k) ，选中了 Component 之后，再单独地考虑从这个 Component 的分布中选取一个点就可以了──这里已经回到了普通的 Gaussian 分布，转化为了已知的问题。

那么如何用 GMM 来做 clustering 呢？其实很简单，现在我们有了数据，假定它们是由 GMM 生成出来的，那么我们只要根据数据推出 GMM 的概率分布来就可以了，然后 GMM 的 K 个 Component 实际上就对应了 K 个 cluster 了。根据数据来推算概率密度通常被称作 density estimation ，特别地，当我们在已知（或假定）了概率密度函数的形式，而要估计其中的参数的过程被称作“参数估计”。


2. 参数与似然函数：

现在假设我们有 N 个数据点，并假设它们服从某个分布（记作 p(x) ），现在要确定里面的一些参数的值，例如，在 GMM 中，我们就需要确定 影响因子pi(k)、各类均值pMiu(k) 和 各类协方差pSigma(k) 这些参数。 我们的想法是，找到这样一组参数，它所确定的概率分布生成这些给定的数据点的概率最大，而这个概率实际上就等于  ，我们把这个乘积称作似然函数 (Likelihood Function)。通常单个点的概率都很小，许多很小的数字相乘起来在计算机里很容易造成浮点数下溢，因此我们通常会对其取对数，把乘积变成加和 \sum_{i=1}^N \log p(x_i)，得到 log-likelihood function 。接下来我们只要将这个函数最大化（通常的做法是求导并令导数等于零，然后解方程），亦即找到这样一组参数值，它让似然函数取得最大值，我们就认为这是最合适的参数，这样就完成了参数估计的过程。

下面让我们来看一看 GMM 的 log-likelihood function ：


由于在对数函数里面又有加和，我们没法直接用求导解方程的办法直接求得最大值。为了解决这个问题，我们采取之前从 GMM 中随机选点的办法：分成两步，实际上也就类似于K-means 的两步。



3. 算法流程：

1.  估计数据由每个 Component 生成的概率（并不是每个 Component 被选中的概率）：对于每个数据 x_i 来说，它由第 k 个 Component 生成的概率为


其中N（xi | μk,Σk）就是后验概率。


2. 通过极大似然估计可以通过求到令参数=0得到参数pMiu，pSigma的值。具体请见这篇文章第三部分。


其中 N_k = \sum_{i=1}^N \gamma(i, k) ，并且 \pi_k 也顺理成章地可以估计为 N_k/N 。


3. 重复迭代前面两步，直到似然函数的值收敛为止。



4. matlab实现GMM聚类代码与解释：


说明：fea为训练样本数据，gnd为样本标号。算法中的思想和上面写的一模一样，在最后的判断accuracy方面，由于聚类和分类不同，只是得到一些 cluster ，而并不知道这些 cluster 应该被打上什么标签，或者说。由于我们的目的是衡量聚类算法的 performance ，因此直接假定这一步能实现最优的对应关系，将每个 cluster 对应到一类上去。一种办法是枚举所有可能的情况并选出最优解，另外，对于这样的问题，我们还可以用 Hungarian algorithm 来求解。具体的Hungarian代码我放在了资源里，调用方法已经写在下面函数中了。


注意：资源里我放的是Kmeans的代码，大家下载的时候只要用bestMap.m等几个文件就好~


1. gmm.m，最核心的函数，进行模型与参数确定。

[cpp] view plaincopyprint?
  1. function varargout = gmm(X, K_or_centroids)  
  2. % ============================================================  
  3. % Expectation-Maximization iteration implementation of  
  4. % Gaussian Mixture Model.  
  5. %  
  6. % PX = GMM(X, K_OR_CENTROIDS)  
  7. % [PX MODEL] = GMM(X, K_OR_CENTROIDS)  
  8. %  
  9. %  - X: N-by-D data matrix.  
  10. %  - K_OR_CENTROIDS: either K indicating the number of  
  11. %       components or a K-by-D matrix indicating the  
  12. %       choosing of the initial K centroids.  
  13. %  
  14. %  - PX: N-by-K matrix indicating the probability of each  
  15. %       component generating each point.  
  16. %  - MODEL: a structure containing the parameters for a GMM:  
  17. %       MODEL.Miu: a K-by-D matrix.  
  18. %       MODEL.Sigma: a D-by-D-by-K matrix.  
  19. %       MODEL.Pi: a 1-by-K vector.  
  20. % ============================================================  
  21. % @SourceCode Author: Pluskid (http://blog.pluskid.org)  
  22. % @Appended by : Sophia_qing (http://blog.csdn.net/abcjennifer)  
  23.       
  24.   
  25. %% Generate Initial Centroids  
  26.     threshold = 1e-15;  
  27.     [N, D] = size(X);  
  28.    
  29.     if isscalar(K_or_centroids) %if K_or_centroid is a 1*1 number  
  30.         K = K_or_centroids;  
  31.         Rn_index = randperm(N); %random index N samples  
  32.         centroids = X(Rn_index(1:K), :); %generate K random centroid  
  33.     else % K_or_centroid is a initial K centroid  
  34.         K = size(K_or_centroids, 1);   
  35.         centroids = K_or_centroids;  
  36.     end  
  37.    
  38.     %% initial values  
  39.     [pMiu pPi pSigma] = init_params();  
  40.    
  41.     Lprev = -inf; %上一次聚类的误差  
  42.       
  43.     %% EM Algorithm  
  44.     while true  
  45.         %% Estimation Step  
  46.         Px = calc_prob();  
  47.    
  48.         % new value for pGamma(N*k), pGamma(i,k) = Xi由第k个Gaussian生成的概率  
  49.         % 或者说xi中有pGamma(i,k)是由第k个Gaussian生成的  
  50.         pGamma = Px .* repmat(pPi, N, 1); %分子 = pi(k) * N(xi | pMiu(k), pSigma(k))  
  51.         pGamma = pGamma ./ repmat(sum(pGamma, 2), 1, K); %分母 = pi(j) * N(xi | pMiu(j), pSigma(j))对所有j求和  
  52.    
  53.         %% Maximization Step - through Maximize likelihood Estimation  
  54.           
  55.         Nk = sum(pGamma, 1); %Nk(1*k) = 第k个高斯生成每个样本的概率的和，所有Nk的总和为N。  
  56.           
  57.         % update pMiu  
  58.         pMiu = diag(1./Nk) * pGamma' * X; %update pMiu through MLE(通过令导数 = 0得到)  
  59.         pPi = Nk/N;  
  60.           
  61.         % update k个 pSigma  
  62.         for kk = 1:K   
  63.             Xshift = X-repmat(pMiu(kk, :), N, 1);  
  64.             pSigma(:, :, kk) = (Xshift' * ...  
  65.                 (diag(pGamma(:, kk)) * Xshift)) / Nk(kk);  
  66.         end  
  67.    
  68.         % check for convergence  
  69.         L = sum(log(Px*pPi'));  
  70.         if L-Lprev < threshold  
  71.             break;  
  72.         end  
  73.         Lprev = L;  
  74.     end  
  75.    
  76.     if nargout == 1  
  77.         varargout = {Px};  
  78.     else  
  79.         model = [];  
  80.         model.Miu = pMiu;  
  81.         model.Sigma = pSigma;  
  82.         model.Pi = pPi;  
  83.         varargout = {Px, model};  
  84.     end  
  85.    
  86.     %% Function Definition  
  87.       
  88.     function [pMiu pPi pSigma] = init_params()  
  89.         pMiu = centroids; %k*D, 即k类的中心点  
  90.         pPi = zeros(1, K); %k类GMM所占权重（influence factor）  
  91.         pSigma = zeros(D, D, K); %k类GMM的协方差矩阵，每个是D*D的  
  92.    
  93.         % 距离矩阵，计算N*K的矩阵（x-pMiu）^2 = x^2+pMiu^2-2*x*Miu  
  94.         distmat = repmat(sum(X.*X, 2), 1, K) + ... %x^2, N*1的矩阵replicateK列  
  95.             repmat(sum(pMiu.*pMiu, 2)', N, 1) - ...%pMiu^2，1*K的矩阵replicateN行  
  96.             2*X*pMiu';  
  97.         [~, labels] = min(distmat, [], 2);%Return the minimum from each row  
  98.    
  99.         for k=1:K  
  100.             Xk = X(labels == k, :);  
  101.             pPi(k) = size(Xk, 1)/N;  
  102.             pSigma(:, :, k) = cov(Xk);  
  103.         end  
  104.     end  
  105.    
  106.     function Px = calc_prob()   
  107.         %Gaussian posterior probability   
  108.         %N(x|pMiu,pSigma) = 1/((2pi)^(D/2))*(1/(abs(sigma))^0.5)*exp(-1/2*(x-pMiu)'pSigma^(-1)*(x-pMiu))  
  109.         Px = zeros(N, K);  
  110.         for k = 1:K  
  111.             Xshift = X-repmat(pMiu(k, :), N, 1); %X-pMiu  
  112.             inv_pSigma = inv(pSigma(:, :, k));  
  113.             tmp = sum((Xshift*inv_pSigma) .* Xshift, 2);  
  114.             coef = (2*pi)^(-D/2) * sqrt(det(inv_pSigma));  
  115.             Px(:, k) = coef * exp(-0.5*tmp);  
  116.         end  
  117.     end  
  118. end  


2. gmm_accuracy.m调用gmm.m，计算准确率：

[cpp] view plaincopyprint?
  1. function [ Accuracy ] = gmm_accuracy( Data_fea, gnd_label, K )  
  2. %Calculate the accuracy Clustered by GMM model  
  3.   
  4. px = gmm(Data_fea,K);  
  5. [~, cls_ind] = max(px,[],1); %cls_ind = cluster label   
  6. Accuracy = cal_accuracy(cls_ind, gnd_label);  
  7.   
  8.     function [acc] = cal_accuracy(gnd,estimate_label)  
  9.         res = bestMap(gnd,estimate_label);  
  10.         acc = length(find(gnd == res))/length(gnd);  
  11.     end  
  12.   
  13. end  


3. 主函数调用

gmm_acc = gmm_accuracy(fea,gnd,N_classes);








写了本文进行总结后自己很受益，也希望大家可以好好YM下上面pluskid的gmm.m，不光是算法，其中的矩阵处理代码也写的很简洁，很值得学习。

另外看了两份东西非常受益，一个是pluskid大牛的《漫谈 Clustering (3): Gaussian Mixture Model》，一个是JerryLead的EM算法详解，大家有兴趣也可以看一下，写的很好。

 

http://blog.csdn.net/abcjennifer/article/details/8198352


设计一款人工智能产品 设计一款人工智能产品
原创粉丝点击
热门IT博客
香港虚拟主机 知乎香港虚拟主机 知乎
pageadmin 源码pageadmin 源码
男朋友礼物 知乎
朱雀投资怎么样知乎
淘宝心级怎么算
夏老师c语言视频教程
凌波多媒体教室软件
手机直播源码教程
画江湖之换世门生知乎
55开淘宝店地址
2017全球社交网络排名
统计局数据报送
淘宝流量入口有哪些
奔跑吧兄弟4网络直播
设计模型软件下载
风月知相思蓝语凌云
快速傅里叶变换算法
童自荣事件知乎
袜子推荐 知乎
陕西广电网络是国企吗
热门问题 老师的惩罚 人脸识别 我在镇武司摸鱼那些年 重生之率土为王 我在大康的咸鱼生活 盘龙之生命进化 天生仙种 凡人之先天五行 春回大明朝 姑娘不必设防，我是瞎子 老旧小区改造下水一楼不同意怎么办 老旧小区下水改造没改怎么办 替公司租房子中介不退押金怎么办 想在昆山找合租房的该怎么办 链家二手房价钱买贵了怎么办 拿私人房产证抵押借钱不还怎么办 在借贷宝里借钱不还怎么办 出租屋的大门感应钥匙弄丢了怎么办 法院拍卖的房子房主不配合怎么办 租的房子如果房主卖了怎么办 房东把门锁换了里面的东西怎么办 房租没到期房东把门锁换了怎么办 租了三年店面房东违反了合同怎么办 学生登录教务系统的密码忘记怎么办 铜陵无牌助力车被交警查到怎么办 福州超标电动车被交警抓到怎么办 单位自管公租房承租人去世怎么办 取得房产证后贷款批不下来怎么办 租店面遇到难搞的房东怎么办 在拆违通知书上签字了该怎么办 单位没交公职金的退休后怎么办 公租房合同到期后没有续签怎么办 租房合同没到期不想租了怎么办 租的房子是人家公租房怎么办 五险合一软件口令忘记了怎么办 计生办把婚育状况登记错了怎么办? 医院发票法院要保险也要怎么办 上海社保里生育险暂停参保怎么办? 痔疮手术后大姨妈来了怎么办 微创痔疮术后第五天涨出血怎么办? 肚子胀疼大便拉不出来怎么办 得痔疮了该怎么办昆明东大治 下体痒还没去检查就来月经了怎么办 直肠造口手术后造口肠子突出怎么办 痔疮pph手术瘢痕两年了该怎么办 肛瘘挂线术后六天腹泻了怎么办 刚满月的孩子鼻子不通气怎么办 齐鲁医院挂的号晚了怎么办 手机微信安装后注册失败怎么办 舞蹈基本功胸怎么都转不动怎么办 饥荒手机版第10天遇到的狗怎么办

程序博客网，程序员的互联网技术博客家园。csdn论坛精品 msdn技术资料都在这里