cs Python ML 集体智慧编程

No.2 聚类-发现群组 Discovering Groups

Programming Collective Intelligence No.2

Posted by kamzero on 2019-08-28

Programming Collective Intelligence No.2 聚类-发现群组

前言

前天花了大半天时间实现了推荐系统之后，昨天主要在看书、学英语，虽然纯看书效率不高，但是保持一天输入一天输出的节奏还是相对比较合理的。既然有了想做的事，就保持专注吧。

昨天浏览了大半本《集体智慧编程》，对书中所要讲述的内容已经有了全面的了解。第三章是讲聚类，第四章搭建搜索引擎，都是比较有趣的。这不单单是一本讲机器学习的书，正如副标题“构建Web2.0应用”所言，这是一本教大家如何利用数据挖掘出信息、创建更为成功的网站的书。

虽然我很想很想上手搜索引擎，很想一步跳到后面几章的决策树、非负矩阵因式分解的实践上去，但是以我现在对Python的熟悉程度和代码能力，还是乖巧一点，先把聚类搞定了再说叭。开始~

今天要实现什么？

  • [x] 从各种不同的数据来源中构造算法所需要的数据

  • [x] 实现两种不同的聚类算法

    • [x] 分级聚类
    • [x] K-均值聚类

  • [x] 更多有关距离度量的知识

  • [x] 简单的图形可视化

  • [x] 将复杂的数据投影到二维空间中

概述

聚类也适用于数据量比较大的时候，顾名思义，就是按照特性把相近的对象聚在一起分成若干类。比如对于商家而言，年龄、收入相仿的人可能有完全不一样的穿衣风格，但是我们可以通过聚类算法分出若干种时装风格，以此为依据打造商品。令人惊叹的是，聚类在计量生物学中也有应用，科学家们用它来寻找具有相似行为的基因组，这些基因组中的基因会以相同的方式相应外界活动。

监督学习和无监督学习

第一次听到这两个概念是在去年中控杯的培训上，形象地讲，监督学习就是有老师指导的学习，我们知道正确答案，我们知道期望，我们利用样本和期望输出来学习如和预测。常见的监督学习法有：神经网络、决策树、向量支持机、贝叶斯过滤。而无监督学习就是，我们自己也不知道正确答案，但知道怎样的结果是更优的，比如AlphaGo和他的弟弟后期练习围棋时，就是两个机器对弈，没有老师指导，但我们知道赢了棋的更优，这就是无监督学习。更为精准的表述是，无监督学习并不是利用带有正确答案的样本数据进行训练，它的目的是在一组数据中找寻某种结构，而数据本身并不是我们我们要找的答案。

在时装的例子中，聚类的结果并不会告诉商家每一位顾客可能买什么，也不会预测新的顾客适合哪种时装，而是从采集的数据中发现不同的群组。

任务一：博客聚类分组

准备-构造数据

我们想做的第一件事是对博客中的词汇出现频度进行分析，可以根据用户撰写的文字对博客用户进行分组，也可以根据内容对一个博客下的文章进行分组，还可以根据大家对词汇的用法将词汇分组。

在这个问题上，英文的书写模式有着很大的优势，多为一词一意，便于对语句进行分割、分析。而中文词与词之间没有明显的间隔，当然可以按单个汉字来分析，但若想要进行词汇的分析就较为复杂。我在网上了解到了几个优秀的中文词分割算法，似乎以贪心为主，留下待之后实践。https://blog.csdn.net/dancegreen/article/details/79953387

单词向量

为聚类算法准备数据的常见做法是定义一组公共的数值型属性。**我们可以利用这些属性对数据项进行比较，比如第二章用1-5的数来代表对影片的评分情况，在这里，我们可以使用单词出现的次数，也可以归一化之后使用单词出现的频率。为此，我们需要下载一系列博客订阅源，从中提取文本，并建立起一个单词频度表。

对订阅源中的单词进行计数

几乎所有的博客都可以在线阅读，或者通过 RSS 订阅源进行阅读。RSS订阅源是个包含博客及其所有文章条目信息的简单的XML文档。为了给每个博客中的单词计数，首先要解析这些订阅源。有一个非常不错的库能够完成这项工作，它就是feedparser，我们可以从https://pypi.org/project/feedparser/上下载到它。

我们可以用这个库轻松地从任何RSS或Atom订阅源中得到标题、链接和文章条目，这是一个从订阅源中提取所有单词的函数。

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
import feedparser
import re
#re是正则表达式的缩写

#返回一个RSS订阅源的标题和包含单词计数情况的字典
def getwordcounts(url):
    #根据url解析订阅源
    d = feedparser.parse(url)
    wc={}

    #循环遍历所有的文章条目。可大概猜测出每个函数和变量的意义。
    for e in d.entries:
        if 'summary' in e: summary = e.summary
        else: summary=e.description

        #提取单词列表
        words=getwords(e.title+' '+summary)
        for word in words:
            wc.setdefault(word,0)
            wc[word]+=1
    return d.feed.title,wc

def getwords(html):
    #去除所有的html标记
    txt=re.compile(r'<[^>]+>').sub('',html)

    #利用所有非字母字符拆分出单词
    #此处存疑,不能拆分出如C++这样的词
    words=re.compile(r'[^A-Z^a-z]+').split(txt)

    #转化为小写
    return [word.lower() for word in words if word!='']

主体代码，开始吧

现在，我们需要一个订阅源的列表，在工作目录下，名为feedlist.txt，里面多是一些文字为主的热门博客。将下列代码加入generatefeedvector.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
#遍历文件中每一行，生成针对每个博客的单词统计，以及出现这些单词的博客数目apcount
apcount={}
wordcounts={}
feedlist=[line for line in open('feedlist.txt')]
for feedurl in feedlist:
    title,wc=getwordcounts(feedurl)
    wordcounts[title]=wc
    for word,count in wc.items():
        apcount.setdefault(word,0)
        if count>1:
            apcount[word]+=1

#只保留出现频率在每个百分比范围内的单词，减少需要考察的单词总量，我们认为这些单词是有特色的
wordlist=[]
for w,bc in apcount.items():
    frac=float(bc)/len(feedlist)
    if frac>0.1 and frac<0.5: wordlist.append(w)

#生成文件，并输出一个大矩阵
out=open('blogdata.txt','w')
out.write('Blog')
for word in wordlist: out.write('\t%s' % word)
out.write('\n')
for blog,wc in wordcounts.items():
    out.write(blog)
    for word in wordlist:
        if word in wc: out.write('\t%s' % wc[word])
        else: out.write('\t0')
    out.write('\n')

运行花了很长时间，少说得有十几分钟，在我背了两个list的单词之后，终于生成了blogdata.txt!

分级聚类

分级积累通过连续不断地将最为相似的群组两两合并，来构造出一个群组的层级结构。在每次迭代的过程中，分级聚类算法会计算每两个群组间的距离，并将距离最近的两个群组合并成一个新的群组。不知道为什么，我很自然地联系起了离散数学中学到的霍夫曼编码。

对于分级聚类有一种可视化形式，就是树状图。树状图不仅可以利用连线来表达每个聚类的构成情况，而且还可以利用距离来体现构成聚类的各元素间相隔的远近，这能帮助我们有效地确定一个聚类中各元素间的相似程度，并以此来指示聚类的紧密程度。

加载文件

首先，我们需要加载数据文件。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
def readfile(filename):
    #逐行读取，这里示例代码有误，已更正
    lines=[line for line in open(filename)]

    #第一行是列标题
    colnames=lines[0].strip().split('\t')[1:]
    rownames=[]
    data=[]
    for line in lines[1:]:
        p=line.strip().split('\t')
        #每行的第一列是行名
        rownames.append(p[0])
        #剩余部分就是该行对应的数据
        data.append([float(x) for x in p[1:]]) #给最后加一个list
    return rownames,colnames,data

相关度

下一步，我们来定义紧密度。为了纠正文章长度对结果的影响，我们采用皮尔逊相关度，它判断两组数据与某条直线的拟合程度。我们以两个数字列表为参数，返回他们的皮尔逊相关度。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
from math import sqrt
def pearson(v1,v2):
    #简单求和
    sum1=sum(v1)
    sum2=sum(v2)

    #求平方和
    sum1Sq=sum([pow(v,2) for v in v1])
    sum2Sq=sum([pow(v,2) for v in v2])

    #求乘积之和
    pSum=sum([v1[i]*v2[i] for i in range(len(v1))])

    #计算r
    num = pSum-(sum1*sum2/len(v1))
    den = sqrt((sum1Sq-pow(sum1,2)/len(v1))*(sum2Sq-pow(sum2,2)/len(v1)))
    if den==0 : return 0

    return 1.0-num/den

聚类之类

我们要创建一个聚类的类，把一个聚类具有的属性放入其中，并以此来描述这颗层级树。每一个聚类，可以是树中的枝节点，也可以是与数据集中实际数据行相对应的叶节点。除此之外，每一个聚类还包括了指示其位置的信息，这个信息可以来自叶节点的行数据，也可以来自枝节点经合并后的数据。

1
2
3
4
5
6
7
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

聚类合并

分级聚类算法以一组对应于原始数据项的聚类开始。函数的主循环部分会尝试每一组可能的配对并计算他们的相关度，以此来找出最佳配对。最佳配对的两个聚类会被合并成一个聚类，新聚类中包含的各项数据，等于将两个旧聚类中的数据求均值之后得到的结果。这一过程一直重复，直到只剩下一个聚类。

由于整个计算过程非常耗时，所以不妨将每个配对的相关度计算结果保存起来，因为这样的计算会反复发生，直到配对中的某一项被合并到一个新的聚类中为止。

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
def hcluster(rows,distance=pearson):
    distances={}
    currentclustid=-1

    #最初的聚类是数据集中的行
    clust=[bicluster(rows[i],id=i) for i in range(len(rows))]

    while len(clust)>1:
        lowestpair=(0,1)
        closest=distance(clust[0].vec,clust[1].vec)

        #遍历每一个配对，寻找最小值
        #用distances来缓存距离的计算值
        for i in range(len(clust)):
            for j in range(i+1,len(clust)):
                if (clust[i].id,clust[j].id) not in distances:
                    distances[(clust[i].id,clust[j].id)]=distance(clust[i].vec,clust[j].vec)

                d=distances[(clust[i].id,clust[j].id)]
                
                if d<closest:
                    closest=d
                    lowestpair=(i,j)

        #计算两个聚类的平均值
        mergevec=[
            (clust[lowestpair[0]].vec[i]+clust[lowestpair[1]].vec[i])/2.0
            for i in range(len(clust[0].vec))
        ]

        #建立新的聚类 
        newcluster=bicluster(
            mergevec,left=clust[lowestpair[0]],right=clust[lowestpair[1]],
            distance=closest,id=currentclustid
        )

        #不在原始集合中的聚类，id为负数？
        currentclustid-=1
        del clust[lowestpair[0]]
        del clust[lowestpair[1]]
        clust.append(newcluster)

    return clust[0]

最最最简单的可视化

其实就是递归遍历聚类树，以类似文件系统结构的形式把它打印出来啦。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
def printclust(clust,labels=None,n=0):
    #利用缩进来建立层级布局
    for i in range(n): print(' ')
    
    if clust.id<0:
        #负数标记 - 分支
        print('-')
    else:
        #正数标记 - 叶节点
        if labels==None:print(clust.id)
        else: print(labels[clust.id])

    if clust.left!=None: printclust(clust.left,labels=labels,n=n+1)
    if clust.right!=None: printclust(clust.right,labels=labels,n=n+1)

打开一个python会话，输入以下命令

1
2
3
4
import clusters
blognames,words,data=clusters.readfile('blogdata.txt')
clust=clusters.hcluster(data)
clusters.printclust(clust,labels=blognames)

绘制树状图

由于刚才模仿文件层级输出的可视化效果不太好，我决定采用绘制树状图的方法。这里用到了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
from PIL import Image,ImageDraw

#递归确定聚类的总体高度
def getheight(clust):
    #叶节点 高度为1
    if clust.left==None and clust.right==None: return 1

    #否则高度为两分支之和
    return getheight(clust.left)+getheight(clust.right)

#求根节点的总体误差，以此为依据调整线条长度
def getdepth(clust):
    #叶节点距离为0
    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
        #线的长度
        ll=clust.distance*scaling
        #聚类到其子节点的垂直线
        draw.line((x,top+h1/2,x,bottom-h2/2),(255,0,0))
        #连接左侧节点的水平线
        draw.line((x,top+h1/2,x+ll,top+h1/2),(255,0,0))
        #连接右侧节点的水平线
        draw.line((x,bottom-h2/2,x+ll,bottom-h2/2),(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))
1
clusters.drawdendrogram(clust,blognames,jpeg='blogclust.jpg')

1566914526883

列聚类

就像第二章基于用户的推荐和基于物品的推荐一样，在这个聚类中，我们只要把单词和博客反转，把大矩阵转置，就可以对单词进行聚类，这样可以知晓哪些单词经常结合在一起使用。对于商家来说，他们可以知道哪些商品应该捆绑销售。

1
2
3
4
5
6
7
#矩阵转置
def rotatematrix(data):
    newdata=[]
    for i in range(len(data[0])):
        newrow=[data[j][i] for j in range(len(data))]
        newdata.append(newrow)
    return newdata
1
2
3
4
reload(clusters)
rdata=clusters.rotatematrix(data)
wordclust=clusters.hcluster(rdata)
clusters.drawdendrogram(wordclust,labels=words,jpeg='wordclust.jpg')

此处如预期般生成了一张长图，我没有调整合适的单词频范围，数据项又很多，因此无意义的聚类比较多，不如对博客聚类合理。

对同一博客文章的聚类

此处留白，在之前的基础上这里很好实现，只需要稍微更改构建数据集的部分，而且我觉得这比较有意义，因为我也面临给博文加tag分类的问题，这也许是一个有效的方案。

K-均值聚类

分级聚类有两个缺点：一是结果是一棵树，并不会自动拆成分组；二是计算量非常大。另一种可供选择的聚类方法是K-均值聚类，它完全不同于分级聚类，是在已知分组数K的情况下根据数据的结构确定聚类的大小。

算法细则

  1. 随机确定k个中心位置
  2. 将每个数据项分配给最邻近的中心点
  3. 聚类中心点移动到分配给该中心点的所有节点的平均位置处
  4. 重复这一过程，直到分配过程不再产生变化

一些疑问

这里很自然地生出了疑问，显然，聚类结果是依赖于初始随机点位置的，比如：一个中心点离所有数据项都不是最近的，那么聚类就会减少，出现空聚类。那么，随机确定K个中心位置真的可以吗？有没有改进方法？

代码实现

K-均值聚类算法思想非常容易理解，代码实现相比分级聚类还是略麻烦一些的。

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
import random

def kcluster(rows,distance=pearson,k=4):
    #确定每个点的最小值和最大值？
    ranges=[(min([row[i] for row in rows]),max([row[i] for row in rows])) for i in range(len(rows[0]))]

    #随机创立K个中心点
    clusters=[
        [random.random()*(ranges[i][1]-ranges[i][0])+ranges[i][0]
        for i in range(len(rows[0]))]for j in range (k)
    ]

    lastmatches = None
    for t in range(100):
        print('Iteration %d'%t) #第t次迭代
        bestmatches=[[]for i in range(k)]

        #在每一行中寻找最近的中心点
        for j in range(len(rows)):
            row=rows[j]
            bestmatch=0 #初始假设第一个中心点距离最近
            for i in range(k):
                d=distance(clusters[i],row) 
                if d<distance(clusters[bestmatch],row):bestmatch=i
            bestmatches[bestmatch].append(j)

        #如果结果与上一次相同，迭代结束
        if bestmatches==lastmatches:break
        lastmatches=bestmatches

        #把中心点移到其所有成员的平均位置处
        for i in range(k):
            avgs=[0.0]*len(rows[0])
            if len(bestmatches[i])>0:
                for rowid in bestmatches[i]:
                    for m in range(len(rows[rowid])):
                        avgs[m]+=rows[rowid][m]
                for j in range(len(avgs)):
                    avgs[j]/=len(bestmatches[i])
                clusters[i]=avgs

    return bestmatches
1
2
kclust=clusters.kcluster(data,k=10)
[blognames[r] for r in kclust[0]]

输出

[‘Unknown title’, ‘NB Blog Feed’, ‘ThinkProgress’, ‘TechEBlog’, ‘TMZ.com’]

多维缩放计数

在真实生活的例子中，我们所要聚类的内容不只包含两个数值，因此不能以二维的形式绘制出来。但是，为了弄明白不同物品之间的关系，把他们绘制在一张纸上，并且用距离的远近来表达相似程度，是一种非常有效的做法。

算法细则

  1. 计算出所有数据项两两之间的目标距离
  2. 将所有数据项随机放置在二维图上
  3. 针对每两两构成的数据项，我们将他们的目标距离与当前距离进行比较，求出误差值，根据误差的情况，按比例将每个数据项的位置移动一定距离
  4. 每一个节点的移动，都是其他节点施加在该节点上推或拉的综合效应。重复该过程，直到我们再也无法通过移动节点来减少总体误差为止。

代码实现

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
def scaledown(data,distance=pearson,rate=0.01):
    n=len(data)

    #每一对数据项之间的距离
    realdist=[
        [distance(data[i],data[j]) for j in range(n)] for i in range(0,n)
    ]
    outersum=0.0

    #随机初始化节点在二维空间中的起始位置
    loc=[[random.random(),random.random()] for i in range(n)]
    fakedist=[[0.0 for j in range(n)] for i in range(n)]

    lasterror=None
    for m in range(0,1000):
        #寻找投影后的距离
        for i in range(n):
            for j in range(n):
                fakedist[i][j]=sqrt(
                    sum([pow(loc[i][x]-loc[j][x],2) for x in range(len(loc[i]))
                    ]))

        #移动节点
        grad=[[0.0,0.0] for i in range(n)]

        totalerror = 0 
        for k in range(n):
            for j in range(n):
                if j==k: continue
                #误差值等于目标距离与当前距离之间差值的百分比
                errorterm=(fakedist[j][k]-realdist[j][k])/realdist[j][k]

                #每一个节点根据误差按比例移动
                grad[k][0]+=((loc[k][0]-loc[j][0])/fakedist[j][k])*errorterm
                grad[k][1]+=((loc[k][1]-loc[j][1])/fakedist[j][k])*errorterm

                #记录总的误差值
                totalerror+=abs(errorterm)
        
        print(totalerror)

        #如果节点移动之后误差更大，结束
        if lasterror and lasterror<totalerror:break
        lasterror=totalerror

        #根据rate与grad移动每一个节点
        for k in range(n):
            loc[k][0]-=rate*grad[k][0]
            loc[k][1]-=rate*grad[k][1]

    return loc

可视化-生成二维图

1
2
3
4
5
6
7
8
def draw2d(data,labels,jpeg='mds2d.jpg'):
    img=Image.new('RGB',(2000,2000),(255,255,255))
    draw=ImageDraw.Draw(img)
    for i in range(len(data)):
        x=(data[i][0]+0.5)*1000
        y=(data[i][1]+0.5)*1000
    draw.text((x,y),labels[i],(0,0,0))
    img.save(jpeg,'JPEG')
1
2
3
4
5
from imp import reload
reload(clusters)
blognames,words,data=clusters.readfile('blogdata.txt')
coords=clusters.scaledown(data)
clusters.draw2d(coords,blognames,jpeg='blogs2d.jpg')

1566962538896

遇到的bug

  • 被零除的问题，特判一下
  • 只画了一个点出来，考虑画图函数问题，draw被放在了循环外面

后记

那些不经意间学到的东西，待整理

  • feedparser 的用法

  • 字典与setdefault

  • python2/3 file与open

  • python3 逐行读取文件https://blog.csdn.net/zhengxiangwen/article/details/55148287

  • AttributeError: object has no attribute ‘title’

    • The URL you are trying to parse with feedparser is either not a valid feed (check it with feedvalidator), but a web page, or the feed is empty, or the title is empty.

      As a workaround, use getattr():

      1
      return getattr(d.feed, 'title', 'Unknown title'), wc

  • s.strip(rm) 删除s字符串中开头、结尾处，位于 rm删除序列的字符 当rm为空时，默认删除空白符（包括’\n’, ‘\r’, ‘\t’, ’ ')

  • str_split ``= str``.split(``'||'``,``2``)按照||进行分割，并且只分割前两个

  • beautiful soup

  • reload前import sys

可以做的有趣的事

  • 可以用beautiful soup搞一些事情

  • 昨天和cby交谈时有了图助内部书籍借阅网站的想法，web2.0的很多东西可以在图助网站上实践

  • 以后再想读这种书做这种整理的时候，把英文翻译成中文

  • 可以用它来分析很多有趣的事情，比如足球数据…😂

写在后面

第二次写学习总结长博客了，慢慢有了一种总结、梳理、提出问题、解决问题的感觉。第一次大多是抄书，很多代码只是复制样例，然后改一改，这一次基本自己敲过一边。这一次梳理的时间比想象中长很多，几乎写了一天半，昨天上午就开始动手实践了，下午和晚上也都花了些时间整理，结果睡觉之前只完成了分级聚类。今天上午终于写完了，也实现了一点事情，可以休息啦。

FEATURED TAGS
cs Python ML 集体智慧编程
FRIENDS