上一篇文章描述了EM算法的原理和详细推导过程。本文着重阐述如何用EM算法来求解高斯混合模型。高斯混合模型依问题的随机变量是一维还是多维分为单变量高斯混合模型和多变量高斯混合模型。这一篇主要讨论单变量高斯混合模型。一、高斯混合模型简介高斯混合模型GMM也叫混合高斯模型是一种‌参数化的概率统计方法‌核心思想是假设复杂数据分布由多个高斯分布钟形正态曲线叠加组成能灵活拟合多峰、非球形的复杂数据分布。一般的GMM由多个高斯分布线性组合而成每个分布拥有独立的‌均值或均值向量‌决定分布中心、‌标准差或协方差矩阵‌决定分布的形状、分散程度和‌混合权重‌该分布在总数据中的占比所有权重之和为1。所谓求解或训练一个GMM模型就是在给定一个样本集的情况下用一定的方法来估计GMM中的参数(均值、标准差和混合权重) 通常的求解方法是EM算法[1]‌。EM算法先计算每个数据点属于各分布的后验概率E步再基于该概率更新所有高斯分布的参数M步反复迭代直到模型收敛。在机器学习中GMM模型主要用于非监督学习过程属于‌软聚类‌方法。不像K-means强制把数据点归为某一个簇而是给出每个点属于不同簇的概率支持模糊分类。可适配椭圆形等非球形的簇形状比K-means更贴合真实场景的复杂数据分布。GMM的主要应用场景有1.计算机视觉动态背景下的运动目标检测、图像分割、视频火焰识别)2. 语音领域(建模语音信号特征分布用于区分不同音素或说话人); 3. 数据处理(异常流量检测、金融客户分群、复杂概率密度估计等场景)。GMM模型的核心数学公式是多个高斯分布的加权求和用来描述数据的整体分布情况 其一般形式如下(按随机变量是一维或多维分为单变量GMM和多变量GMM)单变量GMM,这里和分别是第j个单变量高斯分布的均值和标准差。多变量GMM这里和分别是第j个多变量高斯分布的均值向量和协方差矩阵。其中参数或是待估参数。二、EM算法求解单变量高斯混合模型假设有一个样本集, 这些数据看起来有M个峰类于是可以用如下高斯混合模型拟合这些数据则似然函数为对应的对数似然函数是如果直接对上式的参数(权重均值及标准差)求偏导并令其为 0由于对数与求和无法交换顺序导数方程中参数相互耦合‌很难得到参数的显式解析解。如果采用EM 算法通过引入‌适当的隐变量‌表示样本属于哪个高斯分量可以将复杂的优化问题分解为两个可解步骤E 步和M 步)使难解问题转化为一系列易解的子问题‌保证每次迭代后对数似然函数单调递增直至收敛到局部最优解。为此引入隐变量, 它表示当前样本数据是属于哪个类其所有的取值为:, 其中表示当前数据是属于第个高斯分布。设在EM算法中的参数估计已计算出的情况下Q函数为记隐变量的后验概率为利用贝叶斯公式:由联合密度函数与边缘密度函数条件概率密度函数的关系有于是Q函数为因为正态分布密度函数为Q函数进一步整理为以下求函数的最值点。对于参数因为有约束条件由拉格朗日乘数法构造如下拉格朗日函数。求关于参数的偏导令则在上式的基础上再考虑有于是 因此得到参数的更新值对于参数, 求关于参数的偏导令则因此得到参数的更新值为对于参数的估计问题为便于计算这里直接考虑的估计问题。于是记。求关于参数的偏导令则因此得到参数的更新值为于是可以得到GMM模型单变量情况下的EM算法的基本步骤-----------------------------------单变量GMM模型EM算法-----------------------------------------------------------1 初始化模型参数2 对于每次迭代 t 0, 1, 2, 3, ... 直至收敛(2.1) E-step固定当前模型参数计算隐变量 z 的后验概率:(2.2) M-step基于E-step的结果 用下式更新模型参数3当参数更新的变化小于预设阈值或达到最大迭代次数时算法停止输出最终的参数估计。-------------------------------------------------------------------------------------------------------------------------------三、两个简单例子先给一个验证性的例子例 1设一个单变量高斯混合模型由三个高斯模型混合而成:假设其真实的参数为;;。其权重为. 由以上模型生成一组数据试编写程序用EM算法拟合以上模型。解MATLAB自带函数fitgmdist用来拟合一个高斯混合模型以及计算隐变量后验概率的函数posterior。这里为了解释算法的细节我们用自己编写的函数来实现来实现高斯混合模型的拟合过程。先给出高斯混合模型的函数function YuniGMMpdf(X,Alpha,Mu,Sigma) %单变量高斯混合模型 %********************************************************* % 输入: % X N*1 样本集(N个样本) % Mu K*1 数组 (每一行为某类的均值,K为高斯分量数) % Sigma K*1 数组 (每一层为一类的协方差矩阵) % Alpha K*1列向量 (每一个数值为一类的权重(占比)) % 输出: % Y N*1 N个样本点对应的概率密度值 %********************************************************** [K,d]size(Mu); Nsize(X,1); Yzeros(N,K); if d1 for k1:K Y(:,k)uniNormpdf(X,Mu(k,:),Sigma(k,:)); end YY*Alpha; else disp(This is univariate GMM!); end end再用拒绝采样算法为高斯混合模型编写一个采样函数这里用区间[a,b]上的均匀分布作为辅助分布function samplesSamplingUniGMM(numberSample,Alpha,Mu,Sigma) %高斯混合模型的样本生成函数这里用拒绝采样方法实现 %这里选择一个均匀分布作为辅助分布。 %********************************************************* %输入 % numberSample: 样本数量 % Alpha: 权重(K*1) % Mu: 均值(K*1) % Sigma 标准差(K*1) %输出 % samples: 样本集 %********************************************************* %确定分布的范围[a,b] bmax(Mu)(max(Sigma))^2/2; % 均匀分布的下限 amin(Mu)-(max(Sigma))^2/2; % 均匀分布的上限 % 样本集 samples zeros(numberSample,1); %存放样本 %一个合适的M M(b-a)/(sqrt(2*pi)*min(Sigma)); accepted0; for i1:numberSample while 1 %从辅助分布采样 xa(b-a)*rand(); % 计算接受概率 gammauniGMMpdf(x,Alpha,Mu,Sigma)/(M*uniformpdf(x,a,b)); %生成一个[0,1]上均匀分布随机数 urand(); if ugamma %接受这个样本 accepted accepted 1; samples(accepted,1)x; break; % 接受后退出循环继续下一个样本的生成 end end end end function pdfuniformpdf(x,a,b) %区间[a,b]上的均匀分布概率密度函数 %********************************************************* % 输入: % x 1*1 样本 % a 1*1 左端点 % b 1*1 右端点 % 输出: % pdf 1*1 概率密度值 %********************************************************** if xa xb pdf1/(b-a); else pdf0; end end然后给出用EM算法训练高斯混合模型的函数function [Alpha,Mu,Sigma]uniEM_GMM(X,K) %用EM算法训练单变量高斯混合模型 % ********************************************************* % 输入: % X N*1数组(1维点坐标集) % K 数值(划分类别数量) %输出: % Mu K*1 数组 (每一行为某类的坐标中心) % Sigma K*1 数组 (每一层为一类的协方差矩阵) % Alpha K*1列向量 (每一个数值为一类的权重(占比)) % ********************************************************* [N,d]size(X); % N:元素个数, d:维数 %一开始设置每一类有相同协方差矩阵和权重 Alpha_est(1:K,:)1/K; Sigma_est(1:K,:)repmat(cov(X),[K,1]); %依据各维度的最大最小值构建参数的初始值 tMinmin(X); tMaxmax(X); Mu_estrepmat(tMin,K,1)linspace(0,1,K)*(tMax-tMin); MaxIter1e5; for step1:MaxIter MuMu_est; SigmaSigma_est; AlphaAlpha_est; %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %E-Step: % %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %计算后验概率 WuniPosteriorGMM(X,Alpha_est,Mu_est,Sigma_est); %%%%%%%%%%%%%%%%%%%%% %M-Step: % %%%%%%%%%%%%%%%%%%%%% % Alpha Swsum(W,1); Alpha_estSw/N; % Mu Mu_est(X*W./Sw); % Sigma YX-repmat(Mu_est,size(X,1),1); YY.*Y; Sigma_est(sum(W.*Y,1)./Sw); Sigma_estsqrt(Sigma_est); % 收敛条件(计算均方根误差) r1sum((Mu_est-Mu).^2,all); r2sum((Sigma_est-Sigma).^2,all); r3sum((Alpha_est-Alpha).^2,all); Rmssqrt((r1r2r3)/(K*dd*d*KK)); if Rms1e-3 stepMaxIter disp() disp([第,num2str(step),步EM算法收敛,均方根误差为,num2str(Rms)]) disp(各分量权重) disp(Alpha_est) disp(当前各类中心点) disp(Mu_est) disp(当前各类标准差) disp(Sigma_est) disp(----------------------------------) break; end disp(---------------------------------------) disp([第,num2str(step),次EM算法参数估计完成]) disp([均方根误差:,num2str(Rms)]) end end function WuniPosteriorGMM(X,Alpha,Mu,Sigma) %高斯混合模型隐变量的后验概率 %********************************************************* %输入 % X N*1 样本集(N个样本) % Alpha: K*1 权重 % Mu: K*1 均值 % Sigma K*1 标准差 %输出 % W: N*K 后验概率 %********************************************************* [N,d]size(X); Ksize(Alpha,1); if d1 YuniGMMpdf(X,Alpha,Mu,Sigma); for k1:K W(:,k)Alpha(k,:)*uniNormpdf(X,Mu(k,:),Sigma(k,:))./Y; end else disp(This is univariate GMM!); end end最后在一个主函数中来验证以上函数function mainForEMuniGMM() clear all clc %高斯混合模型分布的参数三个高斯分布组成 Mu [-10;0;8]; % 均值 Sigma [2;1;5]; % 标准差 Alpha[0.2;0.5;0.3]; % 权重 N1000; %生成样本集 XSamplingUniGMM(N,Alpha,Mu,Sigma); subplot(2,1,1) histogram(X,BinWidth, 0.5,BinLimits,[-20,20]); title(EM算法训练前样本直方图); xlabel(值); ylabel(频率); %用EM算法训练 K3; %高斯分布成分数 [Alpha_est,Mu_est,Sigma_est]uniEM_GMM(X,K); X_estSamplingUniGMM(N,Alpha_est,Mu_est,Sigma_est); subplot(2,1,2) histogram(X_est,BinWidth, 0.5,BinLimits,[-20,20]); title(EM算法训练后GMM样本直方图); xlabel(值); ylabel(频率); end运行的结果为再给一个有关图像分割的简单应用例子。GMM是图像分割领域经典的概率聚类方法核心是用多个高斯分布拟合图像的像素特征分布实现像素的软概率分类。用它进行图像分割的基本原理是假设图像中每个类别的像素特征服从独立的高斯分布通过‌EM算法‌迭代估计模型的均值、协方差和混合权重等参数再结合贝叶斯决策规则将像素分配给概率最高的类别完成分割。例2. 编写程序用单变量高斯混合模型实现图像分割。%使用高斯混合模型Gaussian Mixture Model, GMM进行图像分割 clear all clc I imread(yellowlily.jpg); % 读取图像 if size(I, 3) 3 I rgb2gray(I); % 转换为灰度图像 end subplot(1,2,1) imshow(I); title(Original Image); Xdouble(I(:)); % 用EM算法训练 K3; % 指定要使用的组件数即高斯分布的数量 [Alpha,Mu,Sigma]uniEM_GMM(X,K); % 计算每个像素属于每个高斯分布的概率(后验概率) WuniPosteriorGMM(X,Alpha,Mu,Sigma); % 选择一个阈值来分割图像例如使用最大后验概率准则选择每个像素的所属类别 [~,sImage] max(W, [], 2); sImage reshape(sImage, size(I)); subplot(1,2,2) imshow(label2rgb(sImage, jet, w, shuffle)); % 使用不同颜色显示不同的分割区域 title(Segmented Image);运行效果如下四、参考文献[1] Dempster, A.P., Laird, N.M. and Rubin, D.B., 1977. Maximum likelihood from incomplete data via the EM algorithm. Journal of the royal statistical society. Series B (methodological), pp.1-38.