高斯混合模型笔记matlab

人生在世，感谢这两行代码的救命之恩，哈哈哈哈

gmm=gmdistribution.fit(X,2);

P = posterior(gmm,X);

一：高斯分量数目，参数寻优，利用AIC或者BIC准则

close all;
clear;
clc;
%% 数据
s = rng(5,'v5normal');
mu = round((rand(3,2)-0.5)*19)+1;
sigma = round(rand(3,2)*40)/10+1;
X = [mvnrnd(mu(1,:),sigma(1,:),200); ...
     mvnrnd(mu(2,:),sigma(2,:),300); ...
     mvnrnd(mu(3,:),sigma(3,:),400)];
 %% AIC准则 进行参数寻优
% AIC = zeros(1,4);
% NlogL = AIC;
% GM = cell(1,4);
% for k = 1:4
%     GM{k} = gmdistribution.fit(X,k);
%     AIC(k)= GM{k}.AIC;
%     NlogL(k) = GM{k}.NlogL;
% end
%[minAIC,numComponents] = min(AIC);
%% BIC准则 进行参数寻优
BIC = zeros(1,5);
NlogL = BIC;
GM = cell(1,5);
for k = 1:5
    GM{k} = gmdistribution.fit(X,k);
    BIC(k)= GM{k}.BIC;
    NlogL(k) = GM{k}.NlogL;
end
[minBIC,numComponents] = min(BIC);

　　

　　

正经点，写下详细内容：

close all;
clear;
clc;
%% 使用高斯分布（正态分布）
% 随机生成3个中心以及标准差
s = rng(5,'v5normal');
mu = round((rand(3,2)-0.5)*19)+1;
sigma = round(rand(3,2)*40)/10+1;
X = [mvnrnd(mu(1,:),sigma(1,:),200); ...
     mvnrnd(mu(2,:),sigma(2,:),300); ...
     mvnrnd(mu(3,:),sigma(3,:),400)];
%% 产生二元正态分布随机数作为样本数据
% MU1 = [1 2];
% SIGMA1 = [2 0; 0 .5];
% MU2 = [-3 -5];
% SIGMA2 = [1 0; 0 1];
% X = [mvnrnd(MU1,SIGMA1,1000);mvnrnd(MU2,SIGMA2,1000)];
%% 使用高斯混合模型进行聚类分析（我的主要目的是得到每个高斯模型上数据的概率分布情况）
options = statset('Display','off');
gmm = gmdistribution.fit(X,3,'Options',options);
P = posterior(gmm,X);%输出所有样本在每个高斯分量上的概率
%以上是我最重要的两行代码！！！
idx = cluster(gmm,X);%聚类结果：每个样本属于哪个高斯分布，当然是按照概率最大的那个啦，心情好，说话都怪了，哈哈哈
%以下想把每个类单独筛选出来
cluster1 = X(idx == 1,:);%归到第一个高斯分类的样本集合
cluster2 = X(idx == 2,:);%归到第二个
cluster3 = X(idx == 3,:);%归到第二个
p1=P(idx==1,1);%把归到第一个高斯分类的样本在这个高斯分量上的概率取出来
p2=P(idx==2,2);%把归到第二个高斯分类的样本在这个高斯分量上的概率取出来
p3=P(idx==3,3);

二次改进中。。。

close all;
clear;
clc;
%% 使用高斯分布（正态分布）
% 随机生成3个中心以及标准差
s = rng(5,'v5normal');
mu = round((rand(3,2)-0.5)*19)+1;
sigma = round(rand(3,2)*40)/10+1;
X = [mvnrnd(mu(1,:),sigma(1,:),200); ...
     mvnrnd(mu(2,:),sigma(2,:),300); ...
     mvnrnd(mu(3,:),sigma(3,:),400)];
%% 产生二元正态分布随机数作为样本数据
% MU1 = [1 2];
% SIGMA1 = [2 0; 0 .5];
% MU2 = [-3 -5];
% SIGMA2 = [1 0; 0 1];
% X = [mvnrnd(MU1,SIGMA1,1000);mvnrnd(MU2,SIGMA2,1000)];
%% 使用高斯混合模型进行聚类分析（我的主要目的是得到每个高斯模型上数据的概率分布情况）
options = statset('Display','off');
gmm = gmdistribution.fit(X,3,'Options',options);
P = posterior(gmm,X);%输出所有样本在每个高斯分量上的概率
%以上是我最重要的两行代码！！！
idx = cluster(gmm,X);%聚类结果：每个样本属于哪个高斯分布，当然是按照概率最大的那个啦，心情好，说话都怪了，哈哈哈
%以下想把每个类单独筛选出来->为了节省内存，注释掉，直接赋值给gmmdata
% cluster1 = X(idx == 1,:);%归到第一个高斯分类的样本集合
% cluster2 = X(idx == 2,:);%归到第二个
% cluster3 = X(idx == 3,:);%归到第二个
% p1=P(idx==1,1);%把归到第一个高斯分类的样本在这个高斯分量上的概率取出来
% p2=P(idx==2,2);%把归到第二个高斯分类的样本在这个高斯分量上的概率取出来
% p3=P(idx==3,3);
gmmdata(1)=struct('cluster',X(idx == 1,:),'p',P(idx==1,1));%取值分别就是上面的cluster1和p1  
gmmdata(2)=struct('cluster',X(idx == 2,:),'p',P(idx==2,2));
gmmdata(3)=struct('cluster',X(idx == 3,:),'p',P(idx==3,3)); 
%% 对各个簇根据概率进行分区，再根据每个区内的样本数进行采样
n=length(gmmdata);
for i=1:n
    temp=gmmdata(i).cluster;%取出当前高斯分量里的样本
    tempprob=gmmdata(i).p;%取出当前高斯分量样本的概率
    [q1,~]=find(tempprob(:,1)>0&tempprob(:,1)<=0.4);%返回概率区间在0~0.4之间的样本在tempprob上的行号
    Giq(1).sample=temp(q1,:);%取出区间0~0.4之间的数，依然采用结构体进行依次存储
    Giq(1).number=size(q1,1);%存储该区间样本数,即q1的行数
    [q2,~]=find(tempprob(:,1)>0.4&tempprob(:,1)<=0.7);%返回概率区间在0.4~0.7之间的样本在tempprob上的行号
    Giq(2).sample=temp(q2,:);%取出区间0.4~0.7之间的数
    Giq(2).number=size(q2,1);%存储该区间样本数
    [q3,~]=find(tempprob(:,1)>0.7&tempprob(:,1)<=1);%返回概率区间在0.7~1之间的样本在tempprob上的行号
    Giq(3).sample=temp(q3,:);%取出区间0.7~1之间的数
    Giq(3).number=size(q3,1);%存储该区间样本数
    for j=1:length(Giq)
        rus=randperm(Giq(j).number,fix(Giq(j).number/3));%这里采样数为三分之一
        underSample=Giq(j).sample(rus,:);%取出第j个区间里随机采样的样本
        data_new_temp= [data_new_temp;underSample];%把每个区间的采样样本合并，此循环结束后，data_new_temp即为从第i个簇采样的样本
    end
    data_new=[data_new;data_new_temp];
end
%% 采样得到负类样本的输出

三次改进。。。写成函数，然后多利用结构体数组，以及内存预分配

function [ data_new ] = undSampByGmm(X,K )
%UNTITLED 此处显示有关此函数的摘要
%   此处显示详细说明
%% %记录下总样本数和数据维数
[R,C]=size(X);
%% 使用高斯混合模型进行聚类分析（我的主要目的是得到每个高斯模型上数据的概率分布情况）
options = statset('Display','off');
gmm = gmdistribution.fit(X,K,'Options',options);%得到高斯混合模型
P = posterior(gmm,X);%输出所有样本在每个高斯分量上的概率
%以上是我最重要的两行代码！！！
idx = cluster(gmm,X);%聚类结果：每个样本属于哪个高斯分布，当然是按照概率最大的那个啦，心情好，说话都怪了，哈哈哈
gmmdata=repmat(struct('cluster',[R,C],'p',[R,1]),1,3);%结构体预先分配内存，减少迭代中不断分配内存带来的消耗
for i=1:K
    gmmdata(i)=struct('cluster',X(idx == i,:),'p',P(idx==i,i));%取值分别就是上面的cluster1和p1  
end
%% 对各个簇根据概率进行分区，再根据每个区内的样本数进行采样
n=length(gmmdata);%高斯分量个数
%c=max([size(gmmdata(1).p,1),size(gmmdata(2).p,1),size(gmmdata(3).p,1)]);%规模最大簇的样本数
%data_new_temp(c,C)=0;%预分配内存
%data_new(fix(R/3),C)=0;%预分配内存
data_new=[];
for i=1:n
    temp=gmmdata(i).cluster;%取出当前高斯分量里的样本
    tempprob=gmmdata(i).p;%取出当前高斯分量样本的概率
    [q1,~]=find(tempprob(:,1)>0&tempprob(:,1)<=0.4);%返回概率区间在0~0.4之间的样本在tempprob上的行号
    Giq(1).sample=temp(q1,:);%取出区间0~0.4之间的数，依然采用结构体进行依次存储
    Giq(1).number=size(q1,1);%存储该区间样本数,即q1的行数
    [q2,~]=find(tempprob(:,1)>0.4&tempprob(:,1)<=0.7);%返回概率区间在0.4~0.7之间的样本在tempprob上的行号
    Giq(2).sample=temp(q2,:);%取出区间0.4~0.7之间的数
    Giq(2).number=size(q2,1);%存储该区间样本数
    [q3,~]=find(tempprob(:,1)>0.7&tempprob(:,1)<=1);%返回概率区间在0.7~1之间的样本在tempprob上的行号
    Giq(3).sample=temp(q3,:);%取出区间0.7~1之间的数
    Giq(3).number=size(q3,1);%存储该区间样本数
    data_new_temp=[];
    for j=1:length(Giq)
        rus=randperm(Giq(j).number,fix(Giq(j).number/3));%这里采样数为三分之一
        underSample=Giq(j).sample(rus,:);%取出第j个区间里随机采样的样本
        data_new_temp= [data_new_temp;underSample];%把每个区间的采样样本合并，此循环结束后，data_new_temp即为从第i个簇采样的样本
    end
    data_new=[data_new;data_new_temp];
end
end

调用：

close all;
clear;
clc;
%% 使用高斯分布（正态分布）
% 随机生成3个中心以及标准差
s = rng(5,'v5normal');
mu = round((rand(3,2)-0.5)*19)+1;
sigma = round(rand(3,2)*40)/10+1;
X = [mvnrnd(mu(1,:),sigma(1,:),200); ...
     mvnrnd(mu(2,:),sigma(2,:),300); ...
     mvnrnd(mu(3,:),sigma(3,:),400)];
[ data_new ] = undSampByGmm(X,3 );

拿自己的数据测试就出错了。。。

错误内容：

 错误使用 gmcluster (line 193) Ill-conditioned covariance created at iteration 19.

上网找资料：解决方法参考https://cn.mathworks.com/help/stats/gmdistribution.fit.html   仔细阅读In general, you can avoid getting ill-conditioned covariance matrices by using one of the following precautions:下面的内容

In some cases, gmdistribution may converge to a solution where one or more of the components has an ill-conditioned or singular covariance matrix.

The following issues may result in an ill-conditioned covariance matrix:

The number of dimension of your data is relatively high and there are not enough observations.
Some of the features (variables) of your data are highly correlated.
Some or all the features are discrete.
You tried to fit the data to too many components.
In general, you can avoid getting ill-conditioned covariance matrices by using one of the following precautions:

Pre-process your data to remove correlated features.
Set 'SharedCov' to true to use an equal covariance matrix for every component.
Set 'CovType' to 'diagonal'.
Use 'Regularize' to add a very small positive number to the diagonal of every covariance matrix.
Try another set of initial values.

网上中文翻译：

出现以下情况可能导致协方差矩阵，ill-Condition或者变为奇异阵（行列式为0，无法求逆矩阵，而计算多维高斯密度需要用到逆矩阵）：

1.相对于已有的数据量来说，数据维度相对太高，

2.数据属性之间相关度过高，

3.某些或全部属性为离散，

4.或者试图分成太多的聚类。

解决方案大致为：

1.预处理数据，移除相关度太高的属性（这个一般可以使用mutual information，卡方检验，Pearsion product cofficient来处理）

2.让每个聚类分享同一个协方差矩阵

3.使协方差矩阵变成对角阵

4.给协方差矩阵的对角线加一个很小的正数值（这个在我的程序中进行了试验，发现他的意义在于在某次迭代过后协方差矩阵变成全部元素为0的矩阵，加上一个很小的正数后是过程可以继续）

我尝试了第四种，改为gmdistribution.fit(X,3, 'Regularize', 1e-5)之后便得以解决

————————————————————————————————————

采样前后对比

close all;
clear;
clc;
%% 数据
%采样前数据
% 使用高斯分布（正态分布）
% 随机生成3个中心以及标准差
MU1 = [1 2];
SIGMA1 = [2 0; 0 .5];
MU2 = [-4 -5];
SIGMA2 = [1 0; 0 1];
MU3 = [4 -3];
SIGMA3 = [2 0; 0 .5];
X = [mvnrnd(MU1,SIGMA1,1000);mvnrnd(MU2,SIGMA2,1000);mvnrnd(MU3,SIGMA3,1000)];
%采样后数据
newData= undSampByGmm(X,3);
%% 算下均值和方差
u1=mean(X,1);%我需要的列均值
u2=mean(newData,1);%我需要的列均值
s1=power(std(X,0,1),2);
s2=power(std(newData,0,1),2);
%% 高斯混合模型
GMModel = gmdistribution.fit(X,3);
GMModel2 =gmdistribution.fit(newData,3);
%% 画图idx = cluster(obj,X);
idx = cluster(GMModel,X);
cluster1 = X(idx == 1,:);
cluster2 = X(idx == 2,:);
cluster3 = X(idx == 3,:);
%采样后
idx2 = cluster(GMModel2,newData);
cluster4 = newData(idx2 == 1,:);
cluster5= newData(idx2 == 2,:);
cluster6 = X(idx2 == 3,:);

%h1 = scatter(cluster1(:,1),cluster1(:,2),10,'r.');
%h2 = scatter(cluster2(:,1),cluster2(:,2),10,'g.');
%legend([h1 h2],'Cluster 1','Cluster 2','Location','NW')
subplot(1,2,1);%划分窗口
h1 = plot(cluster1(:,1),cluster1(:,2),'r.');
hold on
h2 = plot(cluster2(:,1),cluster2(:,2),'c.');
hold on
h3 = plot(cluster3(:,1),cluster3(:,2),'b.');
legend([h1 h2 h3],'Cluster 1','Cluster 2','Cluster 3')
hold on
ezcontour(@(x1,x2)pdf(GMModel,[x1 x2]),[-10 10],[-8 6]); %用这个！！！
title('{\bf 采样前数据分布}');   %用这个！！！
grid on%画网格，去掉，用grid off
subplot(1,2,2);%划分窗口
h4 = plot(cluster4(:,1),cluster4(:,2),'r.');
hold on
h5 = plot(cluster5(:,1),cluster5(:,2),'c.');
hold on
h6 = plot(cluster6(:,1),cluster6(:,2),'b.');
legend([h4 h5 h6],'Cluster 1','Cluster 2','Cluster 3')
hold on
ezcontour(@(x1,x2)pdf(GMModel,[x1 x2]),[-10 10],[-8 6]); %用这个！！！
title('{\bf 采样后数据分布}');   %用这个！！！
grid on
hold


%% 多维情况
% MU1 = [1 2 1];
% SIGMA1 = [2 0 0;0 1 0;0 0 .5];
% MU2 = [-3 -5 -4];
% SIGMA2 = [1 0 0;0 1 0;0 0 1];
% X = [mvnrnd(MU1,SIGMA1,1000);mvnrnd(MU2,SIGMA2,1000)];
% obj = gmdistribution.fit(X,2);
% idx = cluster(obj,X);
% cluster1 = X(idx == 1,:);
% cluster2 = X(idx == 2,:);
% h1 = plot3(cluster1(:,1),cluster1(:,2),cluster1(:,3),'r.');
% hold on
% h2 = plot3(cluster2(:,1),cluster2(:,2),cluster2(:,3),'g.');
% legend([h1 h2],'Cluster 1','Cluster 2')
% grid on