第5章聚类问题
本章学习要点
 了解聚类问题的基本概念、概率模型以及与分类问题的异同点。
 了解并掌握含隐变量概率模型及其EM求解算法。
 了解并掌握高斯混合模型、三硬币模型这两种常用含隐变量的概率模型。
 了解k均值聚类、DBSCAN聚类这两种常用的聚类方法。

聚类(clustering)是一种特殊的分类问题。分类是根据有标注数据集来训练模型，学习人们预先设定的类别概念，属于有监督学习； 而聚类则是根据无标注
（仅包含特征）
数据集训练模型，即根据数据自身的分布特性或结构，自动将数据聚集成簇(cluster)，形成类别概念，它属于无监督学习。
聚类问题可形式化描述为： 给定无标注数据集D={x1,x2,…,xm}，其中xi是d维样本特征，m为样本容量，然后设计聚类算法训练模型，将数据集D划分成k个不相交的簇{C1,C2,…,Ck}，其中

D=∪ki=1Ci，且若i≠j，则Ci∩Cj=Φ

今后任给新的样本特征x，也可以通过聚类模型将其划归某一类簇。

5.1聚类问题的提出
在分类问题中，如果因为标注训练集的工作量过大，或问题本身并没有预设的类别概念，则分类问题就演变成了聚类问题。为便于后续讲解，这里先对概率符号做一下简化，将离散型概率分布P(X=x)或连续型概率密度p(x)统一记作P(x)，例如

P(X=x)≡P(x)，P(X=x,Y)≡P(x,Y)，P(Y|X=x)≡P(Y|x)

P(Y=y)≡P(y)，P(X,Y=y)≡P(X,y)，P(Y=y|X)≡P(y|X)

P（X）=∫ΩYp(X,y； θ)dy≡∑YP(X,Y； θ)≡∑y∈ΩYP(X,y； θ)

请读者根据上下文加以理解。下面分别给出分类问题、聚类问题的形式化描述。
5.1.1分类问题概述
在分类问题中，假设有n个类别{c1,c2,…,cn}，将分类特征记作随机变量X，类别记作随机变量Y。类别Y是离散型随机变量，其值域为ΩY={c1,c2,…,cn}，概率分布(先验概率)为多项分布，可记作P(Y； α)，即

P(Y=ck)=αk，记作P(Yk； αk)，k=1,2,…,n(51)





其中α=(α1,α2,…,αn)为未知参数，且∑nk=1αk=1。再假设各类特征具有相同的概率分布形式，但所取参数βk不同，将它们记作

P(X|Yk； βk)，k=1,2,…,n(52)

其中，β=(β1,β2,…,βn)为未知参数。
可以根据类别的先验概率分布P(Y； α)和各类的特征概率分布P(X|Yk； βk)设计贝叶斯分类器。给定特征X，将类别Y判定为后验概率P(Yk|X)最大的ck，其分类判别函数为

k*= argmaxk=1,2,…,nP(Yk|X)(53)

或其等价形式

k*= argmaxk=1,2,…,nP(X,Yk)= argmaxk=1,2,…,nP(X|Yk)P(Yk)(54)

为明确未知参数，可在式（54）中标出α、β，即

k*= argmaxk=1,2,…,nP(X,Yk； αk,βk)= argmaxk=1,2,…,nP(X|Yk； βk)P(Yk； αk)(55)

可以看出，分类问题的关键是如何确定未知参数α、β。
给定数据集Dtrain={(x1,y1),(x2,y2),…,(xm,ym)}，其样本容量为m，可以使用极大似然估计方法分别估计出参数α和β。对于参数α，若Dtrain中类别取值为ck的样本点个数等于mk，则使用极大似然估计方法即可估计出最优参数α*为

α*k=mkm，P(Yk)=α*k，k=1,2,…,n(56)

对于参数β，记Dtrain中类别取值为ck样本点所构成的子集为Dk。基于子集Dk，使用极大似然估计方法可以估计第k类特征概率分布(离散型)或概率密度函数(连续型)P(X|Yk； βk)中的最优参数β*k(例如正态分布中的均值μk和方差σ2k)，其似然函数为

L(βk)=∏mkj=1P(xj|Yk； βk)，k=1,2,…,n(57)

最大化式（57）的对数似然函数lnL(βk)，得

β*k= argmaxβk ln L(βk)= argmaxβk∑mkj=1ln P(xj|Yk； βk)，k=1,2,…,n(58)

在估计出最优参数α*=(α*1,α*2,…，α*n)、β*=(β*1,β*2,…，β*n)之后，给定任意新样本特征x，可以按式（55）的分类判别函数进行分类，即

k*= argmaxk=1,2,…,nP(x|Yk； β*k)P(Yk； α*k)(59)

总结如下，分类问题就是使用包含类别标注的样本数据集Dtrain={(x1,y1),(x2,y2),…,(xm,ym)}来训练概率模型中的参数α、β； 然后按照最大后验概率原则设计贝叶斯分类器； 今后任给新的样本特征x都可以使用该分类器进行分类。

5.1.2聚类问题概述
设计式（55）的贝叶斯分类器，其关键是如何确定概率模型中的未知参数α、β。如果给定数据集做过标注(通常为人工标注)，即数据集Dtrain={(x1,y1),(x2,y2),…,(xm,ym)}中的类别标注Y已知，则可以使用极大似然估计方法分别估计出最优参数α*和β*。
如果数据集未做标注(通常是因为工作量过大或难以标注)，即数据集D={x1,x2,…,xm}只包含各样本点分类特征X，其对应的类别标注Y未知，称分类特征X是可观测的变量，而类别Y则是不可观测的变量(一般称为隐变量，hidden variable； 或潜变量，latent variable)。如果概率模型含有不可观测的隐变量，这时机器学习该如何确定模型中的未知参数呢？因为数据集未做标注，无法采用式（56）、式（58）的极大似然估计方法来估计最优参数α*、β*。
对于n个类别c1,c2,…,cn的分类问题，特征X的概率分布P(X； α,β)是联合概率分布P(X,Y； α,β)关于特征X的边缘概率，即

P(X； α,β)=∑nk=1P(X|Yk； βk)P(Yk； αk)=∑nk=1P(X,Yk； αk,βk)(510)

其中，α=(α1,α2,…,αn)、β=(β1,β2,…,βn)为未知参数。为了进行更一般化的讨论，记θ=(α,β)，类别Y的值域为ΩY，则式（510）可改写为

P(X； θ)=∑YP(X,Y； θ)≡∑y∈ΩYP(X,y； θ)(511)

给定未做标注的数据集D={x1,x2,…,xm}，可将该数据集的似然函数定义为

L(θ)=∏mj=1P(xj； θ)=∏mj=1∑y∈ΩYP(xj,y； θ)(512)

对式(512)取自然对数，则数据集D的对数似然函数为

ln L(θ)=ln∏mj=1P(xj； θ)=∑mj=1ln∑y∈ΩYP(xj,y； θ)(513)

最大化式(513)的对数似然函数，求解最优参数θ*，即

θ*= argmaxθ ln L(θ)(514)

在确定出最优参数θ*=(α*,β*)之后，即可确定问题的概率模型，并根据式(59)的分类判别函数将数据集D={x1,x2,…,xm}中的m个样本数据划分成不同的类，这就是聚类。
聚类与分类的区别在于： 聚类问题中没有预设类别，它是由机器学习算法按照数据自身(即样本特征)的分布特性自动提炼出来的； 而分类问题中的类别通常是由人工指定的。类别是一种“概念”，是人们在长期的生产、生活实践中总结出来的。进入信息化时代，很多问题通常是先有数据或只有数据，可通过聚类方法对问题进行研究，这就是聚类分析。例如，网上消费产生了很多数据，这时可以通过聚类方法对消费者或消费行为进行研究分析。
在聚类问题中，因为数据集D不能提供各样本点的类别标注，即式(514)的对数似然函数ln L(θ)包含隐变量Y，这是一种含隐变量的最优化问题。如何设计求解含隐变量的最优化算法，这是聚类分析的关键。
总结如下，聚类问题就是使用未做标注的样本数据集D={x1,x2,…,xm}来训练概率模型中的参数α、β，即含隐变量的参数估计； 然后再按最大后验概率原则设计贝叶斯分类器，将样本数据集中的样本数据划分成不同的类。

5.1.3混合概率模型及其参数估计问题
对于n个类别{c1,c2,…,cn}的分类或聚类问题，将分类特征记作随机变量X，类别记作随机变量Y，则X、Y的联合概率分布P(X,Y)即为分类或聚类问题的概率模型。联合概率分布P(X,Y)还可表示为

P(X,Y)≡P(X,Y； α,β)=P(X|Y； β)P(Y； α)(515)

其中，α=(α1,α2,…,αn)、β=(β1,β2,…,βn)为概率模型中的未知参数； P(Y； α)表示类别概率分布，即P(Y=ck)=αk，k=1,2,…,n； P(X|Y； β)表示各类的特征概率分布，它们的分布形式相同但所取参数不同，可记作P(X|Yk； βk)，k=1,2,…,n。
在联合概率分布P(X,Y； α,β)中，边缘概率分布P(X； α,β)可看作分类特征X的概率模型，它可表示为

P(X； α,β)=∑nk=1P(X|Yk； βk)P(Yk； αk)
=α1P(X|Y1； β1)+α2P(X|Y2； β2)+…+αnP(X|Yn； βn)(516)

可以看出，分类特征X的概率模型P(X； α,β)是由n个概率模型P(X|Y1； β1)、P(X|Y2； β2)、……、P(X|Yn； βn)组成的一个混合概率模型，概率α1、α2、……、αn也被称作混合系数，∑nk=1αk=1。
如果已知参数α=(α1,α2,…,αn)、β=(β1,β2,…,βn)，则按照式（516）的混合概率模型产生样本x的过程是： 首先按照概率α1，α2，…，αn选择模型，假设选择第k个模型； 然后按第k个模型的概率分布P(X|Yk； βk)生成样本x。反过来，如果给定样本数据集D={x1,x2,…,xm}，希望估计概率模型P(X； α,β)中的未知参数α、β，这就是混合概率模型的参数估计问题。
混合概率模型的参数估计问题与聚类问题非常相似，其概率模型在本质上是一样的，而且它们用于参数估计的样本数据集D={x1,x2,…,xm}都不包含类别标注。所不同的是，聚类问题在估计出模型参数后，还需进一步设计贝叶斯分类器，将数据集D中的样本数据划分成不同的类。聚类问题、混合概率模型的参数估计问题都属于含隐变量的最优化问题，求解这样的问题通常使用EM算法。

5.2EM算法
EM算法是一种迭代算法，主要用于求解含隐变量的最优化问题。任给初始参数θ0，EM算法的关键步骤是： 第i次迭代时如何将参数从θi-1更新到θi，使得对数似然函数
ln l(θ)=ln P(X； θ)的函数值逐步上升，即ln l(θi)≥ln l(θi-1)，最终收敛至最大值。

5.2.1EM算法原理
1. 问题描述

假设随机变量X、Y服从参数为θ的概率分布P(X,Y； θ)，其中Y为不可观测或未被观测的隐变量。将随机变量Y(例如类别)的值域记作ΩY，则随机变量X(例如分类特征)的边缘概率为

P(X； θ)=∑YP(X,Y； θ)≡∑y∈ΩYP(X,y； θ)

或

P(X； θ)=∑YP(X|Y； θ)P(Y； θ)≡∑y∈ΩYP(X|y； θ)P(y； θ)

给定观测样本X=x，其似然函数l(θ)可定义为

l(θ)=P(x； θ)=∑y∈ΩYP(x,y； θ)

其中，P(x； θ)被称作不完全数据(incompletedata)x的似然函数，P(x,y； θ)被称作完全数据(completedata)(x,y)的似然函数。将观测样本X=x的对数似然函数定义为

ln l(θ)=ln P(x； θ)=ln∑YP(x,Y； θ)≡ln∑y∈ΩYP(x,y； θ)(517)

最大化式(517)的对数似然函数，求解最优参数θ*，即

θ*= argmaxθ ln l(θ)(518)

因为对数似然函数lnl(θ)包含隐变量Y，因此式(518)是一种含隐变量的极大似然估计问题，需通过特殊的EM算法进行求解。

2. 算法准备： Jensen不等式
根据Jensen不等式，设X为随机变量，E(X)为X的期望，E(f(X))为f(X)的期望，则
(1) 对任意凸函数f，有f(E(X))≤E(f(X))。
（2) 对任意凹函数f，有f(E(X))≥E(f(X))。
（3) 如果X为常量，则上述两个不等式均取等号，即如果X=c，则f(X)=f(c)也为常量，因此有f(E(X))=f(c)=E(f(X))。
Jensen不等式有多种不同的表述形式，这里再给出其中的一种： 若函数f为凸集D上的凹函数，则对于任意n个点xj∈D及实数pj(0≤pj≤1)，j=1,2,…,n，且∑nj=1pj=1，都有

f∑nj=1pjxj≥∑nj=1pjf(xj)(519)

如果x1=x2=…=xn=c，则上述不等式取等号。
若将式(519)应用于随机变量X、Y的联合概率分布，则
 当取pj=P(Y|X)，xj=P(X,Y)，f为自然对数(凹函数)时

ln∑YP(Y|X)P(X,Y)≥∑YP(Y|X)ln
P(X,Y)(520)

 当取pj=P(Y|X)，x1=x2=…=xn=P(X)，f为自然对数(凹函数)时

ln∑YP(Y|X)P(X)=∑YP(Y|X)ln P(X)(521)

其中，边缘概率P(X)是与Y取值无关的量(相当于常量c)。

3. 设计迭代算法
任给初始参数θ0，设计求解式(518)的迭代算法，其关键步骤是： 第i次迭代时如何将参数从θi-1更新到θi，使得对数似然函数ln l(θ)的函数值上升，即ln l(θi)≥ln l(θi-1)。
1) 根据参数θi-1确定概率模型
因为上次迭代的参数θi-1为已知，据此可确定联合概率分布P(X,Y； θi-1)、边缘概率分布P(X； θi-1)和后验概率P(Y|X； θi-1)。这三者之间的关系为

P(X,Y； θi-1)=P(Y|X； θi-1)P(X； θi-1)(522)


2) 确定对数似然函数ln l(θ)的下界
对式(517)的对数似然函数ln l(θ)做如下等价变换:

ln l(θ)=ln∑YP(X,Y； θ)=ln∑YP(Y|X； θi-1)P(X,Y； θ)P(Y|X； θi-1)(523)

根据式(519)的Jensen不等式，如果取

pj=P(Y|X； θi-1)，xj=P(X,Y； θ)P(Y|X； θi-1)

则有

ln l(θ)≥∑YP(Y|X； θi-1)lnP(X,Y； θ)P(Y|X； θi-1)(524)

不等式(524)右边的项是对数似然函数ln l(θ)的下界，可记作下界函数B(θ,θi-1)，或称作证据下界目标(Evidence Lower Bound Objective，ELBO)，即

B(θ,θi-1)=∑YP(Y|X； θi-1)lnP(X,Y； θ)P(Y|X； θi-1)(525)

将优化目标由最大化式(518)的对数似然函数改为最大化式(525)的下界函数B(θ,θi-1)，求解最优参数θ*，即

θ*= argmaxθ B(θ,θi-1)(526)

其中,θi-1是上一轮迭代的参数，为已知量。将θ*作为迭代算法的新参数θi，则有

θi=θ*，ln l(θi)= maxθ B(θ,θi-1)(527)


3) 证明ln l(θi)≥ln l(θi-1)
根据式(523)，若取θ=θi-1，则

ln l(θi-1)=ln∑YP(X,Y； θi-1)=ln∑YP(Y|X； θi-1)P(X； θi-1)(528)

根据式(519)的Jensen不等式，如果取

pj=P(Y|X； θi-1)，x1=x2=…=xn=P(X； θi-1)

则Jensen不等式取等号，即

ln l(θi-1)=ln∑YP(Y|X； θi-1)P(X； θi-1)
=∑YP(Y|X； θi-1)lnP(X； θi-1)
=∑YP(Y|X； θi-1)lnP(X,Y； θi-1)P(Y|X； θi-1)=B(θi-1,θi-1)(529)

由式(527)可知:

ln l(θi)= maxθ B(θ,θi-1)≥B(θi-1,θi-1)=ln l(θi-1)(530)


4) Q函数
对于式(526)的下界函数B(θ,θi-1)的最优化问题做进一步整理： 

θ*= argmaxθ B(θ,θi-1)= argmaxθ∑YP(Y|X； θi-1)lnP(X,Y； θ)P(Y|X； θi-1)
= argmaxθ∑YP(Y|X； θi-1)lnP(X,Y； θ)-∑YP(Y|X； θi-1)lnP(Y|X； θi-1)
= argmaxθ∑YP(Y|X； θi-1)lnP(X,Y； θ)(531)

记

Q(θ,θi-1)=∑YP(Y|X； θi-1)lnP(X,Y； θ)(532)

其中，函数Q(θ,θi-1)可看作完全数据(X,Y)的对数似然函数ln
P(X,Y； θ)在概率分布P(Y|X； θi-1)下的期望，简称Q函数。
式(526)的下界函数B(θ,θi-1)的最优化问题与函数Q(θ,θi-1)的最优化问题等价，即

θ*= argmaxθ B(θ,θi-1)= argmaxθ Q(θ,θi-1)(533)

式（533）是对式（526）的简化，将其最优参数θ*作为迭代算法的新参数θi，同样有

ln l(θi)≥ln l(θi-1)


4. EM算法步骤
假设随机变量X、Y服从参数为θ的概率分布P(X,Y； θ)，其中Y为不可观测或未被观测的隐变量。如果已知联合概率P(X,Y； θ)的分布形式，或其等价形式，例如已知边缘概率P(Y； θ)、条件概率P(X|Y； θ)的分布形式，但分布中的参数θ未知。这时可以根据不完全数据X的样本，使用EM算法估计参数θ，从而建立起随机变量X、Y的概率模型P(X,Y； θ)。
给定未做标注的数据集D={x1,x2,…,xm}，其对数似然函数为

lnL(θ)=ln∏mj=1P(xj； θ)=∑mj=1ln∑y∈ΩYP(xj,y； θ)(534)

使用EM算法最大化式(534)的对数似然函数lnL(θ)，求解最优参数θ*，其算法步骤为： 首先选择初始参数θ0，然后迭代执行如下的E步和M步(EM算法的名称正是来源于此)，直至收敛(例如迭代后参数θi与θi-1无明显变化)。
E步： 根据上次迭代的参数θi-1，求完全数据(X,Y)的对数似然函数lnP(X,Y； θ)在概率分布P(Y|X； θi-1)下的期望。这里的E指的就是期望(expectation)，也即Q函数。数据集D上的Q函数为

Q(θ,θi-1)=∑mj=1∑y∈ΩYP(y|xj； θi-1)lnP(xj,y； θ)(535)

其中，θi-1为已知参数，θ为待求解参数。
M步： 最大化期望。这里的M指的就是最大化(maximization)期望，即最大化Q函数，然后将其最优参数作为迭代后的新参数θi，即

θi= argmaxθ Q(θ,θi-1)(536)

算法收敛后，将最后一次迭代的参数θi作为最优参数θ*。联合概率P(X,Y； θ)的分布形式是已知的，因此将θ*代入其中，这就建立起了随机变量X、Y的概率模型P(X,Y； θ*)。所建立的概率模型P(X,Y； θ*)今后可以用于聚类、分类或其他用途。

5.2.2高斯混合模型
高斯混合模型(Gaussian Mixture Model，GMM)是一种广泛使用的概率模型。该模型由多个高斯分布(即正态分布)混合而成，需使用EM算法来估计模型参数(或称训练模型)。

1. GMM模型描述
高斯混合模型假设有n个类别{c1,c2,…,cn}，将分类特征记作随机变量X，类别记作随机变量Y。其中，类别Y是不可观测的隐变量且服从多项分布，其概率分布P(Y； α)为

P(Y=ck)=αk，记作P(Yk； αk)，k=1,2,…,n(537)

其中，α=(α1,α2,…,αn)为未知参数，且∑nk=1αk=1。再假设各类的特征X都服从高斯分布(即正态分布)，但所取参数βk=(μk,σ2k)不同，它们的概率密度函数可记作

p(x|Yk； βk)=12πσke-(x-μk)22σ2k，k=1,2,…,n(538)

其中，β=(β1,β2,…,βn)为未知参数。
由类别概率分布(先验概率)P(Y； α)和各类的特征概率分布P(X|Yk； βk)，可以计算出联合概率P(X,Y； α,β)、分类特征X的边缘概率P(X； α,β)、类别的条件概率(后验概率)P(Y|X； α,β)，即

P(X,Yk； αk,βk)=P(X|Yk； βk)P(Yk； αk)，k=1,2,…,n(539)

P(X； α,β)=∑nk=1P(X|Yk； βk)P(Yk； αk)
=α1P(X|Y1； β1)+α2P(X|Y2； β2)+…+αnP(X|Yn； βn)(540)

P(Yk|X； α,β)=P(X,Yk； αk,βk)P(X； α,β)，k=1,2,…,n(541)

给定未做标注的数据集D=x1,x2,…,xm，其中只包含分类特征X，未包含类别标注Y，因此是一个不完全数据集。现在希望根据数据集D建立特征X的概率模型P(X； α,β)，这是一个含隐变量的
混合
概率模型（即GMM），其概率分布形式已知，但参数α、β未知。

2. GMM模型参数估计
给定未做标注的数据集D={x1,x2,…,xm}，其对数似然函数为

ln L(α,β)=ln∏mj=1p(xj； α,β)

=ln∏mj=1∑nk=1p(xj|Yk； βk)P(Yk； αk)(542)

将式（537）、式（538）代入式（542），整理可得

ln L(α,β)=∑mj=1ln∑nk=1αk·12πσke-(xj-μk)22σ2k(543)

最大化式(543)的对数似然函数ln L(α,β)，需使用EM算法求解最优参数α*、β*，从而建立起GMM的概率模型
P（X； α*，β*）。

3. 应用EM算法
记θ=(α,β)，首先选择初始参数θ0=(α0,β0)，然后迭代执行如下的E步和M步，直至收敛(例如迭代后参数无变化)。
E步： 根据上次迭代的参数θi-1，求数据集D上的Q函数(即期望)。

Q(θ,θi-1)=∑mj=1∑nk=1P(Yk|xj； θi-1)lnP(xj,Yk； θ)(544)

其中

P(Yk|xj； θi-1)=P(xj,Yk； θi-1)P(xj； θi-1)=αi-1k·12πσi-1ke-(xj-μi-1k)22(σi-1k)2∑nl=1αi-1l·12πσi-1le-(xj-μi-1l)22(σi-1l)2(545)

lnP(xj,Yk； θ)=lnαk·12πσke-(xj-μk)22σ2k

=lnαk+ln12π-12lnσ2k-12σ2k(xj-μk)2(546)


为便于后续演算，记γjk≡P(Yk|xj； θi-1)，γjk满足∑nk=1γjk=γj1+γj2+…+γjn=1。将γjk和式（546）代入式（544）可得

Q(θ,θi-1)=∑mj=1∑nk=1γjklnαk+ln12π-12lnσ2k-12σ2k(xj-μk)2(547)

E步结束。

M步： 最大化Q函数(即最大化期望)，然后将其最优参数作为迭代后的新参数θi。
对式(547)的Q函数求参数μk的偏导并令其等于0。

Q(θ,θi-1)μk=∑mj=1γjk1σ2k(xj-μk)=0(548)

求得最优参数μk为

μk=∑mj=1γjkxj∑mj=1γjk，k=1,2,…,n(549)

再对Q函数求参数σ2k的偏导并令其等于0。

Q(θ,θi-1)σ2k=∑mj=1γjk-12σ2k+12σ4k(xj-μk)2=0(550)

求解最优参数σ2k为

∑mj=1γjk-1+1σ2k(xj-μk)2=0

σ2k=∑mj=1γjk(xj-μk)2∑mj=1γjk，k=1,2,…,n(551)

再根据约束条件∑nk=1αk=1，使用拉格朗日乘子法求解最优参数αk。首先构造Q函数的拉格朗日函数LQ(θ,θi-1)，则

LQ(θ,θi-1)=∑mj=1∑nk=1γjklnαk+ln12π-12lnσ2k-12σ2k(xj-μk)2+

λ∑nk=1αk-1(552)

然后对LQ(θ,θi-1)求参数αk的偏导并令其等于0，即

LQ(θ,θi-1)αk=∑mj=1γjkαk+λ=0(553)

λαk=-∑mj=1γjk，k=1,2,…,n(554)

将n个λαk累加起来，有

λα1+λα2+…+λαn=-∑mj=1γj1+∑mj=1γj2+…+∑mj=1γjn

λ(α1+α2+…+αn)=-∑mj=1(γj1+γj2+…+γjn)

因为∑nk=1αk=1，∑nk=1γjk=1，所以λ=-m，将λ代入式(554)可得

αk=1m∑mj=1γjk，k=1,2,…,n(555)

M步结束。综合式(549)、式(551)和式(555)，取新的迭代参数θi为

(αk,μk,σ2k)，k=1,2,…,n

检查迭代条件，如果参数θi与θi-1无明显变化，或i达到最大迭代次数，则停止迭代； 否则返回E步，继续下次迭代。迭代结束后，将最后一次迭代的参数θi作为最优参数θ*。

5.2.3三硬币模型
高斯混合模型是连续型混合概率模型的代表，三硬币模型则是离散型混合概率模型的代表。假设有三枚硬币，分别记作A、B、C，它们出现正面的概率分别为α、p、q。将硬币正面记作1，反面记作0，然后进行如下三硬币实验： 先掷硬币A，根据结果选择硬币B或C，正面选B，反面选C； 再掷所选出的硬币B或C，若为正面则记录1，反面则记录0； 独立重复m次实验，所记录结果为一个0、1序列，记作{x1,x2,…,xm}。
三硬币模型假设硬币的投掷过程不可见，只能看到最终的记录结果{x1,x2,…,xm}。

1. 模型描述
将硬币A的投掷结果记作随机变量Y(可看作聚类问题中不可观测的类别)，所选出硬币(B或C)的投掷结果记作随机变量X(可看作聚类问题中可观测的分类特征)。随机变量Y、X均服从01分布，记Y=1为Y1，Y=0为Y0； X=1为X1，X=0为X0； 再记参数θ=(α,p,q)，则

P(Y1)=α，P(Y0)=1-α，记作P(Yk； θ)，k=0,1
P(X1|Y1)=p，P(X0|Y1)=1-p，记作P(X|Y1； θ)
P(X1|Y0)=q，P(X0|Y0)=1-q，记作P(X|Y0； θ)(556)

其中，θ=(α,p,q)为未知参数。
由概率分布P(Yk； θ)、P(X|Y1； θ)和P(X|Y0； θ)，可以计算出联合概率P(X,Y； θ)、X的边缘概率P(X； θ)、Y的条件概率P(Y|X； θ)，即

P(X,Y； θ)=P(X1,Y1； θ)=P(X1|Y1； θ)P(Y1； θ)=αp
P(X0,Y1； θ)=P(X0|Y1； θ)P(Y1； θ)=α(1-p)
P(X1,Y0； θ)=P(X1|Y0； θ)P(Y0； θ)=(1-α)q
P(X0,Y0； θ)=P(X0|Y0； θ)P(Y0； θ)=(1-α)(1-q)(557)

P(X； θ)=P(X1； θ)=∑1k=0P(X1|Yk； θ)P(Yk； θ)=αp+(1-α)q
P(X0； θ)=∑1k=0P(X0|Yk； θ)P(Yk； θ)=α(1-p)+(1-α)(1-q)(558)

P(Y|X； θ)=P(Y1|X1； θ)=P(X1,Y1； θ)P(X1； θ)=αpαp+(1-α)q
P(Y0|X1； θ)=1-P(Y1|X1； θ)
P(Y1|X0； θ)=P(X0,Y1； θ)P(X0； θ)=α(1-p)α(1-p)+(1-α)(1-q)
P(Y0|X0； θ)=1-P(Y1|X0； θ)(559)

将掷硬币实验所记录的0、1序列看作数据集D={x1,x2,…,xm}，xj∈{0,1}，j=1,2,…,m，其中只包含第二枚硬币(B或C)的投掷结果X。现在希望根据数据集D建立三硬币的概率模型P(X,Y； θ)，这是一个含隐变量的概率模型，其概率分布形式已知，但参数θ=(α,p,q)未知。

2. 模型参数估计
给定掷硬币实验的数据集D={x1,x2,…,xm}，其对数似然函数为

ln L(θ)=ln∏mj=1P(xj； θ)=∑mj=1ln P(xj； θ)(560)

最大化式（560）的对数似然函数ln L(θ)，需使用EM算法求解最优参数θ*=(α*,p*,q*)，从而建立起三硬币的概率模型。

3. 应用EM算法
首先选择初始参数θ0=(α0,p0,q0)，然后迭代执行如下的E步和M步，直至收敛(例如迭代后参数无变化)。
E步： 根据上次迭代的参数θi-1，求数据集D上的Q函数(即期望)。
首先，将参数θi-1=(αi-1,pi-1,qi-1)代入式（559），求得P(Y|X； θi-1)。为便于后续演算，这里对式（559）做进一步改写。记P(Y1|X1； θ)为γ1，P(Y1|X0； θ)为γ0，则

P(Y|X； θ)=P(Y1|X1； θ)=γ1=αpαp+(1-α)q
P(Y0|X1； θ)=1-γ1
P(Y1|X0； θ)=γ0=α(1-p)α(1-p)+(1-α)(1-q)
P(Y0|X0； θ)=1-γ0(561)

其中，γ1、γ0分别是第二枚硬币(B或C)投掷结果X为正面或反面时，硬币A投掷结果为正面的概率。代入参数θi-1=(αi-1,pi-1,qi-1)，求得γi-11和γi-10。



γi-11=αi-1pi-1αi-1pi-1+(1-αi-1)qi-1
γi-10=αi-1(1-pi-1)αi-1(1-pi-1)+(1-αi-1)(1-qi-1)(562)

再假设数据集D={x1,x2,…,xm}中正面样本点有m1个，反面样本点有m0个，
m1+m0=m，则数据集D上的Q函数为

Q(θ,θi-1)=∑mj=1∑1k=0P(Yk|xj； θi-1)lnP(xj,Yk； θ)
=m1Q1(θ,θi-1)+m0Q0(θ,θi-1)(563)

其中

Q1(θ,θi-1)=∑1k=0P(Yk|X1； θi-1)lnP(X1,Yk； θ)
=P(Y1|X1； θi-1)lnP(X1,Y1； θ)+P(Y0|X1； θi-1)lnP(X1,Y0； θ)
=γi-11ln(αp)+(1-γi-11)ln((1-α)q)(564)

Q0(θ,θi-1)=∑1k=0P(Yk|X0； θi-1)lnP(X0,Yk； θ)
=(Y1|X0； θi-1)lnP(X0,Y1； θ)+P(Y0|X0； θi-1)lnP(X0,Y0； θ)
=γi-10ln(α(1-p))+(1-γi-10)ln((1-α)(1-q))(565)

E步结束。

M步： 最大化Q函数(即最大化期望)，然后将其最优参数作为迭代后的新参数θi。
对式(563)的Q函数求参数α的偏导并令其等于0，即

Q(θ,θi-1)α=m1γi-11α-(1-γi-11)(1-α)+m0γi-10α-(1-γi-10)(1-α)=0(566)

求得最优参数α为

α=1m1+m0(m1γi-11+m0γi-10)=1m∑1k=0mkγi-1k(567)

再对Q函数求参数p的偏导并令其等于0，即

Q(θ,θi-1)p=m1γi-11p-m0γi-101-p=0(568)

求得最优参数p为

p=m1γi-11m1γi-11+m0γi-10=m1γi-11∑1k=0mkγi-1k(569)

最后对Q函数求参数q的偏导并令其等于0，即

Q(θ,θi-1)q=m11-γi-11q-m01-γi-101-q=0(570)

求得最优参数q为

q=m1(1-γi-11)(m1+m0)-(m1γi-11+m0γi-10)=m1(1-γi-11)m-∑1k=0mkγi-1k(571)

综合式（567）、式（569）和式（571），取新的迭代参数θi为(α,p,q)，M步结束。检查迭代条件，如果参数θi与θi-1无明显变化，或i达到最大迭代次数，则停止迭代； 否则返回E步，继续下次迭代。迭代结束后，将最后一次迭代的参数θi作为最优参数θ*。
三硬币模型是基于01(二项)分布的概率模型，其求解算法可推广至多项分布。

5.3k均值聚类
与之前基于概率的聚类方法不同，k均值聚类(kmeans clustering)是一种基于距离的聚类方法，其中k表示类别的个数。假设有k个类，每个类有一个中心点。直观上看，样本点离哪个类的中心点距离近就应该划归哪个类。
k均值聚类是一种无监督学习算法。给定未做标注的数据集D={x1,x2,…,xm}，k均值聚类需要将其划分成k个不相交的簇，可记作{C1,C2,…,Ck}。首先为每个簇建立一个能够代表该簇的原型(prototype，即中心点)，然后将各样本点划归距离最近原型所代表的簇。k均值聚类以簇中样本的均值作为该簇的原型，每个簇构成一类，共k个。将k个类的均值记作均值向量μ=(μ1,μ2,…,μk)，它们是聚类模型的未知参数，需基于数据集D并设计聚类算法来进行学习。
k均值聚类算法主要包括数据预处理、距离度量、均值初始化、均值迭代和数据集聚类等五步，其中距离度量和均值迭代是算法的关键。

5.3.1k均值聚类算法
1. 数据预处理

数据集D通常包含多个特征项，不同特征项的取值范围可能存在较大差异，这就会造成特征项之间的不平等。为防止这种不平等现象被代入后续的距离度量，k均值聚类算法需通过预处理对数据集D进行标准化，也即消除量纲，统一不同特征项的取值范围。常用的数据标准化方法有MinMax、zscore等。

2. 距离度量
对于数值型特征，度量样本点之间的距离通常采用欧氏距离，即

dist(xi,xj)=
∑dl=1
(xil-xjl)2≡‖xi-xj‖2(572)

其中，d表示特征的维数； 样本点xi、xj之间的欧氏距离通常也被记作‖xi-xj‖2。
对于非数值型特征，目前还没有很好的距离度量方法。例如，如何度量苹果、香蕉、桔子之间的距离，或如何度量马、牛、羊之间的距离？k均值聚类方法主要适用于数值型特征的聚类问题。可以将非数值型特征转换为数值型，然后按数值型特征进行处理。
在确定了距离度量之后，k均值聚类算法将聚类模型在数据集D上的损失函数定义为各样本点与其所属类均值之间距离的平方和，即

L(μ)=∑mj=1dist2(xj,μl)=∑mj=1‖xj-μl‖22(573)

其中，μl表示样本点xj属于类Cl，即xj∈Cl，μl为该类的均值。可以看出损失函数L(μ)正比于各类方差的总和(即类内方差)，即

L(μ)=∑mj=1‖xj-μl‖22∝∑kl=1σ2l(574)

最小化式（573）的损失函数L(μ)，相当于最小化类内方差。损失函数L(μ)中各样本点xj的所属类别Cl为隐变量，求解最优参数μ*需通过类似EM的迭代算法进行求解。

3. 均值初始化
记样本均值μ=(μ1,μ2,…,μk)，从标准化后的数据集D中随机选取k个样本点，分别作为k个类的初始均值，将它们记作μ0=(μ01,μ02,…,μ0k)。

4. 均值迭代
有了初始均值μ0，下面从i=1开始迭代执行如下的标注样本和更新均值两步，直至收敛(例如迭代后均值无变化)。
1) 标注样本
对数据集D中的样本数据进行标注。根据上次迭代的均值μi-1=(μi-11,μi-12,…,μi-1k)，将各样本点标注为距离最近均值所代表的类，即

ci-1j= argminl=1,2,…,k‖xj-μi-1l‖22，j=1,2,…,m(575)

其中，xj表示第j个样本点； ci-1j表示该样本点本次迭代的标注结果，ci-1j∈{1,2,…,k}。然后再根据标注结果，将标注相同的样本点划归一类，这样数据集D就被划分成k个类，记作Ci-1={Ci-11,Ci-12,…,Ci-1k}。

2) 更新均值
根据上一步的分类结果Ci-1={Ci-11,Ci-12,…,Ci-1k}，重新计算各类的均值，即更新均值。

μil=1|Ci-1l|∑xj∈Ci-1lxj，l=1,2,…,k(576)

其中，|Ci-1l|表示类Ci-1l中样本点的个数。更新均值后检查迭代条件，如果参数μi与μi-1无明显变化，或i达到最大迭代次数，则停止迭代，将参数μi作为聚类模型的最优参数μ*； 否则返回上一步，继续迭代。

5. 数据集聚类
通过均值迭代求解出的最优参数μ*=(μ*1,μ*2,…,μ*k)，它们是最终聚类模型中各类的样本均值(即原型)。按照与式（575）相同的距离最近原则，对数据集D中样本数据做最终标注，将样本数据xj标注为yj，即

yj= argminl=1,2,…,k‖xj-μ*l‖22，j=1,2,…,m

这样就完成了对数据集的聚类过程。今后任给新样本x，同样也可以按上述方法进行标注(即分类)。

5.3.2关于k均值聚类的讨论

1. k均值聚类算法的收敛性
k均值聚类算法即最小化式（573）的损失函数L(μ)，其中各样本点xj的所属类别Cl相当于隐变量Y，其求解过程实际上是对常规含隐变量概率模型和EM算法的一种简化。
k均值聚类算法中标注样本的操作(见式（575）)相当于EM算法中的E步： 根据上次迭代的参数μi-1，求数据集D上的Q函数（即期望）。

Q(μ,μi-1)=∑mj=1∑YP(Y|xj； μi-1)lnP(xj,Y； μ)(577)

其中，μi-1为已知参数，μ为待求解参数，且

P(Y|x； μi-1)=P(Yc|x； μi-1)=1，c= argminl=1,2,…,k‖xj-μi-1l‖22

P(Yc|x； μi-1)=0，c≠argminl=1,2,…,k‖xj-μi-1l‖22(578)

lnP(x,Y； μ)=lnP(x,Yc； μ)=-‖x-μc‖22，c= argminl=1,2,…,k‖xj-μi-1l‖22

lnP(x,Yc； μ)=0，c≠argminl=1,2,…,k‖xj-μi-1l‖22

(579)

Q(μ,μi-1)=∑mj=1(-‖xj-μc‖22)=-L(μ)，μc为xj所属类的均值(580)

k均值聚类算法中更新均值的操作(见式（576）)相当于EM算法中的M步： 最大化Q函数，实际上就是最小化损失函数L(μ)，即

μi= argmaxμ Q(μ,μi-1)= argminμ L(μ)(581)

式（581）的最优解就是式（576）的最优解。
由EM算法可知，第i次迭代时将参数从μi-1更新到μi，可使损失函数L(μ)的函数值下降，即L(μi)≤L(μi-1)。如果损失函数L(μ)是凸函数(例如式（573）)，则k均值聚类算法可以收敛至最小值。

2. 超参数k值的选择
k均值聚类中的k是一个需要人工预设的超参数，表示类别的个数。可以通过可视化观察数据分布情况，给出一个相对合理的k值； 也可以为k选取一组候选值，然后使用数据集逐个训练并计算各候选k值最终的损失(见式（573）)，从中选取损失最小的候选值作为最优k值。
一个更加灵活的方法是仅将预设的k值看作预估值，聚类时再根据数据的实际分布进行调整。聚类过程中，当簇内样本点数量过少或两个簇距离过近时，对簇进行合并操作； 当簇内样本点数量过多且方差较大时，对簇进行分裂操作。这种根据数据实际分布动态调整k值的方法被称作ISODATA(Iterative SelfOrganizing Data Analysis Techniques Algorithm)聚类。

3. 初始均值的选择
k均值聚类的k个初始均值是从数据集D中随机选取的。从直观上看，这k个初始均值应尽量散开，这样可以加快收敛速度。随机选取第一个初始均值μ01，然后在选择下一个均值μ02时尽量选择距μ01较远的样本点，或者说距μ01越远的样本点被选作μ02的概率应当越大； 重复这个过程，在选择下一个均值μ0i时，距前i-1个{μ01,μ02,…,μ0i-1}越远的样本点被选作μ0i的概率应当越大，直至选出全部k个均值。按照上述方法选取初始均值的k均值聚类方法被称为
kmeans++聚类。

4. 聚类模型的评价指标
评价聚类模型的好坏，轮廓系数(silhouette coefficient)是一个常用的评价指标。对于聚类模型在数据集D={x1,x2,…,xm}的聚类结果，首先计算各样本点xi距本簇中其他样本的平均距离ai(越小越好)，以及距其他最近一个簇中样本间的平均距离bi(越大越好)，然后定义样本点xi的轮廓系数si为

si=bi-aimax(ai,bi)，-1≤si≤1(582)

最终，聚类模型在整个数据集D上的轮廓系数s被定义为所有样本点轮廓系数的平均，即

s=1m∑mi=1si，-1≤s≤1(583)

轮廓系数s越接近于1，则说明聚类的效果越好。

5.3.3使用scikitlearn库中的k均值聚类模型
练习k均值聚类编程，可以使用scikitlearn库提供的葡萄酒数据集wine recognition dataset和k均值聚类模型KMeans。

1. 葡萄酒数据集wine recognition dataset
wine recognition dataset是一个包含人工标注的葡萄酒数据集，其中包含178个样例(产自3个不同的种植园)，每个样例包含13项特征数据(表示13种化学成分含量，均为数值型特征)和一项人工标注(标注来自哪个种植园，分别用0、1、2表示)。
图51给出一个下载葡萄酒数据集的示例代码，它将数据集及其说明文档分别保存到data目录下的wine.csv文件和wine.txt中。



图51下载葡萄酒数据集的示例代码



仅仅依据产自哪个种植园来标注葡萄酒，这种分类方式是科学的吗？可以根据葡萄酒的化学成分进行聚类，对比一下聚类的标注结果与人工标注是否一致。如果两者比较接近，那就说明根据产地对葡萄酒进行分类是科学的。

2. 使用scikitlearn库中的k均值聚类模型
scikitlearn库将k均值聚类模型封装成一个类(类名为KMeans)，并将其存放在sklearn.cluster模块中。KMeans类实现了k均值聚类(含kmeans++)模型的聚类算法fit()，另外还提供一个预测算法predict()，可用于对新样本的分类。
首先从下载到本地data目录下的wine.csv文件中加载葡萄酒数据集，得到一个DataFrame类的二维表格w3，显示其形状(shape，即表格的行数和列数)。取出数据集中的特征(0~12列)，对其进行标准化(例如使用scikitlearn库提供的MinMaxScaler类实现MinMax标准化)，生成一个仅包含特征的数据集X； 再取出人工标注(第13列，列名为target)，生成标注集Y。图52给出了加载本地葡萄酒数据集wine.csv并进行标准化的示例代码。



图52加载本地葡萄酒数据集并进行标准化的示例代码


图53给出了使用KMeans类对葡萄酒数据集进行聚类的示例代码。其中聚类所使用的数据集是图52生成的特征数据集X，k值取3； 调用fit()方法得到聚类结果，其中包括三个类的样本均值cluster_centers_，以及数据集X中各样例的聚类标注labels_(0、1或2)。可以同时显示图52标注集Y中的人工标注，将其与KMeans类标注进行比对，如图54所示。





图53使用KMeans类对葡萄酒数据集进行聚类的示例代码






图54比对KMeans类标注labels_与人工标注labels_Y


图54中KMeans类标注labels_的2、0、1(仅仅是一种类别编号)分别对应人工标注labels_Y的0、1、2。可以看出，两者的分类结果(即哪些样例是属于同一类的)很接近，这说明不同产地葡萄酒的化学成分确实不一样，以产地来衡量葡萄酒质量是有一定道理的。


3. 使用scikitlearn库中的轮廓系数函数
scikitlearn库为评价聚类模型提供了一个计算轮廓系数的函数silhouette_score()，它被存放在sklearn.metrics模块中。图55给出了一段示例代码，用于计算并显示k均值聚类模型在葡萄酒数据集上的轮廓系数。



图55k均值聚类模型在葡萄酒数据集上的轮廓系数



5.4密度聚类DBSCAN
与k均值聚类基于样本距离的方法不同，DBSCAN(DensityBased Spatial Clustering of Applications with Noise)是一种基于样本密度(density)的聚类方法。从直观上看，同类的样本应集中分布在同一区域内，区域内部的样本密度高，边缘的样本密度低。聚类时可以从某个内部点出发逐步向周边扩展，直至遇到边缘点时停止； 重复这个扩展过程，最终扩展至类的整个分布区域。
如果有多个类，不同类的样本应分布在不同区域，区域之间应该有间隙，或者只有少量被称作噪声(noise)或离群点(outlier)的样本。图56给出一个样本分布示意图，其中有三个类，分别用●、▲、╋表示，另外还用★表示相对孤立的噪声样本。



图56三个类的样本分布示意图


DBSCAN是一种典型的基于密度的聚类算法，它能在具有噪声的数据集中发现任意形状的簇，每个簇构成一类。下面结合图56来具体讲解DBSCAN聚类的主要术语及算法实现。

5.4.1DBSCAN聚类术语
给定未做标注的数据集D={x1,x2,…,xm}，下面定义一下DBSCAN聚类用到的主要术语。

1. ε邻域
对于样本点xj∈D，数据集D中与xj距离不超过ε的样本点集合被称作xj的ε邻域(neighborhood)，记作Nε(xj)，

Nε(xj)={xi|xi∈D,distance(xi,xj)≤ε}(584)

这个ε邻域中样本点的个数|Nε(xj)|被称作xj的样本密度。例如，图56中x的邻域内有5个样本点(不含x)，则x的样本密度为5。式（584）中的ε是一个需要人工预设的超参数； 距离distance(xi,xj)可以使用不同的度量形式，常用的是欧氏距离。

2. 核心对象
若xj的样本密度不小于MinPts个样本点，即|Nε(xj)|≥MinPts，则xj被称作核心对象(core object)。其中，MinPts是核心对象的阈值，它是一个需要人工预设的超参数。如果xj是一个核心对象，则说明xj一定属于某个类且位于该类样本分布区域的内部。例如，假设MinPts=3(下同)，则图56中的x、x1、x3等就是核心对象，但x4、x6、x7等就不是核心对象。

3. 密度直达
若xi位于核心对象xj的邻域内，则称xi可由xj密度直达(directly densityreachable)。例如，图56中的x1、x3可由核心对象x密度直达，另外x2也可由核心对象x1密度直达。若xi、xj均为核心对象，则属于双向密度直达； 若其中一个为核心对象，另一个为非核心对象，则属于单向密度直达； 若两者均为非核心对象，则不可能密度直达。例如，x3与x之间属于双向密度直达(两者均为核心对象)； 而x3与x4之间属于单向密度直达(x4为非核心对象)； x4与x6之间则不可能密度直达(两者均为非核心对象)。

4. 密度相连
若存在样本点序列x,x1,x2,…,xn,y，其中x1可由x密度直达，y可由xn密度直达，xi均为核心对象且与xi+1之间双向密度直达，则称x与y之间经由x1,x2,…,xn密度相连(densityconnected)。例如，x1与x3之间经由x密度相连，x2与x4之间经由x1,x,x3密度相连。密度直达也属于密度相连，是密度相连的特例。

5. 簇
给定超参数(ε,MinPts)，数据集D中每个密度相连的最大子集即构成一个簇。换句话说，假设数据集D中的某个核心对象xj属于某个簇C，如果xi∈D且xi与xj密度相连，则xi∈C。可以选择一个核心对象，然后将所有与之密度相连的样本点聚集在一起形成簇，并将它们标注为一类，这就是DBSCAN聚类。

5.4.2DBSCAN聚类算法
给定需要聚类的数据集D={x1,x2,…,xm}，并设定好超参数(ε,MinPts)，DBSCAN聚类算法的步骤如下。

1. 找出核心对象
遍历数据集D，找出其中所有的核心对象，并记作核心对象集合Ω，则

Ω={xj|xj∈D,|Nε(xj)|≥MinPts}(585)


2. 选取核心对象进行扩展
（1) 初始化当前类别标注k=0。
（2) 遍历集合Ω，从中随机选取一个核心对象作为簇的种子(记作s)，然后进行扩展： 将s标记为visited(已访问或已处理)并分别创建一个新簇Ck=s和一个扩展队列Q=Nε(s)，转下一步。
（3) 遍历扩展队列Q，依次取出其中的样本点(记作q)并进行扩展： 将q标记为visited并将其追加到簇Ck中，然后检查q是否为核心对象，若是核心对象则将Nε(q)中所有未被访问过(无visited标记)的样本点追加到队列Q的队尾。
（4） 重复第（3）步，直到扩展队列Q为空。当扩展队列Q为空时，簇Ck完成聚类，将其中所有核心对象从集合Ω中删除(因为已被处理)，然后将k加1，转第（5）步。
（5) 检查集合Ω，如果所有核心对象都已被处理，则集合Ω为空，聚类结束； 否则返回第（2）步，继续选取核心对象并扩展下一个簇。

3. 标注数据集
聚类结束后，根据所扩展出的k个簇C0,C1,…,Ck-1对数据集D中的样本进行标注。将数据集D中属于簇Ci的样本点全部标注为i，剩余样本点(不属于任何簇)标注为噪声(例如标注为-1)。整个聚类算法结束。
与k均值聚类相比，DBSCAN聚类不需要指定k值，但需指定两个阈值(ε,MinPts)，这两个超参数对聚类结果的影响比较大。另外，DBSCAN聚类时没有建立模型，只能用于当前数据集D的聚类，不能再对其他新样本进行分类。DBSCAN聚类最大的优点是可以对非凸样本集(见图57)进行聚类，而k均值聚类只适用于凸样本集(例如服从高斯分布的样本集)。DBSCAN聚类对凸样本集和非凸样本集都适用。



图57非凸样本集示例



5.4.3使用scikitlearn库中的DBSCAN聚类算法
scikitlearn库将DBSCAN聚类算法封装成一个类(类名也为DBSCAN)，并将其存放在sklearn.cluster模块中。DBSCAN类的最主要方法就是聚类算法fit()。由于DBSCAN在聚类时没有建立模型，因此DBSCAN类也没有对新样本进行分类的预测算法predict()。
图58给出了使用DBSCAN类对葡萄酒数据集进行聚类的示例代码。其中聚类所使用的数据集是图52标准化后的特征数据集X； 调用fit()方法得到聚类结果，其中包括数据集X中的核心对象core_sample_indices_，以及各样例的聚类标注labels_(0、1或-1)。可以同时显示图52标

图58使用DBSCAN类对葡萄酒数据集进行聚类的示例代码

注集Y中的人工标注，将其与DBSCAN聚类标注进行比对，也可以显示DBSCAN聚类在葡萄酒数据集上的轮廓系数，如图59所示。



图59比对DBSCAN聚类标注labels_与人工标注labels_Y



5.5向量量化
给定未做标注的数据集D={x1,x2,…,xm}，k均值聚类希望将其划分成k个不相交的簇，并以簇中样本的均值作为簇的原型。向量量化(Vector Quantization，VQ)与k均值聚类有点类似，所不同的是向量量化在聚类后将k个簇的均值向量(μ1,μ2,…,μk)作为今后对向量进行编码的码本(codebook)，记作C=(μ1,μ2,…,μk)，其中每个均值μj被称作一个码字向量(code vector，简称码字)。编码时，任给向量x，向量量化将码本中距离最近的均值μj作为向量x的编码； 或将距离最近均值μj的编号j作为向量x的编码，这样可以减少向量存储或传输时的数据量。向量量化被广泛应用于向量的离散化编码(称作量化)，或向量数据的压缩。

5.5.1向量量化问题
向量量化分两个环节： 一是训练码本，即给定样本数据集，通过聚类算法建立码本； 二是编码，即任给向量，将码本中距离最近的码字(或其编号)作为向量的编码。

1. 向量量化的术语与符号
假设需要量化的向量为d维向量，下面给出量化过程中常用的术语和符号。
 训练集： 用于训练码本的样本数据集，记作T={x1,x2,…,xm}，其中每个样本数据xi均为一d维向量。
 训练码本： 将训练集T聚类成k个不相交的簇，将簇的集合记作S={S1,S2,…,Sk}，其中Sj表示第j个簇； 簇Sj以其样本均值作为簇的码字，将码字记作cj(也是d维向量)，例如若Sj包含mj个样本数据，则

cj=1mj∑x∈Sjx，j=1,2,…,k

所有码字的集合被称作码本，记作C={c1,c2,…,ck}。对训练集进行聚类并生成码本的算法被称作向量量化算法，简称VQ算法。
 编码： 任给向量x，将码本中距离(例如欧氏距离)最近的码字cj(或其编号j)作为向量x的编码，记作q(x)。上述编码准则被称作最近邻准则，其数学形式可表示为

q(x)= argmincj∈C‖x-cj‖2

 失真度： 聚类通常采用欧氏距离来度量向量之间的距离。在向量量化中，向量与其编码之间不完全相等，例如向量x与其编码q(x)之间会存在误差，q(x)≈x。换句话说，对向量编码会造成失真(distortion)。可以用欧氏距离(或其平方)来度量编码的失真程度，即失真度，记作d(x,q(x))≡‖x-q(x)‖2。
 平均失真度： 如果使用码本C对训练集T={x1,x2,…,xm}进行编码，可以使用平均失真度来度量码本C在训练集T上的失真程度，记作DC。

DC=1m∑mi=1d(xi,q(xi))=1m∑mi=1‖xi-q(xi)‖2

 量化准则： 通过训练集T训练码本，所训练出的码本C*应遵循最小失真准则。其含义是，使用码本C*对训练集T进行编码，其平均失真度应当最小。即

C*= argminC DC(586)


2. 进一步认识向量量化问题
向量量化虽然与k均值聚类有些相似，但向量量化问题远比聚类问题深刻。向量量化的原始问题是： 给定向量空间Ω(这里假设为连续型)，其向量取值服从概率分布p(x)，向量量化希望找到一个最优码本C*Ω，使得向量x∈Ω在编码后失真度的期望最小，即

C*Ω= argminC E［d(x,q(x))］= argminC∫Ωd(x,q(x))p(x)dx(587)

为训练码本，向量量化需要先依概率分布p(x)进行抽样，得到训练集T={x1,x2,…,xm}，然后通过VQ算法求得训练集T上的最优码本C*(见式（586）)，并将其作为整个向量空间Ω上最优码本C*Ω的估计。
相比较而言，k均值聚类算法并不关心训练集T从哪里来，聚类结果会被推广哪里去，它只关心训练集T本身的聚类效果。另外，k均值聚类算法关注的是分类效果，其算法目标是让类内方差最小，类间方差最大； 而VQ算法关注的是编码效果，其算法目标只追求类内方差最小，并不需要类间方差最大。由于上述区别，VQ算法与k均值聚类算法在初始均值的选择上存在较大不同。

5.5.2LBGVQ算法
LBGVQ算法是向量量化中比较常用的一种码本训练算法(或称码本学习算法)，它是1980年由Linde、Buzo和Gray三人联合提出的。LBGVQ算法与k均值聚类算法有点类似，它也是一种迭代算法。所不同的是，k均值聚类算法是一次性选择k个初始均值，然后在此基础上进行迭代； 而LBGVQ算法则是从一个初始均值(即初始码字)开始，然后在此基础上先做码字分裂(split)，再做聚类，重复这个“分裂聚类”过程，直到码字达到指定数量(通常为2N个)。

给定训练集T={x1,x2,…,xm}，假设要训练一个包含2N个码字的码本，则LBGVQ算法就是一个N级“分裂-聚类”的迭代过程，其具体算法步骤如下。
（1） 初始化： 设置码本大小N(包含2N个码字)，例如N=3(2N=8)； 设置失真阈值ε>0，例如ε=0.001； 将初始码字c01设为训练集T的样本均值，并计算出初始平均失真度D0C。

c01=1m∑mi=1xi，D0C=1m∑mi=1‖xi-c01‖2

初始码本只包含一个码字，即C0={c01}，设置码字计数器k=1； 设置迭代计数器t=1，然后进入下面的“分裂聚类”迭代过程。
（2） 分裂(第t轮)： 按如下方式将Ct-1={ct-11,ct-12,…,ct-1k}中的每个码字一分为二。

ctj=(1+ε)ct-1j，ctk+j=(1-ε)ct-1j，j=1,2,…,k

并以新分裂出的2k个码字组成新码本，将k乘以2，即k=2k，将新码本记作Ct。

Ct={ct1,ct2,…,ctk}

（3） 聚类(第t轮)： 使用新码本Ct对训练集T={x1,x2,…,xm}进行编码和聚类，将码本Ct中距离最近(也即失真度最小)的码字作为样本xi的编码，记作q(xi)，即

q(xi)= argminctj∈Ct‖xi-ctj‖2，i=1,2,…,m

并将编码相同的样本聚集成一簇，例如将编码为ctj的样本都聚集到簇Sj中，即

Sj={xi|xi∈T,q(xi)=ctj}，j=1,2,…,k

然后计算码本Ct在训练集T上的平均失真度DtC。

DtC=1m∑mi=1d(xi,q(xi))=1m∑mi=1‖xi-q(xi)‖2

检查聚类迭代条件： 如果第t轮迭代的平均失真度DtC较上一轮改进不大，即

Dt-1C-DtCDtC≤ε

则停止聚类迭代，转第（4）步； 否则重新计算各簇的样本均值并将其作为新的码本Ct+1，例如若Sj包含mj个样本数据，则

ct+1j=1mj∑x∈Sjx，j=1,2,…,k

将迭代计数加1，即t=t+1，返回第（3）步开头继续下一轮聚类迭代。
（4） 算法结束条件： 如果码本已达到指定大小，即码字数量k=2N，则算法结束，将Ct作为最终的码本C*； 否则将迭代计数加1，即t=t+1，返回第（2）步继续下一轮“分裂聚类”迭代，即先分裂，再聚类。

5.6本章习题

一、 单选题
1. 下列关于分类与聚类的描述中，错误的是()。
A.  分类是根据有标注数据集来训练模型，属于有监督学习


B. 聚类是根据无标注数据集来训练模型，属于无监督学习
C.  聚类根据数据自身的分布特性或结构自动聚集成簇，形成类别概念
D.  给定无标注数据集D，聚类算法需将数据集D划分成若干个可重叠的簇
2. 下列概率分布的等价表示中，错误的是()。
A.  P(X=x)≡P(x)B. P(Y|X)≡P(Y|x)
C. P(Y=y|X)≡P(y|X)D. ∑y∈ΩYP(X,y； θ)≡∑YP(X,Y； θ)
3. 在聚类问题中，给定未做标注的数据集D={x1,x2,…,xm}，下列说法中错误的是()。
A.  数据集D只包含样本的分类特征X，未包含样本对应的类别标注Y
B. 分类特征X是可观测的变量，样本类别Y是不可观测
或未被观测
的隐变量
C.  数据集D未包含样本类别Y的原因是实际应用不需要它
D. 聚类问题的关键是如何根据数据集D来估计含隐变量概率模型的参数
4. 关于聚类问题和混合概率模型参数估计问题，下列说法中错误的是()。
A.  它们的模型在本质上都属于含隐变量的概率模型
B. 它们的样本数据集D={x1,x2,…,xm}都不包含类别标注
C.  估计它们的模型参数通常都使用EM算法进行求解
D.  它们最终要求解的都是分类特征X的概率分布P(X)
5. 下列关于EM算法的描述中，错误的是()。
A.  EM算法主要用于求解含隐变量的最优化问题
B. EM算法是一种迭代算法
C.  EM算法的关键步骤是第i次迭代时如何将参数从θi-1更新到θi
D.  EM算法每次迭代应能让对数似然函数逐步下降，即ln l(θi)≤ln l(θi-1)
6. ()是EM算法中的Q函数（即期望）。
A.  ∑YP(Y|X； θi-1)lnP(X,Y； θ)
B. P(Y|X； θi-1)lnP(X,Y； θ)
C.  ∑YP(Y|X； θi-1)lnP(X,Y； θ)P(Y|X； θi-1)  
D. P(Y|X； θi-1)lnP(X,Y； θ)P(Y|X； θi-1)
7. 下列关于EM算法步骤的描述中，错误的是()。
A.  EM算法首先选择初始参数θ0
B. EM算法的E步是根据上次迭代的参数θi-1求Q函数（即期望）
C.  EM算法的M步是最小化Q函数，将最优参数作为迭代后的新参数θi
D.  EM算法需重复执行E步和M步，直至收敛
8. 下列关于混合概率模型的描述中，错误的是()。
A.  高斯混合模型由多个正态分布混合而成的
B. 高斯混合模型是连续型混合概率模型的代表
C.  三硬币模型是离散型混合概率模型的代表
D.  三硬币模型的硬币投掷过程是可观测的
9. 下列关于k均值聚类的描述中，错误的是()。
A.  k均值聚类是一种基于概率的聚类方法
B. k均值聚类是一种基于距离的聚类方法
C.  k均值聚类中的k指的是类别个数
D.  均值聚类是一种无监督学习算法
10. 下列关于k均值聚类的描述中，错误的是()。

A.  k均值聚类中的k是需要人工预设的超参数
B. k均值聚类可通过可视化分析给出一个相对合理的k值
C.  根据数据实际分布动态调整k值的方法被称为ISODATA
D.  让k个初始均值尽量靠近的k均值聚类方法被称为kmeans++
11. 下列关于DBSCAN聚类的描述中，错误的是()。
A.  DBSCAN聚类是一种基于概率的聚类方法
B. DBSCAN聚类是一种基于样本密度的聚类方法
C.  DBSCAN聚类认为同类的样本应集中分布在某个区域内
D.  DBSCAN聚类认为不同类的样本应分布在不同区域，区域之间应该有间隙
12. 在DBSCAN聚类中，若xi位于核心对象xj的邻域内，则称xi可由xj()。
A.  密度直达  B. 密度相连  
C. 密度连通  D. 连通
13. 与k均值聚类相比，DBSCAN聚类()。
A.  不需要指定k值
B. 需要指定两个阈值(ε,MinPts)
C.  除了能对当前数据集D聚类之外，还能对其他新样本进行分类
D.  可以对非凸样本集进行聚类
14. 下列关于向量量化的描述中，错误的是()。
A.  向量量化在聚类后将k个簇的均值向量(μ1,μ2,…,μk)作为码本
B. 均值向量(μ1,μ2,…,μk)中的每个均值μj被称作一个码字
C.  任给向量x，向量量化将码本中距离最近的均值(或其编号)作为其编码
D.  向量量化被广泛应用于向量数据的无失真压缩
15. 下列关于k均值聚类与向量量化的描述中，错误的是()。
A.  k均值聚类只关心训练集本身的聚类效果
B. k均值聚类算法的目标是让类内方差最小，类间方差最大
C.  向量量化算法关注的是编码效果
D.  向量量化算法的目标是让类内方差最大，类间方差最小

二、 讨论题
1. 尝试用概率模型对聚类问题进行形式化描述。
2. 尝试对混合概率模型的参数估计问题进行形式化描述。
3. 尝试推导EM算法中对数似然函数ln l(θ)=ln P(x； θ)的下界函数B(θ,θi-1)。
4. 尝试用EM算法求解三硬币模型的参数估计问题。
5. 简述k均值聚类的算法思想。
6. 简述DBSCAN聚类的算法思想。

三、 编程实践题
使用scikitlearn库提供的鸢尾花数据集(iris plants dataset)设计一个k均值聚类模型。具体的实验步骤如下。
(1) 使用函数sklearn.datasets.load_iris()加载鸢尾花数据集。
(2) 查看数据集说明(https://scikitlearn.org/stable/datasets/toy_dataset.html#irisdataset)，并对数据集进行必要的预处理。
(3) 使用sklearn.cluster.KMeans类建立k均值聚类模型，并对鸢尾花数据集中的特征数据进行聚类处理。
(4) 对比聚类结果与人工标注，观察二者是否一致。