高斯混合模型的matlab实现
原文地址:http://www.crescentmoon.info/?p=463
高斯混合函数实现部分是基本上是转载的的pluskid大神文章里的里的代码,加了一点注释,并根据他给的方法二解决 covariance 矩阵 singular 的问题。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
|
function varargout = gmm(X, K_or_centroids)
% ============================================================
%转载自http://blog.pluskid.org/?p=39
% Expectation-Maximization iteration implementation of
% Gaussian Mixture Model.
%
% PX = GMM(X, K_OR_CENTROIDS)
% [PX MODEL] = GMM(X, K_OR_CENTROIDS)
%
% - X: N-by-D data matrix.%需要注意的是这里的X包括了全部
% - K_OR_CENTROIDS: either K indicating the number of
% components or a K-by-D matrix indicating the
% choosing of the initial K centroids.
%
% - PX: N-by-K matrix indicating the probability of each
% component generating each point.
% - MODEL: a structure containing the parameters for a GMM:
% MODEL.Miu: a K-by-D matrix.
% MODEL.Sigma: a D-by-D-by-K matrix.
% MODEL.Pi: a 1-by-K vector.
% ============================================================
threshold = 1e-15;
[N, D] = size(X);
if isscalar(K_or_centroids)
K = K_or_centroids;
% randomly pick centroids
rndp = randperm(N);
centroids = X(rndp(1:K),:);
else
K = size(K_or_centroids, 1);
centroids = K_or_centroids;
end
% initial values
[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); %对矩阵的理解真是出神入化,
% 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}; %注意!!!!!这里和大神代码不同,他返回的是px,而我是 pGamma
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 对二维向量绘图
|
注意:
pluskid大神这里最后返回的是px,我觉得非常奇怪,因为PRML里对点做hard assignment时是根据后验概率来判别的。于是我在大神博客上问了一下,他的解释是最大似然和最大后验的区别,前者是挑x被各个模型产生的概率最大的那个,而后者加上了先验知识,各有道理。一句话就茅塞顿开,真大神也~
然后调用该函数对数据进行聚类,要是数据是二维的或三维的,顺便画个图。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
|
X=[...
1 1 1.1;
1 2 1;
2 1 1;
2 2 1;
5 5 2;
5 6 2;
6 5 2;
6 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
|
下面是对鸢尾花数据集聚类的结果。选了部分维度画图
二维图:
三维图:
高斯混合模型的matlab实现相关推荐
- matlab构建高斯混合模型,使用matlab创建高斯混合模型及绘图
Matlab提供了根据几个独立的高斯模型创建Gaussian Mixture Model(GMM)的函数,即fitgmdist.关于该模型的具体使用方法以及绘制生成的GMM的图形的方法,如下代码所示: ...
- 高斯混合模型(matlab代码+注释)
这里我学习的是Statistical Patte7rn Recognition Toolbox中的emgmm代码,代码中的主要知识点在之前的GMM文档中基本解释清楚,包括EM算法中的两个步骤.我自己先 ...
- K-means算法、高斯混合模型 matlab
K-means算法.高斯混合模型 简介: 本节介绍STANFORD机器学习公开课中的第12.13集视频中的算法:K-means算法.高斯混合模型(GMM).(9.10.11集不进行介绍,略过了哈) 一 ...
- matlab注释分析高斯混合模型
源码from:http://blog.csdn.net/abcjennifer/article/details/8198352#reply Rachel-Zhang 提供的源码. 高斯混合模型 没有输 ...
- 数字图像处理拓展题目——利用Matlab实现动态目标检测 二帧差法、ViBe法、高斯混合模型法,可应用于学生递东西行为检测
1.二帧差法实现动态目标检测 先上效果图: 利用GUI界面显示出来效果图为: 实现流程 1.利用matlab中的VideoReader函数读取视频流. 2.帧差法:获得视频帧数,用for循环对图像每相 ...
- GMM高斯混合模型聚类的EM估计过程matlab仿真
目录 1.算法概述 2.仿真效果 3.MATLAB源码 1.算法概述 高斯混合模型(Gaussian Mixed Model)指的是多个高斯分布函数的线性组合,理论上GMM可以拟合出任意类型的分布,通 ...
- matlab高斯拟合多峰,MATLAB用“fitgmdist”函数拟合高斯混合模型(一维数据)
MATLAB用"fitgmdist"函数拟合高斯混合模型(一维数据) 在MATLAB中"fitgmdist"的用法及其GMM聚类算法中介绍过"fitg ...
- 高斯混合模型聚类(GMM)matlab实现
Gaussian Mixture Model ,就是假设数据服从 Mixture Gaussian Distribution ,换句话说,数据可以看作是从数个 Gaussian Distributio ...
- GMM高斯混合模型及EM算法(matlab实现)
单元 %绘制男女生身高的GMM clc clear all %男女生共取2000人,女生平均身高163,男声平均身高180 male=180+sqrt(10)*randn(1,1000); %产生均值 ...
最新文章
- Linux下控制环境变量
- 模板 - 输入输出优化
- mysql数据库从删库到跑路之mysql多表查询
- 计算机辅助教学 林筑英,视频教学制作技巧.doc
- 网站如何集成百度UEditor编辑器
- Java算法面试题:编写一个程序,将e:\neck目录下的所有.java文件复制到e:\jpg目录下,并将原来文件的扩展名从.java改为.jpg...
- 9.3. where 优化
- 跨域请求——jsonp与cors
- Java数据类型转换超详解
- VBA函数 find
- cxf框架Demo1
- IntelliJ IDEA 12.0搭建Maven Web SSH2架构项目示例(二)
- adb基础命令学习随笔
- Flash动画设计交互式按钮
- 微信小程序父子组件方法调用方法汇总
- 设备管理 android问号,设备管理查有问号怎么修理
- 哈工大计算机854考研经验分享
- Mysql修改表中字段名称、字段类型
- 壁纸搜索系统/壁纸管理系统的设计与实现
- 【JAVA】500勇士问题,杀掉第三个人