机器学习系列：高斯混合模型(1)

上一篇文章描述了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 Y=uniGMMpdf(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); N=size(X,1); Y=zeros(N,K); if d<=1 for k=1:K Y(:,k)=uniNormpdf(X,Mu(k,:),Sigma(k,:)); end Y=Y*Alpha; else disp('This is univariate GMM!'); end end

再用拒绝采样算法为高斯混合模型编写一个采样函数（这里用区间[a,b]上的均匀分布作为辅助分布）：

function samples=SamplingUniGMM(numberSample,Alpha,Mu,Sigma) %高斯混合模型的样本生成函数（这里用拒绝采样方法实现） %这里选择一个均匀分布作为辅助分布。 %********************************************************* %输入： % numberSample: 样本数量 % Alpha: 权重(K*1) % Mu: 均值(K*1) % Sigma： 标准差(K*1) %输出： % samples: 样本集 %********************************************************* %确定分布的范围[a,b] b=max(Mu)+(max(Sigma))^2/2; % 均匀分布的下限 a=min(Mu)-(max(Sigma))^2/2; % 均匀分布的上限 % 样本集 samples = zeros(numberSample,1); %存放样本 %一个合适的M M=(b-a)/(sqrt(2*pi)*min(Sigma)); accepted=0; for i=1:numberSample while 1 %从辅助分布采样 x=a+(b-a)*rand(); % 计算接受概率 gamma=uniGMMpdf(x,Alpha,Mu,Sigma)/(M*uniformpdf(x,a,b)); %生成一个[0,1]上均匀分布随机数 u=rand(); if u<=gamma %接受这个样本 accepted = accepted + 1; samples(accepted,1)=x; break; % 接受后退出循环，继续下一个样本的生成 end end end end function pdf=uniformpdf(x,a,b) %区间[a,b]上的均匀分布概率密度函数 %********************************************************* % 输入: % x 1*1 样本 % a 1*1 左端点 % b 1*1 右端点 % 输出: % pdf 1*1 概率密度值 %********************************************************** if x>=a & x<=b pdf=1/(b-a); else pdf=0; 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]); %依据各维度的最大最小值构建参数的初始值 tMin=min(X); tMax=max(X); Mu_est=repmat(tMin,K,1)+linspace(0,1,K)'*(tMax-tMin); MaxIter=1e+5; for step=1:MaxIter Mu=Mu_est; Sigma=Sigma_est; Alpha=Alpha_est; %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %E-Step: % %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %计算后验概率 W=uniPosteriorGMM(X,Alpha_est,Mu_est,Sigma_est); %%%%%%%%%%%%%%%%%%%%% %M-Step: % %%%%%%%%%%%%%%%%%%%%% % Alpha Sw=sum(W,1); Alpha_est=Sw'/N; % Mu Mu_est=(X'*W./Sw)'; % Sigma Y=X-repmat(Mu_est',size(X,1),1); Y=Y.*Y; Sigma_est=(sum(W.*Y,1)./Sw)'; Sigma_est=sqrt(Sigma_est); % 收敛条件(计算均方根误差) r1=sum((Mu_est-Mu).^2,'all'); r2=sum((Sigma_est-Sigma).^2,'all'); r3=sum((Alpha_est-Alpha).^2,'all'); Rms=sqrt((r1+r2+r3)/(K*d+d*d*K+K)); if Rms<1e-3 && step<MaxIter 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 W=uniPosteriorGMM(X,Alpha,Mu,Sigma) %高斯混合模型隐变量的后验概率 %********************************************************* %输入： % X N*1 样本集(N个样本) % Alpha: K*1 权重 % Mu: K*1 均值 % Sigma： K*1 标准差 %输出： % W: N*K 后验概率 %********************************************************* [N,d]=size(X); K=size(Alpha,1); if d<=1 Y=uniGMMpdf(X,Alpha,Mu,Sigma); for k=1: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]; % 权重 N=1000; %生成样本集 X=SamplingUniGMM(N,Alpha,Mu,Sigma); subplot(2,1,1) histogram(X,'BinWidth', 0.5,'BinLimits',[-20,20]); title('EM算法训练前样本直方图'); xlabel('值'); ylabel('频率'); %用EM算法训练 K=3; %高斯分布成分数 [Alpha_est,Mu_est,Sigma_est]=uniEM_GMM(X,K); X_est=SamplingUniGMM(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'); X=double(I(:)); % 用EM算法训练 K=3; % 指定要使用的组件数，即高斯分布的数量 [Alpha,Mu,Sigma]=uniEM_GMM(X,K); % 计算每个像素属于每个高斯分布的概率(后验概率) W=uniPosteriorGMM(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.