假设我们有一个数据集 { x 1 , . . . , x N } \left\{x_{1},...,x_{N} \right\} {x1,...,xN},它由D维欧几里得空间中的随机变量 x x x的 N N N次观测组成。引入一组 D D D维向量 u k , k = 1 , . . . , K u_{k},k=1,...,K uk,k=1,...,K,对于每个数据点 x n x_{n} xn,我们引入一组对应的二值指示向量 r n k ∈ { 0 , 1 } r_{nk}\in \left\{0,1 \right\} rnk∈{0,1},只有当 x n x_{n} xn被分配到 u k u_{k} uk中, r n k = 1 r_{nk}=1 rnk=1。那么很容易地我们就得到了这样一个目标函数又是也会被称为目标度量 J = ∑ n = 1 N ∑ k = 1 K r n k ∣ ∣ x n − u k ∣ ∣ 2 J=\sum_{n=1}^{N}\sum_{k=1}^{K}r_{nk}||x_{n}-u_{k}||^{2} J=n=1∑Nk=1∑Krnk∣∣xn−uk∣∣2它表示每个数据点与它被分配的向量 u k u_{k} uk之间的距离的平方和。我们可以用迭代的方法来最小化这个目标函数。
- 我们初始化 u k u_{k} uk
- 保持 u k u_{k} uk不变,我们改变 r n k r_{nk} rnk来最小化 J J J
- 保持 r n k r_{nk} rnk不变,我们改变 u k u_{k} uk来最小化 J J J
- 不断重复2,3步骤直到收敛
整体来看,更新
r
n
k
r_{nk}
rnk和更新
u
k
u_{k}
uk的两个阶段分别对应于EM算法中的E(期望)步骤和M(最大化)步骤。
先来看第二步,很好理解,通俗易懂的来说就是把
x
n
x_{n}
xn分配到离它最近的
u
k
u_{k}
uk上。从公式上看也很好理解,
J
J
J是
r
n
k
r_{nk}
rnk的一个线性函数,求其求导得到
∣
∣
x
n
−
u
k
∣
∣
2
||x_{n}-u_{k}||^{2}
∣∣xn−uk∣∣2。也就是说我们要对每个
x
n
x_{n}
xn分配到离它最近的
u
k
u_{k}
uk上,这样才能使
∣
∣
x
n
−
u
k
∣
∣
2
||x_{n}-u_{k}||^{2}
∣∣xn−uk∣∣2最小。用公式可以这样表示
r
n
k
=
{
1
i
f
k
=
a
r
g
m
i
n
j
∣
∣
x
n
−
u
j
∣
∣
2
0
o
t
h
e
r
w
i
s
e
r_{nk}=\left\{\begin{matrix} 1 \qquad if \quad k=arg \ min_{j}||x_{n}-u_{j}||^{2}\\0 \qquad otherwise \qquad \qquad\quad\quad\quad\quad \end{matrix}\right.
rnk={1ifk=arg minj∣∣xn−uj∣∣20otherwise再来看第三步,对
u
k
u_{k}
uk求导=0即
∑
n
=
1
N
r
n
k
(
x
n
−
u
k
)
=
0
\sum_{n=1}^{N}r_{nk}(x_{n}-u_{k})=0
n=1∑Nrnk(xn−uk)=0
u
k
=
∑
n
r
n
k
x
n
∑
n
r
n
k
u_{k}=\frac{\sum_{n}r_{nk}x_{n}}{\sum_{n}r_{nk}}
uk=∑nrnk∑nrnkxn这个结果也个简单的含义即令
u
k
u_{k}
uk等于类别
k
k
k中所有数据点的均值。
python代码如下
import numpy as np
import matplotlib.pyplot as plt
x=1+np.random.randn(20,1)
y=2+np.random.randn(20,1)
z1=np.concatenate((x,y),axis=1)
x=-3-np.random.randn(20,1)
y=1+np.random.randn(20,1)
z2=np.concatenate((x,y),axis=1)
x=-4-np.random.randn(20,1)
y=-6-np.random.randn(20,1)
z3=np.concatenate((x,y),axis=1)
x=1+np.random.randn(20,1)
y=-4-np.random.randn(20,1)
z4=np.concatenate((x,y),axis=1)
z=np.concatenate((z1,z2,z3,z4),axis=0)
z5=z.T
plt.scatter(z5[0],z5[1])
plt.show()
def distEclud(vecA,vecB):
return np.sqrt(np.sum(np.power(vecA-vecB,2)))
def Randcenter(dataset,k):
data=dataset.T
n=dataset.shape[1]
centoids=np.mat(np.zeros((k,n)))
xmax,xmin=max(data[0]),min(data[0])
ymax,ymin=max(data[1]),min(data[1])
dx,dy=(xmax-xmin)/(k+1),(ymax-ymin)/(k+1)
for i in range(k):
centoids[i,0]=xmin+dx*(i+1)
centoids[i,1]=ymin+dy*(i+1)
return centoids
def KMeans(dataset,k,distMeans=distEclud,creatCen=Randcenter):
m=dataset.shape[0]
n=dataset.shape[1]
clusterAssment = np.mat(np.zeros((m,1)))
centoids=creatCen(dataset,k)
flag=True
count=0
while flag:
count+=1
print('第%d次迭代'%count)
flag=False
clusternum=np.mat(np.zeros((k,1)))
newcentoids=np.mat(np.zeros((k,n)))
for i in range(m):
minDist=np.inf
minIndex=-1
for j in range(k):
temp_dist=distMeans(centoids[j],dataset[i])
if temp_dist<minDist:
minDist=temp_dist
minIndex=j
clusterAssment[i]=minIndex
newcentoids[minIndex]+=dataset[i]
clusternum[minIndex]+=1
print(centoids)
for cnt in range(k):
newcentoids[cnt]/=clusternum[cnt]
if distMeans(centoids[cnt],newcentoids[cnt])>1e-10:
flag=True
centoids=newcentoids
return newcentoids,clusterAssment
centers,asse=KMeans(z,4)