使用python實現kmean演演算法

2023-11-02 22:00:32

1. 簡介

kmean 是無監督學習的一種演演算法，主要是用來進行聚類分析的，他會在資料集中算出幾個點作為簇中心，求這些資料集與這些簇中心的距離，並將距離同一個簇中心距離最近的資料歸為一類。因此，kmean最重要的地方便是關於簇中心的選擇。

他的演演算法流程簡單總結如下

簇個數的選擇；
計算樣本到選取的簇中心距離，劃分樣本，將距離同一個簇中心最近的樣本歸為一類；
設定一個迭代次數，不斷更新簇中心；

2. kmean演演算法過程

這裡初始簇中心的選擇並不是隨機選擇因為隨機選擇的話，有的簇中心可能會重疊，這裡改為了先選擇距離最遠的兩個樣本作為兩個初始簇中心，下一個簇中心就從距離之前的簇中心距離最近的樣本中將距離最遠的作為下一個，如此反覆直到滿足定義的k大小。同時圖中的新的簇中心的更新應該是新老簇中心不相等時才更新才對。其實kmean的效果好壞很大的程度上取決於簇中心數量和初始簇中心的選擇上，並且他適合一些聚整合簇的樣本進行劃分。

2.1 簇個數的選擇

畫圖法：最直接的辦法還是畫圖，從圖中直接看整個樣本大概分為幾大類，這種辦法適合低維（1-3維）的資料，對於高維資料很難直接成圖，不過可以通過特徵提取選擇能體現樣本差異的特徵維（這裡的差異很難定義，畢竟是無監督學習，可以簡單認為成圖後樣本聚集為幾類，這些類間間隔較遠），或者直接降維為低維度資料成圖觀察。
窮舉k法：通過定義一系列不同的簇個數，並通過聚類指標評價這些不同簇個數的效果，選擇效果最好的那個；

2.2 聚類評價指標

還有更多的評價指標，這裡我選擇了一些簡單的，易於實現的。

2.2.1 輪廓係數

參考：https://baike.baidu.com/item/%E8%BD%AE%E5%BB%93%E7%B3%BB%E6%95%B0/17361607?fr=aladdin

補充：如果Si接近0代表該樣本越可能在類的邊界上；越接近-1代表越可能屬於其他類；

2.2.2 緊密性指標

這裡指的是先計算同簇的樣本距離簇中心的平均距離CPi，再計算CPi的平均值；

2.2.3 間隔性指標

3. 程式碼

資料集：經典的鳶尾花資料集，裡面收錄了三個不同品種共四個不同性狀的資料。

#資料集
5.1,3.5,1.4,0.2,Iris-setosa
4.9,3.0,1.4,0.2,Iris-setosa
4.7,3.2,1.3,0.2,Iris-setosa
4.6,3.1,1.5,0.2,Iris-setosa
5.0,3.6,1.4,0.2,Iris-setosa
5.4,3.9,1.7,0.4,Iris-setosa
4.6,3.4,1.4,0.3,Iris-setosa
5.0,3.4,1.5,0.2,Iris-setosa
4.4,2.9,1.4,0.2,Iris-setosa
4.9,3.1,1.5,0.1,Iris-setosa
5.4,3.7,1.5,0.2,Iris-setosa
4.8,3.4,1.6,0.2,Iris-setosa
4.8,3.0,1.4,0.1,Iris-setosa
4.3,3.0,1.1,0.1,Iris-setosa
5.8,4.0,1.2,0.2,Iris-setosa
5.7,4.4,1.5,0.4,Iris-setosa
5.4,3.9,1.3,0.4,Iris-setosa
5.1,3.5,1.4,0.3,Iris-setosa
5.7,3.8,1.7,0.3,Iris-setosa
5.1,3.8,1.5,0.3,Iris-setosa
5.4,3.4,1.7,0.2,Iris-setosa
5.1,3.7,1.5,0.4,Iris-setosa
4.6,3.6,1.0,0.2,Iris-setosa
5.1,3.3,1.7,0.5,Iris-setosa
4.8,3.4,1.9,0.2,Iris-setosa
5.0,3.0,1.6,0.2,Iris-setosa
5.0,3.4,1.6,0.4,Iris-setosa
5.2,3.5,1.5,0.2,Iris-setosa
5.2,3.4,1.4,0.2,Iris-setosa
4.7,3.2,1.6,0.2,Iris-setosa
4.8,3.1,1.6,0.2,Iris-setosa
5.4,3.4,1.5,0.4,Iris-setosa
5.2,4.1,1.5,0.1,Iris-setosa
5.5,4.2,1.4,0.2,Iris-setosa
4.9,3.1,1.5,0.1,Iris-setosa
5.0,3.2,1.2,0.2,Iris-setosa
5.5,3.5,1.3,0.2,Iris-setosa
4.9,3.1,1.5,0.1,Iris-setosa
4.4,3.0,1.3,0.2,Iris-setosa
5.1,3.4,1.5,0.2,Iris-setosa
5.0,3.5,1.3,0.3,Iris-setosa
4.5,2.3,1.3,0.3,Iris-setosa
4.4,3.2,1.3,0.2,Iris-setosa
5.0,3.5,1.6,0.6,Iris-setosa
5.1,3.8,1.9,0.4,Iris-setosa
4.8,3.0,1.4,0.3,Iris-setosa
5.1,3.8,1.6,0.2,Iris-setosa
4.6,3.2,1.4,0.2,Iris-setosa
5.3,3.7,1.5,0.2,Iris-setosa
5.0,3.3,1.4,0.2,Iris-setosa
7.0,3.2,4.7,1.4,Iris-versicolor
6.4,3.2,4.5,1.5,Iris-versicolor
6.9,3.1,4.9,1.5,Iris-versicolor
5.5,2.3,4.0,1.3,Iris-versicolor
6.5,2.8,4.6,1.5,Iris-versicolor
5.7,2.8,4.5,1.3,Iris-versicolor
6.3,3.3,4.7,1.6,Iris-versicolor
4.9,2.4,3.3,1.0,Iris-versicolor
6.6,2.9,4.6,1.3,Iris-versicolor
5.2,2.7,3.9,1.4,Iris-versicolor
5.0,2.0,3.5,1.0,Iris-versicolor
5.9,3.0,4.2,1.5,Iris-versicolor
6.0,2.2,4.0,1.0,Iris-versicolor
6.1,2.9,4.7,1.4,Iris-versicolor
5.6,2.9,3.6,1.3,Iris-versicolor
6.7,3.1,4.4,1.4,Iris-versicolor
5.6,3.0,4.5,1.5,Iris-versicolor
5.8,2.7,4.1,1.0,Iris-versicolor
6.2,2.2,4.5,1.5,Iris-versicolor
5.6,2.5,3.9,1.1,Iris-versicolor
5.9,3.2,4.8,1.8,Iris-versicolor
6.1,2.8,4.0,1.3,Iris-versicolor
6.3,2.5,4.9,1.5,Iris-versicolor
6.1,2.8,4.7,1.2,Iris-versicolor
6.4,2.9,4.3,1.3,Iris-versicolor
6.6,3.0,4.4,1.4,Iris-versicolor
6.8,2.8,4.8,1.4,Iris-versicolor
6.7,3.0,5.0,1.7,Iris-versicolor
6.0,2.9,4.5,1.5,Iris-versicolor
5.7,2.6,3.5,1.0,Iris-versicolor
5.5,2.4,3.8,1.1,Iris-versicolor
5.5,2.4,3.7,1.0,Iris-versicolor
5.8,2.7,3.9,1.2,Iris-versicolor
6.0,2.7,5.1,1.6,Iris-versicolor
5.4,3.0,4.5,1.5,Iris-versicolor
6.0,3.4,4.5,1.6,Iris-versicolor
6.7,3.1,4.7,1.5,Iris-versicolor
6.3,2.3,4.4,1.3,Iris-versicolor
5.6,3.0,4.1,1.3,Iris-versicolor
5.5,2.5,4.0,1.3,Iris-versicolor
5.5,2.6,4.4,1.2,Iris-versicolor
6.1,3.0,4.6,1.4,Iris-versicolor
5.8,2.6,4.0,1.2,Iris-versicolor
5.0,2.3,3.3,1.0,Iris-versicolor
5.6,2.7,4.2,1.3,Iris-versicolor
5.7,3.0,4.2,1.2,Iris-versicolor
5.7,2.9,4.2,1.3,Iris-versicolor
6.2,2.9,4.3,1.3,Iris-versicolor
5.1,2.5,3.0,1.1,Iris-versicolor
5.7,2.8,4.1,1.3,Iris-versicolor
6.3,3.3,6.0,2.5,Iris-virginica
5.8,2.7,5.1,1.9,Iris-virginica
7.1,3.0,5.9,2.1,Iris-virginica
6.3,2.9,5.6,1.8,Iris-virginica
6.5,3.0,5.8,2.2,Iris-virginica
7.6,3.0,6.6,2.1,Iris-virginica
4.9,2.5,4.5,1.7,Iris-virginica
7.3,2.9,6.3,1.8,Iris-virginica
6.7,2.5,5.8,1.8,Iris-virginica
7.2,3.6,6.1,2.5,Iris-virginica
6.5,3.2,5.1,2.0,Iris-virginica
6.4,2.7,5.3,1.9,Iris-virginica
6.8,3.0,5.5,2.1,Iris-virginica
5.7,2.5,5.0,2.0,Iris-virginica
5.8,2.8,5.1,2.4,Iris-virginica
6.4,3.2,5.3,2.3,Iris-virginica
6.5,3.0,5.5,1.8,Iris-virginica
7.7,3.8,6.7,2.2,Iris-virginica
7.7,2.6,6.9,2.3,Iris-virginica
6.0,2.2,5.0,1.5,Iris-virginica
6.9,3.2,5.7,2.3,Iris-virginica
5.6,2.8,4.9,2.0,Iris-virginica
7.7,2.8,6.7,2.0,Iris-virginica
6.3,2.7,4.9,1.8,Iris-virginica
6.7,3.3,5.7,2.1,Iris-virginica
7.2,3.2,6.0,1.8,Iris-virginica
6.2,2.8,4.8,1.8,Iris-virginica
6.1,3.0,4.9,1.8,Iris-virginica
6.4,2.8,5.6,2.1,Iris-virginica
7.2,3.0,5.8,1.6,Iris-virginica
7.4,2.8,6.1,1.9,Iris-virginica
7.9,3.8,6.4,2.0,Iris-virginica
6.4,2.8,5.6,2.2,Iris-virginica
6.3,2.8,5.1,1.5,Iris-virginica
6.1,2.6,5.6,1.4,Iris-virginica
7.7,3.0,6.1,2.3,Iris-virginica
6.3,3.4,5.6,2.4,Iris-virginica
6.4,3.1,5.5,1.8,Iris-virginica
6.0,3.0,4.8,1.8,Iris-virginica
6.9,3.1,5.4,2.1,Iris-virginica
6.7,3.1,5.6,2.4,Iris-virginica
6.9,3.1,5.1,2.3,Iris-virginica
5.8,2.7,5.1,1.9,Iris-virginica
6.8,3.2,5.9,2.3,Iris-virginica
6.7,3.3,5.7,2.5,Iris-virginica
6.7,3.0,5.2,2.3,Iris-virginica
6.3,2.5,5.0,1.9,Iris-virginica
6.5,3.0,5.2,2.0,Iris-virginica
6.2,3.4,5.4,2.3,Iris-virginica
5.9,3.0,5.1,1.8,Iris-virginica

#python3

#----------------------------
# author: little shark
# date: 2022/3/12
"""
Kmean 實現
"""

from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
import random

class Algo:
	def __init__(self):
		self.x = None
		self.y = None
		self.rows,self.cols = None, None
		self.dis = None
		self.dis_i = None
		self.rbf_i = None
		self.G_i = None
		self.GDis_i = True
		
	def loadData(self,file):
		"""
		Read file
		"""
		df=pd.read_csv(file,delimiter=",",header=None)
		self.x=df.iloc[:,0:-1].values
		self.y=df.iloc[:,-1].values
		
		self.rows,self.cols=self.x.shape
		#print(self.x.shape)
		
	def dataNorm(self):
		"""
		歸一化資料
		"""
		assert self.rows != None or self.cols !=None, "First you must loadData!!!"
		
		#get min matrix
		min_=self.x.min(axis=0).reshape(1,self.cols)
		
		#get max matrix
		max_=self.x.max(axis=0).reshape(1,self.cols)
		
		#get range matrix
		range_=max_-min_

		#update self.x
		self.x=(self.x-min_)/range_
	
	def euDistances(self,x):
		"""
		Get EuclideanDistances matrix
		"""
		assert self.rows != None or self.cols !=None, "First you must loadData!!!"
		
		#輸入矩陣大小
		rows,cols = x.shape
		
		#define a empty matrix
		dis=np.zeros(rows**2).reshape(rows,rows)
		
		#fill dis matrix by calculating distance of each point
		for i in range(rows):
			for j in range(rows):
				if i == j:
					dis[i][j] = 0
					continue
				d=sum(pow(x[i]-x[j],2))**0.5
				dis[i][j]=d
				dis[j][i]=d
		
		#get EuclideanDistances matrix
		self.dis=dis.copy()
		self.dis_i=True
		
		return dis
		
		
	def kmean(self,k=3,maxIteration=100,x=None,print_=False):
		"""
		kmean實體
		"""
		assert self.rows != None or self.cols !=None, "First you must loadData!!!"
		assert x.any() != None, "you must import a data to kmean"
		
		#輸入矩陣大小
		rows,cols = x.shape
		
		#定義一個標籤列表，代表所有樣本的類別
		label=np.zeros(rows)
		labelCode=[i for i in range(k)]
		
		#計算輸入矩陣的距離矩陣
		dis=self.euDistances(x)
		
		#------------------------------------------
		#選擇初始簇中心，數量=k
		#最遠的兩個點當作簇中心
		max_=0
		maxI=0
		maxJ=0
		for i in range(rows):
			for j in range(rows):
				if dis[i][j] > max_:
					max_ = dis[i][j]
					maxI = i
					maxJ = j
		k1I = maxI
		k2I = maxJ
		k1=x[k1I].reshape(1,cols)
		k2=x[k2I].reshape(1,cols)
		
		Ci=[k1I,k2I] #儲存簇中心索引
		C=[k1,k2] #儲存簇中心值
		
		#step 3,檢索其他簇中心
		for i in range(2,k):
			tempI=[]
			temp=[]
			for j in range(len(C)):
				
				#step 3.1 距離簇中心最近的點
				minI=0
				min_=99999
				index=0
				for d in dis[Ci[j]]:
					#略過本身和已經檢索過的點
					if d==0 or index in Ci:
						index+=1
						continue
					if d<min_:
						minI=index
						min_=d
					index+=1
				
				
				#step 3.2 新增距離大小及其索引
				temp.append(dis[Ci[j]][minI])
				tempI.append(minI)
				
			#step 3.3 選擇距離已選出的簇中心距離最近的點距離最大的作為下一個簇中心
			nextKI=tempI[temp.index(max(temp))]
			nextK=x[nextKI].reshape(1,cols)
			Ci.append(nextKI)
			C.append(nextK)
		
		#------------------------------------
		#更新label
		for i in range(rows):
			sample=x[i].reshape(1,cols)
			
			#判斷距離最近的簇中心
			minI=0
			min_=99999
			for j in range(k):
				disV=self.getDistance(sample,C[j])
				#print(dis)
				if disV<min_:
					minI=j
					min_=disV
			label[i]=labelCode[minI]
		
		#------------------------------------
		#更新簇中心
		for iter in range(maxIteration):
			for j in range(k):
				#獲取同類別下的資料
				dataSet=x[label == labelCode[j]]
				if len(dataSet) == 0:
					continue
				
				#計算平均值
				mean_=dataSet.mean(axis=0).reshape(1,cols)
				
				#更新簇中心
				C[j]=mean_
				
			#------------------------------------
			#更新label
			for i in range(rows):
				sample=x[i].reshape(1,cols)
				
				#判斷距離最近的簇中心
				minI=0
				min_=99999
				for j in range(k):
					disV=self.getDistance(sample,C[j])
					#print(dis)
					if disV<min_:
						minI=j
						min_=disV
				label[i]=labelCode[minI]
			
			cp = self.CP(label,C,x)
			sp = self.SP(C)
			conCoe = self.contourCoefficient(label,C,x)
			
			if print_:
				print("Iteration %-4d CP = %-7.4f CP = %-7.4f coef = %-7.4f"%(iter,cp,sp,conCoe))
			
		return label,C
	
	def getDistance(self,x1,x2):
		"""
		計算兩個向量的幾何距離
		"""
		return pow(((x1-x2)**2).sum(),0.5)
				
	def CP(self,label,c,x):
		"""
		緊密性評價指標
		"""
		k = len(c)
		labelCode = np.unique(label)
		
		cp = 0
		for i in range(k):
			#同簇矩陣
			x1 = x[label == i]
			
			if len(x1) == 0:continue
			
			#大小
			rows, cols = x1.shape
			
			#距離簇中心的歐式距離和
			sum_ = pow(((x1 - c[i])**2).sum(),0.5)
			
			cp += sum_/rows
		
		return cp/k
		
	def SP(self,c):
		"""
		間隔性評價指標
		"""
		k = len(c)
		
		sp = 0
		for i in range(k):
			for j in range(1,k):
				sp += self.getDistance(c[i],c[j])
		
		return 2*sp/(k**2-k)
		
	
	def calAi(self,v,vs):
		"""
		輪廓係數中計算ai
		v: 某向量
		vs: 該向量的同簇向量
		"""
		rows,cols = vs.shape
		ai = pow((vs - v)**2,0.5).sum()/rows
		
		return ai
		
	def calBi(self,v,vs,vsLabel):
		"""
		輪廓係數中計算bi
		v: 某向量
		vs: 該向量的非同簇向量
		vsCode: 該向量的非同簇向量的label
		"""
		ks = np.unique(vsLabel)
		k = len(ks)
		min_ = 99999
		
		for i in range(k):
			#非同簇中同簇向量
			x1 = vs[vsLabel == i]
			
			if len(x1) == 0:continue
			
			#非同簇中同簇平均距離
			disAver = self.calAi(v,x1)
			
			if disAver < min_:
				min_ = disAver
		
		return min_
		
	def contourCoefficient(self,label,c,x):
		"""
		輪廓係數
		"""
		k = len(c)
		rows,cols = x.shape
		S = np.zeros(rows)
		
		for i in range(rows):
			#同簇樣本
			x1 = x[label == label[i]]
			
			#非同簇樣本
			x2 = x[label != label[i]]
			label2 = label[label != label[i]]
			
			#當前樣本
			sample = x[i]
			
			#ai ,bi 計算
			ai = self.calAi(sample,x1)
			bi = self.calBi(sample,x2,label2)
			
			s = (bi-ai)/max([ai,bi])
			
			S[i] = s
		
		return S.mean()
		
if __name__ == "__main__":
	#-------------------------
	#直接定義簇個數
	a=Algo()
	
	#載入資料
	a.loadData("資料集/iris (2).data")
	
	#執行演演算法
	label,c = a.kmean(x=a.x,k=3)
	kCode = np.unique(label)
	
	print(label)
	
	#畫出資料劃分類
	for k in kCode:
		plt.scatter(a.x[label == k][:,0],a.x[label == k][:,2],alpha=0.6,label="sample %d"%k)
	
	#畫出簇中心
	for i in range(len(c)):
		plt.scatter(x=c[i][:,0],y=c[i][:,2],color="red",marker="*",label="center %d"%i)
	plt.legend()
	plt.show()
	
	
	#----------------------------
	#窮舉k法
	#exit()
	coef = []
	CP = []
	SP = []
	for k in range(2,10):
		label,c = a.kmean(x=a.x,k=k)
		coef.append(a.contourCoefficient(label,c,a.x))
		CP.append(a.SP(c))
		SP.append(a.CP(label,c,a.x))
		print("* When k = %d done *"%k)
	
	plt.subplot(131)
	plt.plot(range(2,10),SP,alpha=0.8,linestyle="--",color="red",label="SP")
	plt.xlabel("K")
	plt.legend()
	
	plt.subplot(132)
	plt.plot(range(2,10),coef,alpha=0.8,linestyle=":",color="green",label="coef")
	plt.xlabel("K")
	plt.legend()
	
	plt.subplot(133)
	plt.plot(range(2,10),CP,alpha=0.8,linestyle="-",color="blue",label="CP")
	plt.xlabel("K")
	plt.legend()
	plt.show()

4. 輸出結果

4.1 命令列輸出

解釋：首先輸出的是所有樣本的聚類劃分，這裡並沒有打亂樣本的順序進行分析，所以同品種的樣本是連續的一大塊，大致來看直接定義k=3的劃分效果還不錯。後面的輸出是窮舉k法，並計算對應的評價指標程序。

4.2 視覺化

解釋：不同顏色代表不同類別的資料，紅色星星代表簇中心，這裡資料的維度展示的是第1列和第3列維度。

4.3 窮舉k法指標

解釋：第一個圖是間隔性指標SP用來衡量聚類中心兩兩之間的平均距離，隨著k值增加，有比較明顯的波動，可以看見當k=3時有一個最小值，表明k=3時類間的聚類最近，理想情況下我們希望這個這個指標比較大，表明類別區分明顯；第二個是輪廓係數，隨著k值增加逐漸下降；第三個是緊密性指標，用來衡量內中資料距離簇中心距離，隨著k值增加呈現一點的波動；從圖中來看任意一個指標都不能準確的表明k=3(真實的簇數)是一個好的簇數量，也許這些指標並不適合用來尋找定義的簇數量。

單從輪廓係數來看，在知道k=3的情況下，越接近1代表聚類越合理，顯然是不現實的，且從公式上來說當k=1時，輪廓係數會非常接近1，難道說“大家都是一家人，四海皆兄弟”嗎？不太可能吧，且k=2的時候係數為0.75，難道劃分為2更好嗎？顯然在知道真實分類情況下也不適合，那麼我們統計分析的結果真的就和真實一樣嗎？

首先從kmean上來說，他是無監督的演演算法，這些資料是沒有先驗的標籤的，他是用來衡量資料相似情況，並把相似度較高的資料歸為一類，體現的是”同種中屬性更接近的思想“，可是他受限於初始k的選擇上；再從分類指標上來看，這些指標主要是用來衡量簇內、簇間劃分程度，如果k越大，越多的樣本可能”各自為政“，那麼他們間的SP、CP和輪廓係數可知是越來越小，從這些指標代表的意義來說並不能準確的反映簇個數。

再者，人類認知上存在缺陷，我們認知我們所認知的，我們並不能直觀的認知到客觀規律，只能通過我們能感知到的一切來觀測進而總結推理出客觀規律，我們目前的知識都是”前人拿著測量工具和放大鏡“一點點積累下來的，我們所思所想的都是實踐的產物，工具、思想都是人類活動下的結晶，也就是說存在不足，可是卻是我們能拿出手的”家當“。

在這裡我主要想說的是在資料分析的時候不要太過迷信指標，特別是對一些不知真實分佈的資料來做分析，如果這些指標能準確並一致的反映真實資料，那我們早就稱霸宇宙，預測未來了。

回到之前的話題：我們統計分析的結果真的就和真實一樣嗎? 就目前來說，這些測量和分析方法是我們能拿出來的，隨著人類的進步這些東西要麼進化，要麼淘汰。這些工具只是在一定程度上反映了真實而已，也許離得很近，也許大相徑庭…至於最終的結果在沒有真正認知到客觀規律前，我們永遠不知道，只是拿著我們人類活動的一系列產物來標榜而已，也許未來人類回首現在，也會覺得有趣吧。

到此這篇關於使用python實現kmean演演算法的文章就介紹到這了,更多相關python實現kmean演演算法內容請搜尋it145.com以前的文章或繼續瀏覽下面的相關文章希望大家以後多多支援it145.com！