f*******a
发帖数: 663 | 1 原帖见
http://www.mitbbs.com/article_t/DataSciences/6761.html
一点分析心得，与大家共享，以期抛砖引玉。
感谢zhaoce的总结一文让我看到这篇文章；也非常感谢f0008朋友在我始终无法下载附
件的情况下把附件发给了我。
===========================================================================
聚类算法能利用的一般是局部特性，如邻域点距离、基于核函数的密度估计。Mean-
shift算法就是一个非常经典的算法，以梯度方法迭代至局部密度峰值点。
这个算法的思路其实与Mean-shift很类似，虽然作者要在文章中反复说和Mean-shift不
一样，但本质上非常相近。MS以梯度寻找峰值点，而这个算法则是直接在点群中搜索峰
值点。这样做是基于一个近似假设：峰值点和点群中的某点距离不远。举个极端的例子
：只有一类，20个点均匀分布在一个圆上。MS算法可以准确聚到圆心，只要核大小足够
。而这个算法，只能聚在这20个点中的某个点上。这个假设在一般情况下可以接受，也
能保持足够的精度。在这个假设下，以简单穷举就可以搜到局部峰值点。
具体做法：
1）对每点A，计算核函数下的密度估计rho(A)。
2）对每点A，寻找密度大于rho(A)且距离最近的点B，定义 delta(A) = dist(A, B)，
并记录映射关系Label(A) = Label(B)。
3）画一个rho- delta图，手工选取门限(rho_min, delta_min)，确定峰值点。根据映
射关系回溯各点的Label。
4）对于每个cluster，计算边界上点的平均rho，低于此值的则认为是outlier
需要注意的是：公式（1）用了一个截断核，把rho定义成邻域点个数，而且特地加以解
释。这是为了区别MS，故意把读者带到坑里去。这样做有两个问题：
- Fig.1B中的rho值明显非离散值；
- 如果一个类中，中心位置有2个点正好有同样多的邻域数，不合并的话那么必然崩
掉。
看了代码之后就很明确了，rho就是基于核函数的密度估计。
这个算法除了原始数据输入外，需要指定3个参数：核尺寸dc以及聚类中心检出参数(
rho_min, delta_min)。后两个参数在确定dc后设置有一定的弹性，所以稍微改改有可
能只要设置dc，而这个一般由数据特性决定比较好确定。
算法总的来说性能不错，参见各结果图，其中有几张dc没有列出来用的是默认值。但也
没有看上去那么美：
1）聚类中心检出这个太灵活，现在是只靠主观判定，碰到难点的数据可能不太好选（
如Fig. 4A）。
2）dc的选择其实很重要，而不是像补充材料中显示的那么鲁棒。补充材料中的Fig. S1
是没有去噪的结果，而且实际跑的结果远没这么好，特别是在核尺寸较大的时候。
3）一个类内出现两个点的时候结果就废了，大部分都被被认为是噪点。
这篇文章的贡献，我觉得是有争议的，从审到接受用了快一年也能看出来。在聚类中心
近似假设的基础上，用搜索替代了梯度求解，这个贡献就见仁见智了。优点是克服了MS
的各向同性问题，而且具备较好的抗噪能力。
感觉审稿人未必很了解聚类算法，除去公式（1）和图1的明显不符，最后的比较结果只
是和Kmeans对比，这个真没什么说服力。
代码功力一般，可能是别的语言转过来的，Matlab的基本操作不熟，运算效率低。 | f*******a
发帖数: 663 | 2 高斯核函数的测试结果
【在 f*******a 的大作中提到】
: 原帖见
: http://www.mitbbs.com/article_t/DataSciences/6761.html
: 一点分析心得，与大家共享，以期抛砖引玉。
: 感谢zhaoce的总结一文让我看到这篇文章；也非常感谢f0008朋友在我始终无法下载附
: 件的情况下把附件发给了我。
: ===========================================================================
: 聚类算法能利用的一般是局部特性，如邻域点距离、基于核函数的密度估计。Mean-
: shift算法就是一个非常经典的算法，以梯度方法迭代至局部密度峰值点。
: 这个算法的思路其实与Mean-shift很类似，虽然作者要在文章中反复说和Mean-shift不
: 一样，但本质上非常相近。MS以梯度寻找峰值点，而这个算法则是直接在点群中搜索峰
| f*******a
发帖数: 663 | 3 截断核函数的测试结果
【在 f*******a 的大作中提到】
: 原帖见
: http://www.mitbbs.com/article_t/DataSciences/6761.html
: 一点分析心得，与大家共享，以期抛砖引玉。
: 感谢zhaoce的总结一文让我看到这篇文章；也非常感谢f0008朋友在我始终无法下载附
: 件的情况下把附件发给了我。
: ===========================================================================
: 聚类算法能利用的一般是局部特性，如邻域点距离、基于核函数的密度估计。Mean-
: shift算法就是一个非常经典的算法，以梯度方法迭代至局部密度峰值点。
: 这个算法的思路其实与Mean-shift很类似，虽然作者要在文章中反复说和Mean-shift不
: 一样，但本质上非常相近。MS以梯度寻找峰值点，而这个算法则是直接在点群中搜索峰
| g*********n
发帖数: 119 | | f***8
发帖数: 571 | | x**l
发帖数: 2337 | 6 哪位能解释一下decision graph那个图啊，看不懂呢。谢谢了。 | f*******a
发帖数: 663 | 7 开始忘贴代码了，有朋友要求，就把修改后的代码贴在这里。改动不多，可以部分提升
效率。原来的也没删，注释掉了。供参考。
=========================================================================
clear;
close all
disp('The only input needed is a distance matrix file')
disp('The format of this file should be: ')
disp('Column 1: id of element i')
disp('Column 2: id of element j')
disp('Column 3: dist(i,j)')
if(0)
% mdist=input('name of the distance matrix file (with single quotes)?\n'
);
mdist = 'example_distances.dat';
disp('Reading input distance matrix')
xx=load(mdist);
N=size(xx,1);
ND=max(xx(:,2));
NL=max(xx(:,1));
if(NL>ND)
ND=NL;
end
dist = zeros(ND, ND);
% for i=1:ND
% for j=1:ND
% dist(i,j)=0;
% end
% end
for i=1:N
ii=xx(i,1);
jj=xx(i,2);
dist(ii,jj)=xx(i,3);
dist(jj,ii)=xx(i,3);
end
percent=2.0;
fprintf('average percentage of neighbours (hard coded): %5.6f\n',
percent);
position=round(N*percent/100);
sda=sort(xx(:,3));
dc=sda(position);
clear xx;
save('PeakCluData.mat');
return;
else
load('PeakCluData.mat');
end;
dc = dc * 3;
dc = 0.1;
% rhomin = 20;
% deltamin = 0.1;
rhomin = 10;
deltamin = 0.15;
percent=2.0;
fprintf('average percentage of neighbours (hard coded): %5.6f\n', percent);
fprintf('Computing Rho with gaussian kernel of radius: %12.6f\n', dc);
rho = zeros(1, ND);
% for i=1:ND
% rho(i)=0.;
% end
% %
% % Gaussian kernel
% %
% for i=1:ND-1
% for j=i+1:ND
% rho(i)=rho(i)+exp(-(dist(i,j)/dc)*(dist(i,j)/dc));
% rho(j)=rho(j)+exp(-(dist(i,j)/dc)*(dist(i,j)/dc));
% end
% end
% Kernel distance
k_dist = exp(-dist.^2/dc^2);
rho = sum(k_dist) - 1;
% %
% % "Cut off" kernel
% %
% for i=1:ND-1
% for j=i+1:ND
% if (dist(i,j)
% rho(i)=rho(i)+1.;
% rho(j)=rho(j)+1.;
% end
% end
% end
maxd=max(max(dist));
% [rho_sorted,ordrho]=sort(rho,'descend');
[rho_sorted,ordrho]=sort(-rho);
rho_sorted = -rho_sorted;
delta(ordrho(1))=-1.;
nneigh(ordrho(1))=0;
% Find local peak among neighbors; mark the clustering trace
for ii=2:ND
delta(ordrho(ii))=maxd;
for jj=1:ii-1
if(dist(ordrho(ii),ordrho(jj))
delta(ordrho(ii))=dist(ordrho(ii),ordrho(jj));
nneigh(ordrho(ii))=ordrho(jj);
end
end
end
delta(ordrho(1))=max(delta(:));
disp('Generated file:DECISION GRAPH')
disp('column 1:Density')
disp('column 2:Delta')
fid = fopen('DECISION_GRAPH', 'w');
for i=1:ND
fprintf(fid, '%6.2f %6.2f\n', rho(i),delta(i));
end
disp('Select a rectangle enclosing cluster centers')
scrsz = get(0,'ScreenSize');
%figure('Position',[6 72 scrsz(3)/4. scrsz(4)/1.3]);
figure('Position',[scrsz(3)/4 scrsz(4)/4 scrsz(3)/2. scrsz(4)/2]);
for i=1:ND
ind(i)=i;
gamma(i)=rho(i)*delta(i);
end
subplot(1,2,1)
tt=plot(rho(:),delta(:),'o','MarkerSize',5,'MarkerFaceColor','k','
MarkerEdgeColor','k');
title(sprintf('Decision Graph: dc = %.3f', dc));
xlabel ('\rho')
ylabel ('\delta')
subplot(1,2,1)
% Read the threshold set manually
% rect = getrect(1);
% rhomin=rect(1);
% deltamin=rect(4);
disp(sprintf('rhomin = %f, deltamin = %f', rhomin, deltamin));
NCLUST=0;
for i=1:ND
cl(i)=-1;
end
for i=1:ND
if ( (rho(i)>rhomin) && (delta(i)>deltamin))
NCLUST=NCLUST+1;
cl(i)=NCLUST;
icl(NCLUST)=i;
end
end
fprintf('NUMBER OF CLUSTERS: %i \n', NCLUST);
disp('Performing assignation')
%assignation
for i=1:ND
if (cl(ordrho(i))==-1)
cl(ordrho(i))=cl(nneigh(ordrho(i)));
end
end
%halo
% for i=1:ND
% halo(i)=cl(i);
% end
halo = cl;
if (NCLUST>1)
% for i=1:NCLUST
% bord_rho(i)=0.;
% end
bord_rho = zeros(1, NCLUST);
% For each cluster, calculate the average rho of border points
for i=1:ND-1
for j=i+1:ND
if ((cl(i)~=cl(j))&& (dist(i,j)<=dc))
rho_aver=(rho(i)+rho(j))/2.;
if (rho_aver>bord_rho(cl(i)))
bord_rho(cl(i))=rho_aver;
end
if (rho_aver>bord_rho(cl(j)))
bord_rho(cl(j))=rho_aver;
end
end
end
end
% Points whose rho is less than the average border rho will be considered
% as noise
for i=1:ND
if (rho(i)
halo(i)=0;
end
end
end
for i=1:NCLUST
nc=0;
nh=0;
for j=1:ND
if (cl(j)==i)
nc=nc+1;
end
if (halo(j)==i)
nh=nh+1;
end
end
fprintf('CLUSTER: %i CENTER: %i ELEMENTS: %i CORE: %i HALO: %i \n', i,icl(
i),nc,nh,nc-nh);
end
cmap=colormap;
for i=1:NCLUST
ic=int8((i*64.)/(NCLUST*1.));
subplot(1,2,1)
hold on
plot(rho(icl(i)),delta(icl(i)),'o','MarkerSize',8,'MarkerFaceColor',cmap(
ic,:),'MarkerEdgeColor',cmap(ic,:));
end
subplot(1,2,2)
disp('Performing 2D nonclassical multidimensional scaling')
Y1 = mdscale(dist, 2, 'criterion','metricstress');
plot(Y1(:,1),Y1(:,2),'o','MarkerSize',2,'MarkerFaceColor','k','
MarkerEdgeColor','k');
%title ('2D Nonclassical multidimensional scaling','FontSize',15.0)
title (['Clu result: \rho_{min} = ', sprintf('%.2f', rhomin) ', \delta_{min}
= ' sprintf('%.2f', deltamin) ]);
xlabel ('X')
ylabel ('Y')
axis equal;
for i=1:ND
A(i,1)=0.;
A(i,2)=0.;
end
for i=1:NCLUST
nn=0;
ic=int8((i*64.)/(NCLUST*1.));
for j=1:ND
if (halo(j)==i)
nn=nn+1;
A(nn,1)=Y1(j,1);
A(nn,2)=Y1(j,2);
end
end
hold on
plot(A(1:nn,1),A(1:nn,2),'o','MarkerSize',2,'MarkerFaceColor',cmap(ic,:),'
MarkerEdgeColor',cmap(ic,:));
end
%for i=1:ND
% if (halo(i)>0)
% ic=int8((halo(i)*64.)/(NCLUST*1.));
% hold on
% plot(Y1(i,1),Y1(i,2),'o','MarkerSize',2,'MarkerFaceColor',cmap(ic,:),
'MarkerEdgeColor',cmap(ic,:));
% end
%end
faa = fopen('CLUSTER_ASSIGNATION', 'w');
disp('Generated file:CLUSTER_ASSIGNATION')
disp('column 1:element id')
disp('column 2:cluster assignation without halo control')
disp('column 3:cluster assignation with halo control')
for i=1:ND
fprintf(faa, '%i %i %i\n',i,cl(i),halo(i));
end
im = getframe(gcf);
imwrite(im.cdata, 'PeakCluRes.png'); | a*******k
发帖数: 261 | | c****t
发帖数: 19049 | | n*****3
发帖数: 1584 | 10 matlab?
【在 c****t 的大作中提到】
: 赞！这么踏实作功课的很少见啊
| N******n
发帖数: 3003 | 11 这个文章发到science上，看着很新颖，不过好像类似的以前发到JMachine learning
Nonparametric statistical approach to clustering via mode identification
主要是计算kernal distensity,来确定mode,通过EM算法实现。 | N******n
发帖数: 3003 | 12 新颖之处可能在：
一个是用图表示出可能几个cluster和区分outlier.
其他的好像没有区别。 | j**********3
发帖数: 3211 | 13 我一直以为zhaoce大牛只是java大牛，想不到ds也有的参与 | T*****u
发帖数: 7103 | |
|
