MATLAB中“fitgmdist”的用法及其GMM聚类算法
作者：凯鲁嘎吉 - 博客园 http://www.cnblogs.com/kailugaji/
高斯混合模型的基本原理：聚类——GMM，MATLAB官方文档中有关于fitgmdist的介绍：fitgmdist。我之前写过有关GMM聚类的算法：GMM算法的matlab程序。这篇文章主要应用MATLAB自带的函数来进行聚类。
1. fitgmdist函数介绍
fitgmdist的使用形式：gmm = fitgmdist(X,k,Name,Value)
输入
‘RegularizationValue’, 0。（取值：0， 0.1， 0.01,....，正则化系数，防止协方差奇异）
'CovarianceType', 'full'。（取值： 'full'，协方差矩阵是非对角阵，'diagonal'，协方差矩阵为对角阵）
‘Start’, 'plus'。 （取值：‘randSample’，随机初始化，‘plus’，k-means++初始化，‘S’，自定义初始化），其中S = struct('mu',init_Mu,'Sigma',init_Sigma,'ComponentProportion',init_Components);
‘Options’，statset('Display', 'final', 'MaxIter', MaxIter, 'TolFun', TolFun)。 （'Display'有三个取值：‘final’ 显示最终的输出结果、‘iter’ 显示每次迭代的结果、‘off’ 不显示优化参数信息；'MaxIter'：默认100，最大迭代次数；'TolFun'：默认1e-6，目标函数的终止误差）
输出
gmm.mu：更新完后的聚类中心（均值）
gmm.Sigma：更新完后的协方差矩阵
gmm.ComponentProportion：更新完后的混合比例
gmm.NegativeLogLikelihood：更新完后的负对数似然函数
gmm.NumIterations：实际迭代次数
gmm.BIC：贝叶斯信息准则，用于模型选择
更多参数，请在命令行输入properties(gmm)
2. 高斯混合模型聚类实例
generate.m
- function data=generate()
- %生成数据
- mu1 = [1 2];
- Sigma1 = [2 0; 0 0.5];
- mu2 = [-1 -2];
- Sigma2 = [1 0;0 1];
- data = [mvnrnd(mu1,Sigma1,400), ones(400,1);mvnrnd(mu2,Sigma2,600), 2*ones(600,1)];
- X=[data(:, 1), data(:, 2)];
- figure(1)
- plot(X(:,1), X(:,2),'bo')
- title('Scatter Plot')
- xlim([min(X(:)) max(X(:))]) % Make axes have the same scale
- ylim([min(X(:)) max(X(:))])
具体数据
GMM_main.m
- function [accuracy,NumIterations]=GMM_main(data, K)
- %主函数
- [~, data_dim]=size(data);
- X=data(:, 1:data_dim-1); %数据
- real_label=data(:, data_dim);
- [label, ~, NumIterations]=Matlab_gmm_2(X, K);
- accuracy=succeed(real_label,K,label);
Matlab_gmm.m
- function [label, NegativeLogLikelihood, NumIterations]=Matlab_gmm(X, K)
- %协方差矩阵为对角阵，数据独立同分布
- [X_num,X_dim]=size(X);
- para_sigma_inv=zeros(X_dim, X_dim, K);
- N_pdf=zeros(X_num, K); %单高斯分布的概率密度函数
- RegularizationValue=0.001; %正则化系数，协方差矩阵求逆
- MaxIter=100; %最大迭代次数
- TolFun=1e-8; %终止条件
- % 自己设置初始化参数
- % init_Mu = [1 1; 2 2];
- % init_Sigma(:,:,1) = [1 1; 1 2];
- % init_Sigma(:,:,2) = 2*[1 1; 1 2];
- % init_Components = [1/2,1/2];
- % S = struct('mu',init_Mu,'Sigma',init_Sigma,'ComponentProportion',init_Components);
- % gmm=fitgmdist(X, K, 'RegularizationValue', RegularizationValue, 'CovarianceType', 'diagonal', 'Start', 'S', 'Options', statset('Display', 'final', 'MaxIter', MaxIter, 'TolFun', TolFun));
- gmm=fitgmdist(X, K, 'RegularizationValue', RegularizationValue, 'CovarianceType', 'diagonal', 'Start', 'plus', 'Options', statset('Display', 'final', 'MaxIter', MaxIter, 'TolFun', TolFun));
- NegativeLogLikelihood=gmm.NegativeLogLikelihood;
- NumIterations=gmm.NumIterations; %迭代次数
- mu=gmm.mu; %均值
- Sigma=gmm.Sigma; %协方差矩阵
- ComponentProportion=gmm.ComponentProportion; %混合比例
- for k=1:K
- sigma_inv=1./Sigma(:,:,k); %sigma的逆矩阵,(X_dim, X_dim)的矩阵
- para_sigma_inv(:, :, k)=diag(sigma_inv); %sigma^(-1)
- end
- for k=1:K
- coefficient=(2*pi)^(-X_dim/2)*sqrt(det(para_sigma_inv(:, :, k))); %高斯分布的概率密度函数e左边的系数
- X_miu=X-repmat(mu(k,:), X_num, 1); %X-miu: (X_num, X_dim)的矩阵
- exp_up=sum((X_miu*para_sigma_inv(:, :, k)).*X_miu,2); %指数的幂，(X-miu)'*sigma^(-1)*(X-miu)
- N_pdf(:,k)=coefficient*exp(-0.5*exp_up);
- end
- responsivity=N_pdf.*repmat(ComponentProportion,X_num,1); %响应度responsivity的分子，（X_num,K）的矩阵
- responsivity=responsivity./repmat(sum(responsivity,2),1,K); %responsivity:在当前模型下第n个观测数据来自第k个分模型的概率，即分模型k对观测数据Xn的响应度
- %聚类
- [~,label]=max(responsivity,[],2);
- figure(2)
- scatter(X(:,1),X(:,2),10,'.') % Scatter plot with points of size 10
- hold on
- gmPDF = @(x,y)reshape(pdf(gmm,[x(:) y(:)]),size(x));
- fcontour(gmPDF,[-6 6])
Matlab_gmm_2.m
- function [label, NegativeLogLikelihood, NumIterations]=Matlab_gmm_2(X, K)
- %协方差矩阵为非对角阵，数据不独立
- [X_num,X_dim]=size(X);
- N_pdf=zeros(X_num, K); %单高斯分布的概率密度函数
- RegularizationValue=0.001; %正则化系数，协方差矩阵求逆
- MaxIter=100; %最大迭代次数
- TolFun=1e-8; %终止条件
- % 自己设置初始化参数
- % init_Mu = [1 1; 2 2];
- % init_Sigma(:,:,1) = [1 1; 1 2];
- % init_Sigma(:,:,2) = 2*[1 1; 1 2];
- % init_Components = [1/2,1/2];
- % S = struct('mu',init_Mu,'Sigma',init_Sigma,'ComponentProportion',init_Components);
- % gmm=fitgmdist(X, K, 'RegularizationValue', RegularizationValue, 'CovarianceType', 'diagonal', 'Start', 'S', 'Options', statset('Display', 'final', 'MaxIter', MaxIter, 'TolFun', TolFun));
- gmm=fitgmdist(X, K, 'RegularizationValue', RegularizationValue, 'CovarianceType', 'full', 'Start', 'plus', 'Options', statset('Display', 'final', 'MaxIter', MaxIter, 'TolFun', TolFun));
- NegativeLogLikelihood=gmm.NegativeLogLikelihood;
- NumIterations=gmm.NumIterations; %迭代次数
- mu=gmm.mu; %均值
- Sigma=gmm.Sigma; %协方差矩阵
- ComponentProportion=gmm.ComponentProportion; %混合比例
- for k=1:K
- X_miu=X-repmat(mu(k,:), X_num, 1); %X-miu: (X_num, X_dim)的矩阵
- sigma_inv=inv(Sigma(:,:,k)); %sigma的逆矩阵,(X_dim, X_dim)的矩阵
- exp_up=sum((X_miu*sigma_inv).*X_miu,2); %指数的幂，(X-miu)'*sigma^(-1)*(X-miu)
- coefficient=(2*pi)^(-X_dim/2)*sqrt(det(sigma_inv)); %高斯分布的概率密度函数e左边的系数
- N_pdf(:,k)=coefficient*exp(-0.5*exp_up);
- end
- responsivity=N_pdf.*repmat(ComponentProportion,X_num,1); %响应度responsivity的分子，（X_num,K）的矩阵
- responsivity=responsivity./repmat(sum(responsivity,2),1,K); %responsivity:在当前模型下第n个观测数据来自第k个分模型的概率，即分模型k对观测数据Xn的响应度
- %聚类
- [~,label]=max(responsivity,[],2);
- figure(2)
- scatter(X(:,1),X(:,2),10,'.') % Scatter plot with points of size 10
- hold on
- gmPDF = @(x,y)reshape(pdf(gmm,[x(:) y(:)]),size(x));
- fcontour(gmPDF,[-6 6])
succeed.m
- function accuracy=succeed(real_label,K,id)
- %输入K：聚的类，id：训练后的聚类结果，N*1的矩阵
- N=size(id,1); %样本个数
- p=perms(1:K); %全排列矩阵
- p_col=size(p,1); %全排列的行数
- new_label=zeros(N,p_col); %聚类结果的所有可能取值，N*p_col
- num=zeros(1,p_col); %与真实聚类结果一样的个数
- %将训练结果全排列为N*p_col的矩阵，每一列为一种可能性
- for i=1:N
- for j=1:p_col
- for k=1:K
- if id(i)==k
- new_label(i,j)=p(j,k); %iris数据库，1 2 3
- end
- end
- end
- end
- %与真实结果比对，计算精确度
- for j=1:p_col
- for i=1:N
- if new_label(i,j)==real_label(i)
- num(j)=num(j)+1;
- end
- end
- end
- accuracy=max(num)/N;
结果
以第二种情况为例，数据不独立，协方差矩阵不是只在对角线上有元素。
- >> [accuracy,NumIterations]=GMM_main(data, 2)
- 32 iterations, log-likelihood = -3449.42
-
- accuracy =
-
- 0.995000000000000
-
-
- NumIterations =
-
- 32
经验首页