作者:mobiledu2502931637 | 来源:互联网 | 2023-10-11 03:07
假设我们有一个数据集{x1,...,xN}\left\{x_{1},...,x_{N} \right\}{x1,...,xN},它由D维欧几里得空间中的随机变量xxx的NNN次观测组成。引入一组DDD维向量uk,k=1,...,Ku_{k},k=1,...,Kuk,k=1,...,K,对于每个数据点xnx_{n}xn,我们引入一组对应的二值指示向量rnk∈{0,1}r_{nk}\in \left\{0,1 \right\}rnk∈{0,1},只有当xnx_{n}xn被分配到uku_{k}uk中,rnk=1r_{nk}=1rnk=1。那么很容易地我们就得到了这样一个目标函数又是也会被称为目标度量J=∑n=1N∑k=1Krnk∣∣xn−uk∣∣2J=\sum_{n=1}^{N}\sum_{k=1}^{K}r_{nk}||x_{n}-u_{k}||^{2}J=n=1∑Nk=1∑Krnk∣∣xn−uk∣∣2它表示每个数据点与它被分配的向量uku_{k}uk之间的距离的平方和。我们可以用迭代的方法来最小化这个目标函数。
- 我们初始化uku_{k}uk
- 保持uku_{k}uk不变,我们改变rnkr_{nk}rnk来最小化JJJ
- 保持rnkr_{nk}rnk不变,我们改变uku_{k}uk来最小化JJJ
- 不断重复2,3步骤直到收敛
整体来看,更新rnkr_{nk}rnk和更新uku_{k}uk的两个阶段分别对应于EM算法中的E(期望)步骤和M(最大化)步骤。
先来看第二步,很好理解,通俗易懂的来说就是把xnx_{n}xn分配到离它最近的uku_{k}uk上。从公式上看也很好理解,JJJ是rnkr_{nk}rnk的一个线性函数,求其求导得到∣∣xn−uk∣∣2||x_{n}-u_{k}||^{2}∣∣xn−uk∣∣2。也就是说我们要对每个xnx_{n}xn分配到离它最近的uku_{k}uk上,这样才能使∣∣xn−uk∣∣2||x_{n}-u_{k}||^{2}∣∣xn−uk∣∣2最小。用公式可以这样表示rnk={1ifk=argminj∣∣xn−uj∣∣20otherwiser_{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再来看第三步,对uku_{k}uk求导=0即∑n=1Nrnk(xn−uk)=0\sum_{n=1}^{N}r_{nk}(x_{n}-u_{k})=0n=1∑Nrnk(xn−uk)=0uk=∑nrnkxn∑nrnku_{k}=\frac{\sum_{n}r_{nk}x_{n}}{\sum_{n}r_{nk}}uk=∑nrnk∑nrnkxn这个结果也个简单的含义即令uku_{k}uk等于类别kkk中所有数据点的均值。
python代码如下
import numpy as np
import matplotlib.pyplot as pltx=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.Tn=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&#61;distEclud,creatCen&#61;Randcenter):m&#61;dataset.shape[0]n&#61;dataset.shape[1]clusterAssment &#61; np.mat(np.zeros((m,1)))centoids&#61;creatCen(dataset,k)flag&#61;Truecount&#61;0while flag:count&#43;&#61;1print(&#39;第%d次迭代&#39;%count)flag&#61;Falseclusternum&#61;np.mat(np.zeros((k,1)))newcentoids&#61;np.mat(np.zeros((k,n)))for i in range(m):minDist&#61;np.infminIndex&#61;-1for j in range(k):temp_dist&#61;distMeans(centoids[j],dataset[i])if temp_dist<minDist:minDist&#61;temp_distminIndex&#61;jclusterAssment[i]&#61;minIndexnewcentoids[minIndex]&#43;&#61;dataset[i]clusternum[minIndex]&#43;&#61;1print(centoids)for cnt in range(k):newcentoids[cnt]/&#61;clusternum[cnt]if distMeans(centoids[cnt],newcentoids[cnt])>1e-10:flag&#61;Truecentoids&#61;newcentoidsreturn newcentoids,clusterAssmentcenters,asse&#61;KMeans(z,4)