博文
三种常见的聚类算法的python实现 kmeans、Hierarchical clusteri
||
聚类是机器学习、数据挖掘相关的一类很常见的问题。关于聚类算法的介绍这里就不多写了，因为无论是教科书还是网络上都有太多的资料了。这里，用一个《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. See more : http://en.wikipedia.org/wiki/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://blog.sciencenet.cn/blog-1271266-855059.html
上一篇：甲基化DNA免疫共沉淀测序原理及简介（MeDIP-seq）
下一篇：全基因组重亚硫酸盐测序（whole genome bisulfite sequencing）
