juan503949201的个人博客分享 http://blog.sciencenet.cn/u/juan503949201

博文

三种常见的聚类算法的python实现 kmeans、Hierarchical clusteri

已有 10661 次阅读 2014-12-30 17:23 |个人分类:Machine Learning|系统分类:科研笔记

聚类是机器学习、数据挖掘相关的一类很常见的问题。关于聚类算法的介绍这里就不多写了,因为无论是教科书还是网络上都有太多的资料了。这里,用一个《Programming Collective Intelligence》中的聚类例子,写几个经典聚类算法的实现,分别是hierachiclaCluster、kmeans、kmedoids。

  另外,最近一直在看数据挖掘、自然语言处理相关的东西,通过看资料发现有些东西很好理解,但是长时间不用的话,过一段时间就忘记得差不多了。其实换个角度想,也是自己对这些东西理解得不深刻。我觉得踏踏实实的实现这些算法是必要的,因为在实现过程中,我们可以发现哪些地方是算法的核心思想,哪些地方是性能瓶颈,这也为进一步优化提供了基础。

 

问题背景及数据集:

    blog的聚类。如何判断两个blog是相似的呢?由于相关的博客的主题内容应该是相同的,即会出现很多相同的词。于是我们为每个blog定义一个向量,该向量的维度为数据集中所有出现的不同词的个数,向量的值为对应词出现的次数。判断两个向量间的相似性,我们使用pearson相似性度量算法。

   测试数据集描述:blogdata.txt数据集,该数据集第一行包含所有数据集中出现的不同词,共有m个,剩下所有行,每一行对应一篇blog,共m+1列,用tab分开。第一列为博客名,接下来m列为一个向量,代码每个词出现的次数。接下来就聚类吧~

Hierachical clustering :

维基百科:http://en.wikipedia.org/wiki/Hierarchical_clustering

kmeans clustering :

维基百科:http://en.wikipedia.org/wiki/Kmeans

kmedoids clustering :

维基百科:http://en.wikipedia.org/wiki/K-medoids

虽然上面三种算法都很好理解,但是这都是基础算法,要想深入,还有很多很多相关问题需要解决,比如k如何设置;随机选取初始点的问题等等,而且如何选取好用的聚类算法也值得商榷。

github代码位置:https://github.com/LixinZhang/bookreviews/tree/master/Programming_Collective_Intelligence/chapter3

clusterBase.py 用来导入数据

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
from math import sqrt
 
def importData(FIFE = 'blogdata.txt') :
   blogwords = []
   blognames = []
   f = open(FIFE, 'r')
   words = f.readline().split()
   #//remove 'rn'
   for line in f:    
       blog = line[:-2].split('t')
       blognames.append(blog[0])        
       blogwords.append([int(word_c) for word_c in blog[1:]]       )
   return blogwords,blognames
 
 
def pearson_distance(vector1, vector2) :
   """
   Calculate distance between two vectors using pearson method
   """
   sum1 = sum(vector1)
   sum2 = sum(vector2)
 
   sum1Sq = sum([pow(v,2) for v in vector1])
   sum2Sq = sum([pow(v,2) for v in vector2])
 
   pSum = sum([vector1[i] * vector2[i] for i in range(len(vector1))])
 
   num = pSum - (sum1*sum2/len(vector1))
   den = sqrt((sum1Sq - pow(sum1,2)/len(vector1)) * (sum2Sq - pow(sum2,2)/len(vector1)))
 
   if den == 0 : return 0.0
   return 1.0 - num/den

hierachiclaCluster.py

注意distances字典的使用,可以减少大量重复的计算

另外这个聚类算法,最终生成的是一颗树形结构。因此打印结果,需要递归的进行

在集体智慧编程书里,给出了一个利用PIL模块画树形图的算法,挺有参考价值的。

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
#/usr/bin/python
from clusterBase import importData,pearson_distance
 
class bicluster:
   def __init__(self, vec, left=None,right=None,distance=0.0,id=None) :
       self.left = left
       self.right = right
       self.vec = vec
       self.id = id
       self.distance = distance
 
def hcluster(blogwords,blognames) :
   biclusters = [ bicluster(vec = blogwords[i], id = i ) for i in range(len(blogwords)) ]
   distances = {}
   flag = None;
   currentclusted = -1
   while(len(biclusters) > 1) :
       min_val = 2;
       biclusters_len = len(biclusters)
       for i in range(biclusters_len-1) :
           for j in range(i + 1, biclusters_len) :
               if distances.get((biclusters[i].id,biclusters[j].id)) == None:
                   distances[(biclusters[i].id,biclusters[j].id)] = pearson_distance(biclusters[i].vec,biclusters[j].vec)
               d = distances[(biclusters[i].id,biclusters[j].id)]
               if d < min_val :
                   min_val = d
                   flag = (i,j)
       bic1,bic2 = flag
       newvec = [(biclusters[bic1].vec[i] + biclusters[bic2].vec[i])/2 for i in range(len(biclusters[bic1].vec))]
       newbic = bicluster(newvec, left=biclusters[bic1], right=biclusters[bic2], distance=min_val, id = currentclusted)
       currentclusted -= 1
       del biclusters[bic2]
       del biclusters[bic1]
       biclusters.append(newbic)
   return biclusters[0]
 
'''
Print the tree structure, save as a jpeg image file.
'''
 
from PIL import Image, ImageDraw
 
def getheight (clust) :
   if clust.left == None and clust.right == None :return 1
   return getheight(clust.left) + getheight(clust.right)
 
def getdepth(clust) :
   if clust.left == None and clust.right == None :return 0
   return max(getdepth(clust.left),getdepth(clust.right)) + clust.distance
 
def drawdendrogram(clust,labels,jpeg='clusters.jpg') :
   h = getheight(clust) * 20
   w = 1200
   depth = getdepth(clust)
   scaling = float(w-150) / depth
    
   img = Image.new('RGB',(w,h),(255,255,255))
   draw = ImageDraw.Draw(img)
   draw.line((0,h/2,10,h/2),fill=(255,0,0))
 
   drawnode(draw,clust,10,(h/2),scaling,labels)
   img.save(jpeg,'JPEG')
 
def drawnode(draw, clust, x, y, scaling, labels) :
   if clust.id < 0 :
       h1 = getheight(clust.left) * 20
       h2 = getheight(clust.right) * 20
       top = y - (h1+h2)/2
       bottom = y + (h1+h2)/2
       #line length
       ll = clust.distance * scaling
       draw.line((x,top+h1/2,x,bottom-h2/2),fill=(255,0,0))
       draw.line((x,top+h1/2,x+ll,top+h1/2),fill=(255,0,0))
       draw.line((x,bottom-h2/2,x+ll,bottom-h2/2),fill=(255,0,0))
       drawnode(draw,clust.left, x+ll, top+h1/2,scaling, labels)
       drawnode(draw,clust.right,x+ll, bottom-h2/2,scaling,labels)
   else :
       draw.text((x+5,y-7),labels[clust.id],(0,0,0))
 
if __name__ == '__main__' :
   #pearson_distance
   blogwords,blognames = importData()
   clust = hcluster(blogwords,blognames)
   print clust
   #drawdendrogram(clust,blognames,jpeg='blogclust.jpg')

 

kmeans.py

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
from clusterBase import importData, pearson_distance
import random
 
def print_matchs(matchs) :
   for i in range(len(matchs)) :
       print i , '---->',
       for item in matchs[i] :
           print item,
       print
   print '-'*20
 
def kmeans(blogwords, k) :
   min_max_per_word = [ [min([row[i] for row in blogwords]), max([row[i] for row in blogwords])]  for i in range(len(blogwords[0]))]
   # generate k clusters randomly
   clusters = []
   for i in range(k) :
       cluster = []
       for min_, max_ in min_max_per_word :
           cluster.append(random.random() * (max_ - min_) + min_)
       clusters.append(cluster)
   lables = []
   matchs = [ [] for i in range(k)]
   lastmatchs = [ [] for i in range(k)]
 
   rounds = 100
   while rounds > 0 :
       matchs = [ [] for i in range(k)]
       print 'round t',rounds
       for i in range(len(blogwords)) :
           bestmatch_cluster = None
           min_distance = 2.1
           for j in range(k) :
               dis = pearson_distance(clusters[j], blogwords[i])
               if dis < min_distance :
                   min_distance = dis
                   bestmatch_cluster = j
           matchs[bestmatch_cluster].append(i)
       print_matchs(matchs)
       print_matchs(lastmatchs)
       if matchs == lastmatchs : break
       lastmatchs = [[ item for item in matchs[i] ] for i in range(k)]
       #move the centroids to the average of their members
       for j in range(k) :
           avg = [0.0 for i in range(len(blogwords[0])) ]
           for m in matchs[j] :
               vec = blogwords[m]
               for i in range(len(blogwords[0])) :
                   avg[i] += vec[i]
           avg = [ item / len(matchs[j]) for item in avg]
           clusters[j] = avg
       rounds -= 1
 
        
if __name__ == '__main__' :
   blogwords,blognames = importData()
   kmeans(blogwords,5)

k-medoids.py

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
'''
Another clusting method , k-medoids.
The most common realisation of k-medoid clustering is the Partitioning Around Medoids (PAM) algorithm and is as follows:[2]
1. Initialize: randomly select k of the n data points as the medoids
2. Associate each data point to the closest medoid. ("closest" here is defined using any valid distance metric, most commonly Euclidean distance, Manhattan distance or Minkowski distance)
3. For each medoid m
    For each non-medoid data point o
        Swap m and o and compute the total cost of the configuration
4. Select the configuration with the lowest cost.
5. repeat steps 2 to 4 until there is no change in the medoid.
'''
from clusterBase import importData, pearson_distance
import random
 
distances_cache = {}
 
def totalcost(blogwords, costf, medoids_idx) :
   size = len(blogwords)
   total_cost = 0.0
   medoids = {}
   for idx in medoids_idx :
       medoids[idx] = []
   for i in range(size) :
       choice = None
       min_cost = 2.1
       for m in medoids :
           tmp = distances_cache.get((m,i),None)
           if tmp == None :
               tmp = pearson_distance(blogwords[m],blogwords[i])
               distances_cache[(m,i)] = tmp
           if tmp < min_cost :
               choice = m
               min_cost = tmp
       medoids[choice].append(i)
       total_cost += min_cost
   return total_cost, medoids
    
 
def kmedoids(blogwords, k) :
   size = len(blogwords)
   medoids_idx = random.sample([i for i in range(size)], k)
   pre_cost, medoids = totalcost(blogwords,pearson_distance,medoids_idx)
   print pre_cost
   current_cost = 2.1 * size # maxmum of pearson_distances is 2.    
   best_choice = []
   best_res = {}
   iter_count = 0
   while 1 :
       for m in medoids :
           for item in medoids[m] :
               if item != m :
                   idx = medoids_idx.index(m)
                   swap_temp = medoids_idx[idx]
                   medoids_idx[idx] = item
                   tmp,medoids_ = totalcost(blogwords,pearson_distance,medoids_idx)
                   #print tmp,'-------->',medoids_.keys()
                   if tmp < current_cost :
                       best_choice = list(medoids_idx)
                       best_res = dict(medoids_)
                       current_cost = tmp
                   medoids_idx[idx] = swap_temp
       iter_count += 1
       print current_cost,iter_count
       if best_choice == medoids_idx : break
       if current_cost <= pre_cost :
           pre_cost = current_cost
           medoids = best_res
           medoids_idx = best_choice
        
    
   return current_cost, best_choice, best_res
 
def print_match(best_medoids, blognames) :
   for medoid in best_medoids :
       print blognames[medoid],'----->',
       for m in best_medoids[medoid] :
           print '(',m,blognames[m],')',
       print
       print '---------' * 20
 
if __name__ == '__main__' :
   blogwords, blognames = importData()
   best_cost,best_choice,best_medoids = kmedoids(blogwords,8)
   print_match(best_medoids,blognames)




https://wap.sciencenet.cn/blog-1271266-855059.html

上一篇:甲基化DNA免疫共沉淀测序原理及简介(MeDIP-seq)
下一篇:全基因组重亚硫酸盐测序(whole genome bisulfite sequencing)

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2021-12-2 09:07

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部