高斯混合模型聚类算法

原文地址:http://blog.csdn.net/hjimce/article/details/45244603

作者:hjimce

高斯混合算法是EM算法的一个典型的应用,EM算法的推导过程这里不打算详解,直接讲GMM算法的实现。之前做图像分割grab cut 算法的时候,只知道把opencv中的高斯混合模型代码复制下来,然后封装成类使用,学的比较浅。结果没过几天发现高斯混合算法又忘了差不多了,于是用matlab去亲自写过一遍,终于发现了高斯混合模型的奥义。我的理解是高斯混合模型其实是进化版的k均值算法,因此学习高斯混合模型,最好还是把k均值算法写过一遍。高斯混合与k均值的本质区别在于权值问题,k均值采用的是均匀权值,而高斯混合的权值需要根据高斯模型的概率进行确定。

开始学习高斯混合模型,需要先简单复习一下单高斯模型的参数估计方法,描述一个高斯模型其实就是要计算它的均值、协方差矩阵(一维空间为方差,二维以上称之为协方差矩阵):

假设有数据集X={x1,x2,x3……,xn},那么用这些数据来估计单高斯模型参数的计算公式为

OK,开始写代码前,先用matlab生成数据集,然后在进行聚类:

利用matlab的生成高斯模型数据集X:

mu = [2 3]; SIGMA = [1 0; 0 2]; r1 = mvnrnd(mu,SIGMA,1000); plot(r1(:,1),r1(:,2),'r+');

然后利用上面的估计方法计算均值,和协方差是否满足均值为[2 3],协方差为[1 0; 0 2];测试代码如下,r2、covmat即为计算结果

[m n]=size(r1); center=sum(r1)./m; r2(:,1)=r1(:,1)-center(1); r2(:,2)=r1(:,2)-center(2); covmat=1/m*r2'*r2;

先把单高斯模型的函数写好,因为高斯混合模型是它的进化版,计算高斯混合模型过程中需要调用单高斯模型参数估计,写好代码后面才不会乱掉。开始高斯混合建模之前,我先用matlab生成一个测试数据集data,如下图,然后再进行算法测试。

生成数据集代码如下:

%生成测试数据 mu = [2 3];%测试数据1 SIGMA = [1 0; 0 2]; r1 = mvnrnd(mu,SIGMA,100); plot(r1(:,1),r1(:,2),'.'); hold on; mu = [10 10];%测试数据2 SIGMA = [ 1 0; 0 2]; r2 = mvnrnd(mu,SIGMA,100); plot(r2(:,1),r2(:,2),'.'); mu = [5 8];%测试数据3 SIGMA = [ 1 0; 0 2]; r3= mvnrnd(mu,SIGMA,100); plot(r3(:,1),r3(:,2),'.'); data=[r1;r2;r3];

ok,数据生成完毕,接着我们正式开始高斯混合算法解析,先看一下高斯混合模型的建模求参步骤:

高斯混合模型的求解,说得简单一点就是要求解高斯模型中的均值与协方差,现在我们要把上述的数据分成3类,那么我们就是要求解3个均值及其对应的3个协方差矩阵。先讲一下总体步骤,高斯混合模型包含3个步骤:

a.初始化各个高斯模型的参数,及每个高斯模型的权重;

b.根据各个高斯模型的参数及其权重,计算每个点属于各个高斯模型的权重,计算公式为:

其中:,Wj是每个高斯模型在这个模型所占用得权重。这个公式说的简单一点就是每个高斯模型的权重与其概率的乘积,这样计算出来就相当于每个高斯模型在每个数据点中的所占用的比例。

c.更新各个高斯模型的均值与方差,计算公式如下:

d.更新各个高斯模型的总权重,计算公式如下:

其实第c、d两个步骤,无所谓顺序,你完全可以总权重更新放在各个模型参数更新之前。迭代过程就是b、c、d三个步骤进行更新就可以了。OK,接着结合上面的公式写一写代码。

(1)初始化高斯模型参数。

这一步初始化,在实际应用中一般是先通过k均值算法进行初始聚类,然后根据聚类结果进行计算初始化参数。不过这里我为了测试,我们选择随机初始化,这样才能看出GMM算法到底能不能实现聚类。

我这里各个高斯模型初始均值(中心)的初始化方法选择跟k均值的初始化方法一样,也就是随机选择k个点位置作为k个高斯模型的初始均值。然后协方差矩阵的初始化,我选择单位矩阵,具体代码如下:

[m n]=size(data); kn=3; countflag=zeros(1,kn); tdata=cell(1,kn);%建立3个空矩阵 mu=cell(1,kn);%建立3个空矩阵 sigma=cell(1,kn);%建立3个空矩阵 %方案2 随机初始化参数 for i=1:knmu{1,i}=data(i*10,:);sigma{1,i}=eye(2,2);weightp(i)=1/kn; end

(2)计算各个模型在各个点的权重值

这一步是计算每个数据点属于各个高斯混合的概率,说白了就是计算权值:

pro_ij=zeros(m,kn);%存储每个点属于每个类的概率for i=1:msumpk=0;for j=1:knpk(j)=weightp(j)*GSMPro(mu{1,j},sigma{1,j},data(i,:));sumpk=sumpk+pk(j);endfor j=1:knpro_ij(i,j)=pk(j)/sumpk;endend

(3)步骤c 更新参数

for j=1:kn[mu{1,j},sigma{1,j}]=WeightGSM(data,pro_ij(:,j)); end

(4)步骤d 更新各个模型的总权重

for j=1:knweightp(j)=sum(pro_ij(:,j))/m;end

然后把步骤2、3、4的代码放在循环语句中进行迭代就ok了。最后贴一下整份代码:

1、脚本文件:

close all; clear; clc; %生成测试数据 mu = [2 3];%测试数据1 SIGMA = [1 0; 0 2]; r1 = mvnrnd(mu,SIGMA,100); plot(r1(:,1),r1(:,2),'.'); hold on; mu = [10 10];%测试数据2 SIGMA = [ 1 0; 0 2]; r2 = mvnrnd(mu,SIGMA,100); plot(r2(:,1),r2(:,2),'.');mu = [5 8];%测试数据3 SIGMA = [ 1 0; 0 2]; r3= mvnrnd(mu,SIGMA,100); plot(r3(:,1),r3(:,2),'.');data=[r1;r2;r3];[m n]=size(data); kn=3; countflag=zeros(1,kn); tdata=cell(1,kn);%建立10个空矩阵 mu=cell(1,kn);%建立10个空矩阵 sigma=cell(1,kn);%建立10个空矩阵 % 方案1 初始化采用kmeans,做参数的初步估计 % Idx=kmeans(data,kn); % figure(2);%绘制初始化结果 % hold on; % for i=1:m % if Idx(i)==1 % plot(data(i,1),data(i,2),'.y'); % elseif Idx(i)==2 % plot(data(i,1),data(i,2),'.b'); % end % end % for i=1:m % tdata{1,Idx(i)}=[tdata{1,Idx(i)};data(i,:)]; % end % for i=1:kn % [mu{1,i},sigma{1,i}]=GSMData(tdata{1,i}); % end % for i=1:kn % [trow,tcol]=size(tdata{1,i}); % weightp(i)=trow/m; % end %方案2 随机初始化 for i=1:knmu{1,i}=data(i*10,:);sigma{1,i}=eye(2,2);weightp(i)=1/kn; endit=1;while it<1000%E步 计算每个点处于每个类的概率pro_ij=zeros(m,kn);%存储每个点属于每个类的概率for i=1:msumpk=0;for j=1:knpk(j)=weightp(j)*GSMPro(mu{1,j},sigma{1,j},data(i,:));sumpk=sumpk+pk(j);endfor j=1:knpro_ij(i,j)=pk(j)/sumpk;endend %M步 for j=1:kn[mu{1,j},sigma{1,j}]=WeightGSM(data,pro_ij(:,j)); end%更新权值for j=1:knweightp(j)=sum(pro_ij(:,j))/m;endsumw=sum(weightp);it=it+1; end for i=1:m[value index]=max(pro_ij(i,:));Idx(i)=index; end figure(2); hold on; for i=1:mif Idx(i)==1plot(data(i,1),data(i,2),'.y');elseif Idx(i)==2plot(data(i,1),data(i,2),'.b');elseif Idx(i)==3plot(data(i,1),data(i,2),'.r');end end% figure(3); % %px=gmmstd(data,3); % for i=1:m % [value index]=max(px(i,:)); % Idx(i)=index; % end % hold on; % for i=1:m % if Idx(i)==1 % plot(data(i,1),data(i,2),'.y'); % elseif Idx(i)==2 % plot(data(i,1),data(i,2),'.b'); % elseif Idx(i)==3 % plot(data(i,1),data(i,2),'.r'); % end % end %单高斯模型参数估计 % [m n]=size(r1); % center=sum(r1)./m; % r2(:,1)=r1(:,1)-center(1); % r2(:,2)=r1(:,2)-center(2); % covmat=1/m*r2'*r2;

2、相关函数

function [ mu ,sigma ] = WeightGSM(data,weight)%计算加权均值[m n]=size(data);sumweight=sum(weight);weightdata=[];for i=1:mweightdata(i,:)=weight(i)*data(i,:);endcenter=sum(weightdata)/sumweight;%计算加权协方差for i=1:nr2(:,i)=data(:,i)-center(i);endfor i=1:mr1(i,:)=weight(i)*r2(i,:);endsigma=1/sumweight*r1'*r2;mu=center; end function [pro] = GSMPro(mu ,sigma,x)pro=exp(-0.5*(x-mu)*inv(sigma)*(x-mu)');pro=1/sqrt(2*pi*det(sigma))*pro; end

看以下最后的测试结果:

*******************作者:hjimce     联系qq:1393852684 更多资源请关注我的博客:http://blog.csdn.net/hjimce                  原创文章,转载请注明出处 *******************

机器学习(四)高斯混合模型相关推荐

  1. 【机器学习】高斯混合模型详解

    目录 1 引言 2 高斯混合模型 2.1 高斯分布 2.2 高斯混合模型 3 高斯混合模型的求解 4 参考文献 1 引言   高斯混合模型(Gaussian Mixture Model, GMM)是单 ...

  2. 机器学习 - GMM高斯混合模型

    高斯混合模型:一个数据集可以由1 - n个高斯模型加权求和来生成.采用概率模型来表达数据分布. 由于涉及的符号太多,因此算法详细推导过程在这份PDF里面: http://download.csdn.n ...

  3. 【机器学习之高斯混合模型(Gaussian Mixed Model,GMM) 】

    文章目录 前言 一.高斯混合模型(Gaussian Mixed Model,GMM) 是什么? 二.详解GMM 1.初步原理 2.EM算法 3.深读原理 3.GMM(高斯混合模型)和K-means的比 ...

  4. 机器学习 : 高斯混合模型及EM算法

    Mixtures of Gaussian 这一讲,我们讨论利用EM (Expectation-Maximization)做概率密度的估计.假设我们有一组训练样本 x(1),x(2),...x(m) { ...

  5. 百面机器学习—7.K均值算法、EM算法与高斯混合模型要点总结

    文章目录 一.总结K均值算法步骤 二.如何合理选择K值? 三.K均值算法的优缺点是什么? 四.如何对K均值算法进行调优? 五.EM算法解决什么问题? 六.EM算法流程是什么? 六.EM算法能保证收敛嘛 ...

  6. 【机器学习笔记11】高斯混合模型(GMM)【上篇】原理与推导

    文章目录 推荐阅读 前言 高斯混合模型简介 GMM与K-mean 高斯混合模型的概率密度函数 几何角度 混合模型角度 可能会弄混的地方 隐变量的分布与隐变量的后验概率分布 极大似然估计 EM算法求近似 ...

  7. 机器学习笔记之高斯混合模型(一)模型介绍

    机器学习笔记之高斯混合模型--模型介绍 引言 高斯混合模型介绍 示例介绍 从几何角度观察高斯混合模型 从混合模型的角度观察 概率混合模型的引出 从概率生成模型的角度观察高斯混合模型 引言 上一系列介绍 ...

  8. 机器学习教程 之 EM算法 :高斯混合模型聚类算法 (python基于《统计学习方法》实现,附数据集和代码)

    之前写过一篇博客讲述极大似然方法, 这一方法通常适用于知道观测数据 Y Y Y,求解模型参数 θ \theta θ的场合,即 P ( Y ∣ θ ) P(Y|\theta) P(Y∣θ). 但是,在更 ...

  9. 机器学习基础专题:高斯混合模型和最大期望EM算法以及代码实现

    高斯混合模型 混合模型是潜变量模型的一种,是最常见的形式之一.而高斯混合模型(Gaussian Mixture Models, GMM)是混合模型中最常见的一种.zzz代表该数据点是由某一个高斯分布产 ...

最新文章

  1. java项目性能测试过程记录
  2. 快速幂(二进制,十进制)
  3. JAVA笔记(运算符)
  4. oracle10 监听日志,windows 清空oracle的监听日志listener.log
  5. 【算法设计与分析】最长公共子序列问题 动态规划算法 超详细
  6. mysql录入foreigen错误_编译安装MySQL5.6失败的相关问题解决方案
  7. li span 时间向右排版问题
  8. java图形界面颜色随机变换,JavaScript实现鼠标移入随机变换颜色
  9. 9.25-CSS样式以及结构布局
  10. Cocos2d BMFont解析
  11. 轨迹压缩文献阅读: Similarity-Based Compression of GPS Trajectory Data
  12. TFRecord —— tensorflow 下的统一数据存储格式
  13. ufw防火墙配置命令
  14. 怎样手工清除autorun病毒
  15. Python接口自动化
  16. ORA-00257错误解决方法
  17. 删除页眉页脚中横线的方法
  18. 简谈浅层拷贝和深层拷贝
  19. AD生成BOM表/元器件表
  20. 传感器是指纹识别产品的数据入口

热门文章

  1. python comprehension_python list comprehension在一次迭代中产生两个值
  2. linux配置本地yum(CentOS7)
  3. oracle打开文件模式无效,oracle expdp导入时 提示“ORA-39002: 操作无效 ORA-39070: 无法打开日志文件 ”...
  4. python socket多线程 获取朋友列表_python中的(多线程)套接字列表/数组
  5. 第二章 寄存器基础概念
  6. python-9:nonlocal,指定上一级变量
  7. Maven常用的命令
  8. Javascript IE 内存释放
  9. Jquery全选单选功能
  10. 关于指向堆的指针内涵