高斯混合模型
1、定义：
高斯混合模型就是用高斯概率密度函数（正态分布曲线）精确地量化事物，它是一个将事物分解为若干的基于高斯概率密度函数（正态分布曲线）形成的模型。
2、数学推导：
主要更新参数：均值μ 方差∑
设有一组数据x，其中隐含一随机过程z，x与z之间有join概率，
且满足联合分布：
Z服从多项式分布:
在Z给定的情况下，X服从正态分布：
这样的模型称为高斯混合模型
我们用该模型似然函数来描述已知分类结果，来求未知参数μ 和∑，似然函数如下：
如果z已知，可利用偏导来求参数：
但在实际中，高斯混合模型是EM算法的一个特例，迭代算法分为E-Step（估计z的概率）和M-Step（似然函数最大化）定义
3、代码编写：
X=[1 1 4; 1 2 1; 1 5 2; 2 3 1; 5 0 2; 5 6 2; 6 5 2; 4 6 2.1]%数据X自己改 K=2 %以上设置数据点 [Px,model]=gmm(X,K);%这里得到结果 [~,belong]=max(Px,[],2)%选概率最大的那个数，输出聚类结果 %以下绘图 z=zeros(size(X,1),3); if isscalar(K)%获得类别个数 K_number=K; else K_number = size(K, 1); end for k=1:K_number color=[rand rand rand];%建立颜色矩阵，随机给个颜色 csize=size(z(belong==k,:),1);%数一数有几行 z(belong==k,:)=repmat(color,csize,1);%对属于某一类下的点染色 end if size(X,2)==2 %二维 figure('color','w');%把背景改成白的 scatter(X(:,1),X(:,2),30,z) %axis off;%关掉坐标系显示 else if size(X,2)==3 %3维的情况 figure('color','w');%把背景改成白的 scatter3(X(:,1),X(:,2),X(:,3),30,z,'filled') end end
function varargout = gmm( X,K_or_centroids ) X=[1 1 0; 1 2 1; 2 1 1; 2 2 1; 5 5 2; 5 6 2; 6 5 2; 6 6 2.1];%数据X自己改 K_or_centroids=2 ; nargout=0; threshold = 1e-15; [N,D] = size(X);%行和列 if isscalar(K_or_centroids)%标量返回TRUE否则返回false K=K_or_centroids; rndp = randperm(N);%1到n这些数随机打乱得到的一个数字序列 centroids = X(rndp(1:K),:); else K = size(K_or_centroids,1);%行数 centroids = K_or_centroids; end %初始化 [pMiu pPi pSigma] = init_params(); Lprev = -inf; while true Px = calc_prob();%计算N(x|mu,sigma) % new value for pGamma pGamma = Px .* repmat(pPi, N, 1);%估计 gamma 是个N*K的矩阵 pGamma = pGamma ./ repmat(sum(pGamma, 2), 1, K);%%求每个样本由第K个聚类，也叫“component“生成的概率 % new value for parameters of each Component Nk = sum(pGamma, 1);%N_K pMiu = diag(1./Nk) * pGamma' * X; % 数字 *( K-by-N * N-by-D)加个括号有助理解 pPi = Nk/N; for kk = 1:K Xshift = X-repmat(pMiu(kk, : ), N, 1);%x-u pSigma(:, :, kk) = (Xshift' * ... (diag(pGamma(:, kk)) * Xshift)) / Nk(kk);%更新sigma end % check for convergence L = sum(log(Px*pPi')); if L-Lprev < threshold break; end Lprev = L; end if nargout == 1 varargout = {Px}; else model = []; model.Miu = pMiu; model.Sigma = pSigma; model.Pi = pPi; varargout = {pGamma, model}; end function [pMiu pPi pSigma] = init_params()%初始化参数 pMiu = centroids;% K-by-D matrix pPi = zeros(1, K);%1-by-K matrix pSigma = zeros(D, D, K);% % hard assign x to each centroids distmat = repmat(sum(X.*X, 2), 1, K) + ... % X is a N-by-D data matrix. repmat(sum(pMiu.*pMiu, 2)', N, 1) - ...% X->K列 U->N行 XU^T is N-by-K 2*X*pMiu';%计算每个点到K个中心的距离 [~, labels] = min(distmat, [], 2);%找到离X最近的pMiu，[C,I] labels代表这个最小值是从那列选出来的 for k=1:K Xk = X(labels == k, : );% Xk是所有被归到K类的X向量构成的矩阵 pPi(k) = size(Xk, 1)/N;% 数一数几个归到K类的 pSigma(:, :, k) = cov(Xk); %???计算协方差矩阵，D-by-D matrix,最小方差无偏估计 end end function Px = calc_prob() Px = zeros(N, K); for k = 1:K Xshift = X-repmat(pMiu(k, : ), N, 1);%x-u lemda=1e-5; conv=pSigma(:, :, k)+lemda*diag(diag(ones(D)));%这里处理singular问题，为协方差矩阵加上一个很小lemda*I inv_pSigma = inv(conv);%协方差的逆 tmp = sum((Xshift*inv_pSigma) .* Xshift, 2);%(X-U_k)sigma.*(X-U_k),tmp是个N*1的向量 coef = (2*pi)^(-D/2) * sqrt(det(inv_pSigma));%前面的参数 Px(:, k) = coef * exp(-0.5*tmp);%把数据点 x 带入到 Gaussian model 里得到的值 end end end %repmat 通过拓展向量到矩阵 %inv 求逆 %min 求矩阵最小值，可以返回标签 %X(labels == k, : ) 对行做筛选 % size(Xk, 1) 求矩阵的长或宽 %scatter 对二维向量绘图
4、输出结果：
X = 1.0000 1.0000 4.0000 1.0000 2.0000 1.0000 1.0000 5.0000 2.0000 2.0000 3.0000 1.0000 5.0000 0 2.0000 5.0000 6.0000 2.0000 6.0000 5.0000 2.0000 4.0000 6.0000 2.1000 K =
2 belong = 2 2 2 2 1 1 1 1
5、算法优缺点：
主要优点：投影后样本点不是得到一个确定的分类标记，而是得到每个类的概率
主要缺点：计算复杂度高；当前景与背景颜色相似时，会将前景错误的处理为背景。
6、高斯混合模型与K-means异同点
相同点： （1）需要指定聚类数目K （2）需要指定初始值，K-means的中心点，GMM的各个参数 （3）都是含有EM算法思想不同点： （1）优化目标函数不同，K-means:最短距离，GMM：最大化log似然估计 （2）E步的指标不同，K-means:点到中心的距离（硬指标），GMM：求解每个观测数据的每个component的概率（软指标）
