编写层次聚类算法

编写层次聚类算法 - 图1

我们可以使用优先队列来实现这个聚类算法。

什么是优先队列呢?

普通的队列有“先进先出”的规则,比如向队列先后添加Moa、Suzuka、Yui,取出时得到的也是Moa、Suzuka、Yui:

编写层次聚类算法 - 图2

而对于优先队列,每个元素都可以附加一个优先级,从队列中取出时会得到优先级最高的元素。比如说,我们定义年龄越小优先级越高,以下是插入过程:

编写层次聚类算法 - 图3

取出的第一个元素是Yui,因为她的年龄最小:

编写层次聚类算法 - 图4

我们看看Python中如何使用优先队列:

  1. >>> from Queue import PriorityQueue # 加载优先队列类
  2. >>> singersQueue = PriorityQueue() # 创建对象
  3. >>> singersQueue.put((16, 'Suzuka Nakamoto')) # 插入元素
  4. >>> singersQueue.put((15, 'Moa Kikuchi'))
  5. >>> singersQueue.put((14, 'Yui Mizuno'))
  6. >>> singersQueue.put((17, 'Ayaka Sasaki'))
  7. >>> singersQueue.get() # 获取第一个元素,即最年轻的歌手Yui。
  8. (14, 'Yui Mizuno')
  9. >>> singersQueue.get()
  10. (15, 'Moa Kikuchi')
  11. >>> singersQueue.get()
  12. (16, 'Suzuka Nakamoto')
  13. >>> singersQueue.get()
  14. (17, 'Ayaka Sasaki')

在进行聚类时,我们将分类、离它最近的分类、以及距离插入到优先队列中,距离作为优先级。比如上面的犬种示例,Border Collie最近的分类是Portuguese WD,距离是0.232:

编写层次聚类算法 - 图5

我们将优先队列中距离最小的两个分类取出来,合并成一个分类,并重新插入到优先队列中。比如下图是将Border Collie和Portuguese WD合并后的结果:

编写层次聚类算法 - 图6

重复这个过程,直到队列中只有一个元素为止。当然,我们插入的数据会复杂一些,请看下面的讲解。

从文件中读取数据

数据文件是CSV格式的(以逗号分隔),第一行是列名,第一列是犬种,第二列之后是特征值:

编写层次聚类算法 - 图7

我们用Python的列表结构来存储这些数据,data[0]用来存放所有记录的分类,如data[0][0]是Border Collie,data[0][1]是Boston Terrier。data[1]则是所有记录的高度,data[2]是重量。

特征列的数据都会转换成浮点类型,如data[1][0]是20.0,data[2][0]是45.0等。在读取数据时就需要对其进行标准化。此外,我们接下来会使用“下标”这个术语,如第一条记录Border Collie的下标是0,第二条记录Boston Terrier下标是1等。

初始化优先队列

以Border Collie为例,我们需要计算它和其它犬种的距离,保存在Python字典里:

  1. {1: ((0, 1), 1.0244), # Border Collie(下标为0)和Boston Terrier(下标为1)之间的距离为1.0244
  2. 2: ((0, 2), 0.463), # Border Collie和Brittany Spaniel(下标为2)之间的距离为0.463
  3. ...
  4. 10: ((0, 10), 2.756)} # Border Collie和Yorkshire Terrier的距离为2.756

此外,我们会记录Border Collie最近的分类及距离:这对犬种是(0, 8),即下标为0的Border Collie和下标为8的Portuguese WD,距离是0.232。

距离相等的问题以及为何要使用元组

你也许注意到了,Portuguese WD和Standard Poodle的距离是0.566,Boston Terrier和Brittany Spaniel的距离也是0.566,

如果我们通过最短距离来取,很可能会取出Standard Poodle和Boston Terrier进行组合,这显然是错误的,所以我们才会使用元组来存放这对犬种的下标,以作判断。比如说,Portuguese WD的记录是:

  1. ['Portuguese Water Dog', 0.566, (8, 9)]

它的近邻Standard Poodle的记录是:

  1. ['Standard Poodle', 0.566, (8, 9)]

我们可以通过这个元组来判断这两条记录是否是一对。

距离相等的另一个问题

在介绍优先队列时,我用了歌手的年龄举例,如果他们的年龄相等,取出的顺序又是怎样的呢?

编写层次聚类算法 - 图8

可以看到,如果年龄相等,优先队列会根据记录中的第二个元素进行判断,即歌手的姓名,并按字母顺序返回,如Avaka会比Moa优先返回。

在犬种示例中,我们让距离成为第一优先级,下标成为第二优先级。因此,我们插入到优先队列的一条完整记录是这样的:

编写层次聚类算法 - 图9

重复下述步骤,直到仅剩一个分类

我们从优先队列中取出两个元素,对它们进行合并。如合并Border Collie和Portuguese WD后,会形成一个新的分类:

  1. ['Border Collie', 'Portuguese WD']

然后我们需要计算新的分类和其它分类之间的距离,方法是对取出的两个分类的距离字典进行合并。如第一个分类的距离字段是distanceDict1,第二个分类的是distanceDict2,新的距离字段是newDistanceDict:

  1. 初始化newDistanceDict
  2. 对于distanceDict1的每一个键值对:
  3. 如果这个键在distanceDict2中存在:
  4. 如果这个键在distanceDict1中的距离要比在distanceDict2中的距离小:
  5. distanceDict1中的距离存入newDistanceDict
  6. 否则:
  7. distanceDict2中的距离存入newDistanceDict

编写层次聚类算法 - 图10

经过计算后,插入到优先队列中的新分类的完整记录是:

编写层次聚类算法 - 图11

代码实践

你能将上面的算法用Python实现吗?你可以从hierarchicalClustererTemplate.py这个文件开始,完成以下步骤:

  1. 编写init方法,对于每条记录:
    1. 计算该分类和其它分类之间的欧几里得距离;
    2. 找出该分类的近邻;
    3. 将这些信息放到优先队列的中。
  2. 编写cluster方法,重复以下步骤,直至剩下一个分类:
    1. 从优先队列中获取两个元素;
    2. 合并;
    3. 将合并后的分类放回优先队列中。

编写层次聚类算法 - 图12

解答

注意,我的实现并不一定是最好的,你可以写出更好的!

  1. from queue import PriorityQueue
  2. import math
  3. """
  4. 层次聚类示例代码
  5. """
  6. def getMedian(alist):
  7. """计算中位数"""
  8. tmp = list(alist)
  9. tmp.sort()
  10. alen = len(tmp)
  11. if (alen % 2) == 1:
  12. return tmp[alen // 2]
  13. else:
  14. return (tmp[alen // 2] + tmp[(alen // 2) - 1]) / 2
  15. def normalizeColumn(column):
  16. """计算修正的标准分"""
  17. median = getMedian(column)
  18. asd = sum([abs(x - median) for x in column]) / len(column)
  19. result = [(x - median) / asd for x in column]
  20. return result
  21. class hClusterer:
  22. """该聚类器默认数据的第一列是标签,其它列是数值型的特征。"""
  23. def __init__(self, filename):
  24. file = open(filename)
  25. self.data = {}
  26. self.counter = 0
  27. self.queue = PriorityQueue()
  28. lines = file.readlines()
  29. file.close()
  30. header = lines[0].split(',')
  31. self.cols = len(header)
  32. self.data = [[] for i in range(len(header))]
  33. for line in lines[1:]:
  34. cells = line.split(',')
  35. toggle = 0
  36. for cell in range(self.cols):
  37. if toggle == 0:
  38. self.data[cell].append(cells[cell])
  39. toggle = 1
  40. else:
  41. self.data[cell].append(float(cells[cell]))
  42. # 标准化特征列(即跳过第一列)
  43. for i in range(1, self.cols):
  44. self.data[i] = normalizeColumn(self.data[i])
  45. ###
  46. ### 数据已经读入内存并做了标准化,对于每一条记录,将执行以下步骤:
  47. ### 1. 计算该分类和其他分类的距离,如当前分类的下标是1,
  48. ### 它和下标为2及下标为3的分类之间的距离用以下形式表示:
  49. ### {2: ((1, 2), 1.23), 3: ((1, 3), 2.3)... }
  50. ### 2. 找出距离最近的分类;
  51. ### 3. 将该分类插入到优先队列中。
  52. ###
  53. # 插入队列
  54. rows = len(self.data[0])
  55. for i in range(rows):
  56. minDistance = 99999
  57. nearestNeighbor = 0
  58. neighbors = {}
  59. for j in range(rows):
  60. if i != j:
  61. dist = self.distance(i, j)
  62. if i < j:
  63. pair = (i,j)
  64. else:
  65. pair = (j,i)
  66. neighbors[j] = (pair, dist)
  67. if dist < minDistance:
  68. minDistance = dist
  69. nearestNeighbor = j
  70. nearestNum = j
  71. # 记录这两个分类的配对信息
  72. if i < nearestNeighbor:
  73. nearestPair = (i, nearestNeighbor)
  74. else:
  75. nearestPair = (nearestNeighbor, i)
  76. # 插入优先队列
  77. self.queue.put((minDistance, self.counter,
  78. [[self.data[0][i]], nearestPair, neighbors]))
  79. self.counter += 1
  80. def distance(self, i, j):
  81. sumSquares = 0
  82. for k in range(1, self.cols):
  83. sumSquares += (self.data[k][i] - self.data[k][j])**2
  84. return math.sqrt(sumSquares)
  85. def cluster(self):
  86. done = False
  87. while not done:
  88. topOne = self.queue.get()
  89. nearestPair = topOne[2][1]
  90. if not self.queue.empty():
  91. nextOne = self.queue.get()
  92. nearPair = nextOne[2][1]
  93. tmp = []
  94. ## 我从队列中取出了两个元素:topOne和nextOne,
  95. ## 检查这两个分类是否是一对,如果不是就继续从优先队列中取出元素,
  96. ## 直至找到topOne的配对分类为止。
  97. while nearPair != nearestPair:
  98. tmp.append((nextOne[0], self.counter, nextOne[2]))
  99. self.counter += 1
  100. nextOne = self.queue.get()
  101. nearPair = nextOne[2][1]
  102. ## 将不处理的元素退回给优先队列
  103. for item in tmp:
  104. self.queue.put(item)
  105. if len(topOne[2][0]) == 1:
  106. item1 = topOne[2][0][0]
  107. else:
  108. item1 = topOne[2][0]
  109. if len(nextOne[2][0]) == 1:
  110. item2 = nextOne[2][0][0]
  111. else:
  112. item2 = nextOne[2][0]
  113. ## curCluster即合并后的分类
  114. curCluster = (item1, item2)
  115. ## 对于这个新的分类需要做两件事情:首先找到离它最近的分类,然后合并距离字典。
  116. ## 如果item1和元素23的距离是2,item2和元素23的距离是4,我们取较小的那个距离,即单链聚类。
  117. minDistance = 99999
  118. nearestPair = ()
  119. nearestNeighbor = ''
  120. merged = {}
  121. nNeighbors = nextOne[2][2]
  122. for (key, value) in topOne[2][2].items():
  123. if key in nNeighbors:
  124. if nNeighbors[key][1] < value[1]:
  125. dist = nNeighbors[key]
  126. else:
  127. dist = value
  128. if dist[1] < minDistance:
  129. minDistance = dist[1]
  130. nearestPair = dist[0]
  131. nearestNeighbor = key
  132. merged[key] = dist
  133. if merged == {}:
  134. return curCluster
  135. else:
  136. self.queue.put( (minDistance, self.counter,
  137. [curCluster, nearestPair, merged]))
  138. self.counter += 1
  139. def printDendrogram(T, sep=3):
  140. """打印二叉树状图。树的每个节点是一个二元组。这个方法摘自:
  141. http://code.activestate.com/recipes/139422-dendrogram-drawing/"""
  142. def isPair(T):
  143. return type(T) == tuple and len(T) == 2
  144. def maxHeight(T):
  145. if isPair(T):
  146. h = max(maxHeight(T[0]), maxHeight(T[1]))
  147. else:
  148. h = len(str(T))
  149. return h + sep
  150. activeLevels = {}
  151. def traverse(T, h, isFirst):
  152. if isPair(T):
  153. traverse(T[0], h-sep, 1)
  154. s = [' ']*(h-sep)
  155. s.append('|')
  156. else:
  157. s = list(str(T))
  158. s.append(' ')
  159. while len(s) < h:
  160. s.append('-')
  161. if (isFirst >= 0):
  162. s.append('+')
  163. if isFirst:
  164. activeLevels[h] = 1
  165. else:
  166. del activeLevels[h]
  167. A = list(activeLevels)
  168. A.sort()
  169. for L in A:
  170. if len(s) < L:
  171. while len(s) < L:
  172. s.append(' ')
  173. s.append('|')
  174. print (''.join(s))
  175. if isPair(T):
  176. traverse(T[1], h-sep, 0)
  177. traverse(T, maxHeight(T), -1)
  178. filename = '/Users/raz/Dropbox/guide/data/dogs.csv'
  179. hg = hClusterer(filename)
  180. cluster = hg.cluster()
  181. printDendrogram(cluster)

运行结果和我们手算的一致:

编写层次聚类算法 - 图13

动手实践

这里提供了77种早餐麦片的营养信息,包括以下几项:

  • 麦片名称
  • 热量
  • 蛋白质
  • 脂肪
  • 纤维
  • 碳水化合物
  • 维生素

编写层次聚类算法 - 图14

请对这个数据集进行层次聚类:

  • 哪种麦片和Trix最相近?
  • 与Muesli Raisins & Almonds最相近的是?

数据集来自:http://lib.stat.cmu.edu/DASL/Datafiles/Cereals.html

结果

我们只需将代码中的文件名替换掉就可以了,结果如下:

编写层次聚类算法 - 图15

因此Trix和Fruity Pebbles最相似(你可以去买这两种麦片尝尝)。Muesli Raisins & Almonds和Muesli Peaches & Pecans最相似。

编写层次聚类算法 - 图16

好了,这就是层次聚类算法,很简单吧!