学习记录
聚类
聚类任务
背景
在无监督学习(密度估计、异常检测等)中,训练样本的标记信息是未知的(即不人为指定),旨在发现数据之间的内在联系和规律,为进一步的数据分析提供基础。
此类学习任务中研究最多、应用最广的是聚类。
定义
聚类将数据集中的样本划分为若干个通常是不相交(有的时候会相交)的子集(簇,也叫类),即将样本分类,但是分类后得到的簇的标记是未指定的,这里需要使用者自行定义。
作用
聚类既能作为一个单独过程,用于寻找数据内在的分布结构,可以作为分类等其他学习任务的前驱过程。
主要概括起来,就是寻找数据内在的规律和联系,进行整合,分类,得到的一些数据集(类)。
性能度量(有效性指标)
类似于监督性学习的性能度量。
聚类理想情况就是,簇内相似度高,簇间相似度低。
聚类的性能度量大致有两类,分为外部指标(具有参考模型)和内部指标(直接考察聚类结果)。
外部指标举例的话,比如给你一群人的特征,一种理想情况就是聚类为男|女。
外部指标
书上根据上面的数据给出了如下四个式子
其中 第一个S 表示C中属于相同簇, D 表示C中不属于相同簇, 第二个S 表示
C∗
中属于相同簇,D表示
C∗
中不属于相同簇。
很明显,只有这四种情况,所有得出 a + b + c + d = m (m - 1) / 2,其中m (m - 1) / 2 为所有数据两两相对的情况数。
基于上面的四个式子可以导出下面这些常用的聚类性能度量外部指标:
- Jaccard系数(JC)
JC=aa+b+c - FM指数(FMI)
FMI=aa+b⋅aa+c−−−−−−−−−−√ - Rand指数(RI)
RI=2(a+d)m(m−1)
显然,上述的三个性能度量的结果值均在[0, 1]区间,值越大越好。
内部指标
接着,又能得出下面这四个式子(基于样本之间的距离dist来得出: avg簇内平均,diam簇内最远,min簇间最近,cen簇间中心点间的距离)
基于上面四个式子,我们可以导出下面这些常用的聚类性能度量的内部指标
- DB指数(DBI)
DBI=1k∑i=1kmaxj≠i(avg(Ci)+avg(Cj)dcen(ui,uj))
分析可知,上式求得的就是簇之间两两比较,计算出二个簇的平均距离之和除以二个簇之间的中心距离,取最大值,并求和,最后得出的结果代表着簇间相似度,为什么呢?
这里就要说一下,dist,是一个函数,能够计算二个样本之间的距离,那么,这个得出的距离是什么呢?就是二者的相似性,距离越大,代表相似性越低,差异性越高,那么,cen就代表着二个簇之间的相似性,越小越相似,即上述式子就越大。
但是为什么是簇间相似度而不是簇内相似度,他计算的是cen簇内中间距离,得出的是二个簇之间的相似度。
- Dunn指数(DI)
DI=min1⩽i⩽k{minj≠i(dmin(Ci,Cj)max1⩽l⩽kdiam(Cl))}
类似于DBI的分析,也可以得出DI代表的是簇内相似度,比较的是diam簇内最远。
综上所述,DBI越小越号,DI越大越好(不懂分子那里为什么)。
距离计算
这里主要说的就是上面说的dist函数,若他是一个距离度量,则需满足一些基本性质
- 非负性(即算出的差异性不能为负数)
dist(xi,xj)⩽0 - 同一性(即自己跟自己距离要为0)
dist(xi,xj)=0当且仅当xi=xj - 对称性(即二点之间距离从这走和从对面走是一样的)
dist(xi,xj)=dist(xj,xi) - 直递性(三角不等式)
dist(xi,xj)⩽dist(xi,xk)+dist(xk,xj)
闵可夫斯基距离(MinkovDM)
预备知识:范数
上面这个式子就是范数 ||xi−xj||p
当p=2时,计算的就是欧式距离。
当p=1时,计算的就是曼哈顿距离(街区距离)。
属性分为连续属性(数值属性)和离散属性(列名属性),连续属性一般为有序属性,离散属性一般为无序属性,闵可夫斯基距离可用于有序属性的计算, 如果需要对无序属性进行计算,则需要使用VDM。
VDM(不同值的标量)
看到这里,感觉就是在计算二个值(a,b)在不同簇被引用的比例。
接着,我们就可以使用以下这些小连招- -,处理混合属性的时候
有的时候,距离度量可能不满足直递性,这样的距离成为非度量距离。
原型聚类(基于原型的聚类)
此类算法假设聚类结构能通过一组原型刻画,在现实聚类任务中极为常用。
原理
先对原型进行初始化,然后对原型进行迭代更新求解,采用不同的原型表示、不同的求解方式,将产生不同的算法。
k均值算法(kmeans)
算法流程
kmeans(X,k)
idx = kmeans(X,k)
idx = kmeans(X,k,Name,Value)
[idx,C] = kmeans(_)
[idx,C,sumd] = kmeans(_)
[idx,C,sumd,D] = kmeans(_)
X :N*P的数据矩阵
K: 表示将X划分为几类,为整数
Idx :N*1的向量,存储的是每个点的聚类标号
C: K*P的矩阵,存储的是K个聚类质心位置
sumd 1*K的和向量,存储的是类间所有点与该类质心点距离之和
D N*K的矩阵,存储的是每个点与所有质心的距离
[…]=Kmeans(…,’Param1’,Val1,’Param2’,Val2,…)
这其中的参数Param1、Param2等,主要可以设置为如下:
- ‘Distance’(距离度量)
‘sqEuclidean’ 欧式距离(默认时,采用此距离方式)
‘cityblock’ 绝度误差和,又称:L1
‘cosine’ 针对向量
‘correlation’ 针对有时序关系的值
‘Hamming’ 只针对二进制数据 - ‘Start’(初始质心位置选择方法)
‘sample’ 从X中随机选取K个质心点
‘uniform’ 根据X的分布范围均匀的随机生成K个质心
‘cluster’ 初始聚类阶段随机选择10%的X的子样本(此方法初始使用’sample’方法)
matrix 提供一K*P的矩阵,作为初始质心位置集合 - ‘Replicates’(聚类重复次数)
整数
%第一类数据
mu1=[0 0 0]; %均值
S1=[0.3 0 0;0 0.35 0;0 0 0.3]; %协方差
data1=mvnrnd(mu1,S1,100); %产生高斯分布数据
%%第二类数据
mu2=[1.25 1.25 1.25];
S2=[0.3 0 0;0 0.35 0;0 0 0.3];
data2=mvnrnd(mu2,S2,100);
%第三个类数据
mu3=[-1.25 1.25 -1.25];
S3=[0.3 0 0;0 0.35 0;0 0 0.3];
data3=mvnrnd(mu3,S3,100);
%显示数据
figure, hold on;
plot3(data1(:,1),data1(:,2),data1(:,3),'b+');
hold on;
plot3(data2(:,1),data2(:,2),data2(:,3),'r+');
plot3(data3(:,1),data3(:,2),data3(:,3),'g+');
grid on;
%三类数据合成一个不带标号的数据类
data=[data1;data2;data3]; %这里的data是不带标号的,即没有分好类
k = 3;
[idx, c] = kmeans(data, k); %分成3份
figure, hold on;
for i=1:300
if idx(i) == 1
plot3(data(i,1),data(i,2),data(i,3), 'b+');
elseif idx(i) == 2
plot3(data(i,1),data(i,2),data(i,3), 'r+');
else
plot3(data(i,1),data(i,2),data(i,3), 'g+');
end
end
for i=1:k
plot3(c(i, 1), c(i, 2), c(i, 3), 'ko'); % 画出中心点
end
grid on;
- 输入:
- 输出:
图中原点就是簇的中心。
附赠一个plot的point style 图
学习向量量化(LVQ)
与k均值算法类似,LVQ也是基于原型的聚类,但与一般的不同,LVQ假设数据样本带有类别标记,学习过程利用样本的这些监督信息来辅助聚类。
这个是用来建立神经网络的- -。
算法描述
newlvq(P,S1,PC,LR,LF) 建立LVQ向量化神经网络
P为一个Rx2维的网络输入向量取值范围的矩阵[Pmin Pmax];
Sl表示隐含层神经元的数目;
PC表示在第二层的权值中列所属类别的百分比;
LR表示学习速率,默认值为0.01;
Lf表示学习函数,默认值为learnlv1。
P = [-3 -2 -2 0 0 0 0 +2 +2 +3; ...
0 +1 -1 +2 +1 -1 -2 +1 -1 0];
Tc = [1 1 1 2 2 2 2 1 1 1];
T = ind2vec(Tc);
net = newlvq(P,4,[.6 .4]);
net = train(net,P,T);
- 输入:
- 输出:
高斯混合聚类(EM算法)
利用概率模型来表达聚类原型。
这个算法会算法每个样本可能在每个簇的概率。
算法流程
%第一类数据
mu1=[0 0 0]; %均值
S1=[0.3 0 0;0 0.35 0;0 0 0.3]; %协方差
data1=mvnrnd(mu1,S1,100); %产生高斯分布数据
%%第二类数据
mu2=[1.25 1.25 1.25];
S2=[0.3 0 0;0 0.35 0;0 0 0.3];
data2=mvnrnd(mu2,S2,100);
%第三个类数据
mu3=[-1.25 1.25 -1.25];
S3=[0.3 0 0;0 0.35 0;0 0 0.3];
data3=mvnrnd(mu3,S3,100);
%显示数据
figure, hold on;
plot3(data1(:,1),data1(:,2),data1(:,3),'b+');
hold on;
plot3(data2(:,1),data2(:,2),data2(:,3),'r+');
plot3(data3(:,1),data3(:,2),data3(:,3),'g+');
grid on;
%三类数据合成一个不带标号的数据类
data=[data1;data2;data3]; %这里的data是不带标号的,即没有分好类
k = 3;
[px, model] = gmm(data, k); %分成3份
figure, hold on;
for i=1:300
[m, idx] = max(px(i, :));
if idx == 1
plot3(data(i,1),data(i,2),data(i,3), 'b+');
elseif idx == 2
plot3(data(i,1),data(i,2),data(i,3), 'r+');
else
plot3(data(i,1),data(i,2),data(i,3), 'g+');
end
end
c = model.Miu;
for i=1:k
plot3(c(i, 1), c(i, 2), c(i, 3), 'ko'); % 画出中心点
end
grid on;
- 输入:
- 输出:
密度聚类(基于密度的聚类)
此类算法假设聚类结构能通过样本分布的紧密程度确定。
原理
从样本密度的角度来考察样本之间的可连接性,并基于可连接样本不断扩展聚类簇以获得最终的聚类结果。
DBSCAN
基于上面的条件,这个算法给出了下面几个概念
算法流程
层次聚类
层次聚类视图在不同层次对数据集进行划分,从而形成属性的聚类结构。
数据集的划分可以采用自底向上的聚合策略,也可采用自顶向下的分拆策略。
AGNES
AGNES是采用自底向上聚合策略的层次聚类算法。
算法流程
clusterdata(X, cutoff)
X = [gallery('uniformdata',[10 3],12);...
gallery('uniformdata',[10 3],13)+1.2;...
gallery('uniformdata',[10 3],14)+2.5];
subplot(1, 2, 1), plot3(X(:, 1), X(:, 2), X(:, 3), '*');
grid on;
T = clusterdata(X,'maxclust',3);
subplot(1, 2, 2), scatter3(X(:,1),X(:,2),X(:,3),100,T,'filled');
- 输入:
- 输出:
降维与度量学习
k近邻学习(kNN)
k-Nearest Neighbor,k近邻学习是一种常用的监督学习方法。
工作机制
给定测试样本,基于某种距离度量找出训练集中与其最靠近的k个训练样本,然后基于这k个邻居的信息来进行预测,通常,在分类任务中可使用投票法(选择这k个样本中出现最多的类别标记作为预测结果),在回归任务中可使用平均法(将这k个样本的实值输出标记的平均值作为预测结果)。
另外,还可以基于距离远近进行加权或加权投票,距离越近的样本权重越大。
这里的距离,跟聚类的距离是一个距离,感觉用距离的话,更加形象一点,但有的时候转不过来弯- -。
这里将k近邻学习概括起来,就一句话,近朱者赤,近墨者黑,你跟测试样本越近,就越可能是一个类。
注意,这里要理解,测试样本就是一个样本,我们需要得出这个样本的标记(或者输出标记的平均值),就靠与ta靠近的这k个训练样本,他们决定了测试样本的标记到底是什么。
看到这里,你一下子就能知道这个算法的特点,就是根本不需要进行显式的训练,就是不需要花费时间来进行训练,有这样的特点的学习算法,就称之为懒惰学习,反之,则成为急切学习。
看上面这个图,从近到远,依次选择1 3 5,根据投票法进行分类,得到的结果竟然不同,所以呢,这个算法,k值很重要,起到决定性的影响。
不过我感觉,这个算法的好坏,错误率的高低,取决于距离度量算法- -。
计算错误率
再看一下下面这个讨论
根据这个,k近邻学习,如果是临时性创建分类或者回归,并且精度不要求太精准的话(用的次数越多,精度越高),可以使用,不过有一点需要注意,就是在维度高,数据量小,密度小的时候,错误率会变高。
低维嵌入
当数据的维度变高时,计算量几何增长,使得运算变得困难,不易得出结果(这种情况称为维数灾难)。
一个重要的解决途径就是降维(维数化简)。
降维
通过某种数学变换将原始高维属性空间转变为一个低维子空间,在这个子空间中样本密度大幅提高,距离计算也变得更为容易。
当然,这个降维其实也有他的理论依据,不然能随随便便降么- -
很多时候,人们或者计算机搜集到的数据样本虽然是高维的,但是跟学习任务相关的也许仅仅是某个低维分布,比如你在做图像处理的时候,我们只需要判断有几个独立的物体,这个时候只需要将图片计算阈值,进行二值处理即可,即降维(忽略颜色)。
多维缩放(MDS)
在现实应用中为了有效降维,往往仅需要降维后的距离与原始空间中的距离尽可能接近,而不必严格相等。
算法流程
主成分分析(PCA)
算法流程
复习训练
直方图均衡化
直方图均衡化主要用于增强灰度值动态范围偏小的图像的对比度。
什么叫做动态范围,表示图像中所包含的从“最暗”至“最亮”的范围,动态范围越大,所能表现的层次越丰富,所包含的色彩空间也越广。
直方图
直方图是图像的一种统计表达形式。
直方图可以在一定程度上反映图像的概貌性描述,包括图像的灰度范围、灰度分布、整幅图像的亮度均值和阴暗对比度等
f = imread('coins.png');
g = imhist(f, 25);
subplot(1, 2, 1), imshow(f);
subplot(1, 2, 2), bar(g);
上面这个直方图的代表意义就是将灰度分为25份,这张图片中,每份所占的数量(可以看出来,这个图偏暗)。
直方图均衡化
把原始图像的灰度统计直方图变换为均匀分布的形式,这样就增加了像素灰度值的动态范围,从而达到增强图像整体对比度的效果(就是将直方图的每个小柱子均匀分布,所以也叫直方图均匀化)。
主要步骤
- 预处理
输入图像,计算该图像直方图。 - 灰度变换表
根据输入图像的直方图计算灰度值变化表。 - 查表变换
将预处理得到的直方图使用灰度变换表进行查表变换操作,通过遍历每一个像元(像素点),将原始图像灰度值x放入交换表中,可得到该点的新灰度值。
以上面那个图为例
f = imread('coins.png');
g = imhist(f);
subplot(2, 2, 1), imshow(f);
subplot(2, 2, 2), bar(g);
h = histeq(f); % 全局直方图均衡化
g = imhist(h);
subplot(2, 2, 3), imshow(h);
subplot(2, 2, 4), bar(g);
效果还是很明显的,但是在这里,我们可以看到,图片变得更加亮,显得更加真实,细节更多,但是有部分细节已经消失了,主要是这张图,不应该用这个方法进行处理,这张图应该用对比度增强adapthisteq(这个函数其实是局部直方图均衡化-使用CLAHE算法,也叫限制对比度自适应直方图均衡)进行处理(或者用线性滤波器)。
f = imread('coins.png');
subplot(1, 2, 1), imshow(f), title('原图');
g = adapthisteq(f);
subplot(1, 2, 2), imshow(g), title('增强后');
上面这个图就是使用了对比度增强,有点跑题了- -。
f = imread('天安门.jpg');
r = f(:, :, 1);
g = f(:, :, 2);
b = f(:, :, 3);
g = cat(3, histeq(r), histeq(g), histeq(b));
subplot(1, 2, 1), imshow(f);
subplot(1, 2, 2), imshow(g);
效果明显的看下面这张图
f = imread('天安门.jpg');
hsi = rgb2hsi(f);
h = hsi (:, :, 1);
s = hsi (:, :, 2);
i = hsi (:, :, 3);
g = cat(3, h, s, histeq(i));
g = hsi2rgb(g);
subplot(1, 2, 1), imshow(f);
subplot(1, 2, 2), imshow(g);
可以看出来,这张图的下半部分,去雾效果就非常明显了(是不是很鲜活,很神奇,其实就是将细节放大),但是上半部分出现了光圈- -,
我们只需要使用adapthisteq局部直方图均衡化即可,直接看效果
f = imread('天安门.jpg');
hsi = rgb2hsi(f);
h = hsi (:, :, 1);
s = hsi (:, :, 2);
i = hsi (:, :, 3);
g = cat(3, h, s, adapthisteq(i));
g = hsi2rgb(g);
subplot(1, 2, 1), imshow(f);
subplot(1, 2, 2), imshow(g);
可以看到,光圈没有了- -
adapthisteq这个函数还有很多参数,可以直接help adapthisteq查看帮助文档,不过直接用即可adapthisteq(f)。
不过呢,上面的全局均衡化,局部均衡化,只是简单的进行处理,得到的图像偏暗,存在模糊的地方。
我们可以使用imfilter线性滤波来进行处理
f = imread('天安门.jpg');
hsi = rgb2hsi(f);
h = hsi (:, :, 1);
s = hsi (:, :, 2);
i = hsi (:, :, 3);
g = cat(3, h, s, adapthisteq(imfilter(i, fspecial('average'))));
g = hsi2rgb(g);
subplot(1, 2, 1), imshow(f);
subplot(1, 2, 2), imshow(g);
效果不是很明显,这个时候我们就可以使用Retinex增强处理的MSR。
function in = retinex(f)
% =========================
% 使用Retinex增强算法来对图像进行增强
% f = retinex(f);
% 注意传入的f必须是rgb格式的
% =========================
% 先提取rgb通道
fr = f(:, :, 1);
fg = f(:, :, 2);
fb = f(:, :, 3);
% 将数据归一化
mr = mat2gray(im2double(fr));
mg = mat2gray(im2double(fg));
mb = mat2gray(im2double(fb));
alpha = 15 ^ 2 / 2;
n = 125;
n1 = floor((n + 1) / 2);
for i = 1 : n
for j = 1 : n
b(i, j) = exp(-((i - n1) ^ 2 + (j - n1) ^ 2) / (4 * alpha)) / (pi * alpha);
end
end
% 卷积滤波
nr1 = imfilter(mr, b, 'conv', 'replicate');
ng1 = imfilter(mg, b, 'conv', 'replicate');
nb1 = imfilter(mb, b, 'conv', 'replicate');
ur1 = log(nr1);
ug1 = log(ng1);
ub1 = log(nb1);
tr1 = log(mr);
tg1 = log(mg);
tb1 = log(mb);
yr1 = (tr1 - ur1) / 3;
yg1 = (tg1 - ug1) / 3;
yb1 = (tb1 - ub1) / 3;
alpha = 80 ^ 2 / 2;
% n = 31;
for i = 1 : n
for j = 1 : n
a(i, j) = exp(-((i - n1) ^ 2 + (j - n1) ^ 2) / (4 * alpha)) / (6 * pi * alpha);
end
end
% 卷积滤波
nr2 = imfilter(mr, a, 'conv', 'replicate');
ng2 = imfilter(mg, a, 'conv', 'replicate');
nb2 = imfilter(mb, a, 'conv', 'replicate');
ur2 = log(nr2);
ug2 = log(ng2);
ub2 = log(nb2);
tr2 = log(mr);
tg2 = log(mg);
tb2 = log(mb);
yr2 = (tr2 - ur2) / 3;
yg2 = (tg2 - ug2) / 3;
yb2 = (tb2 - ub2) / 3;
imshow(yr2);
alpha = 120 ^ 2 / 2;
% n = 501;
for i = 1 : n
for j = 1 : n
e(i, j) = exp(-((i - n1) ^ 2 + (j - n1) ^ 2) / (4 * alpha)) / (4 * pi * alpha);
end
end
% 卷积滤波
nr3 = imfilter(mr, e, 'conv', 'replicate');
ng3 = imfilter(mg, e, 'conv', 'replicate');
nb3 = imfilter(mb, e, 'conv', 'replicate');
ur3 = log(nr3);
ug3 = log(ng3);
ub3 = log(nb3);
tr3 = log(mr);
tg3 = log(mg);
tb3 = log(mb);
yr3 = (tr3 - ur3) / 3;
yg3 = (tg3 - ug3) / 3;
yb3 = (tb3 - ub3) / 3;
dr = yr1 + yr2 + yr3;
dg = yg1 + yg2 + yg3;
db = yb1 + yb2 + yb3;
cr = im2uint8(dr);
cg = im2uint8(dg);
cb = im2uint8(db);
in = cat(3, cr, cg, cb);
end
上面这个是Retinex的MSR函数
f = imread('天安门.jpg');
subplot(1, 2, 1), imshow(f), title('原图');
subplot(1, 2, 2), imshow(retinex(f)), title('MSR处理后');
怎么样,是不是超级明显- -,如果嫌太亮,可以将滤大小调高。
参考文献
1. 《MATLAB计算机视觉与深度学习实战》 刘衍琦、詹福宇、蒋献文、周华英. 电子工业出版社2017年6月出版