用Python實現PCA和MDA降維和聚類

藥山浪子 2020-09-29 09:48:22



數據挖掘入門與實戰 ?公眾號: datadw


降維和聚類算是無監督學習的重要領域,還是那句話,不論是PCA、MDA還是K-means聚類,網上大??偨Y的杠杠的,給幾個參考鏈接:?


http://www.cnblogs.com/jerrylead/archive/2011/04/18/2020209.html?
http://bbezxcy.iteye.com/blog/2090591?
http://www.tuicool.com/articles/7nIvum?
http://www.cnblogs.com/python27/p/MachineLearningWeek08.html?
http://blog.pluskid.org/?p=407?
http://www.cnblogs.com/Key-Ky/archive/2013/11/24/3440684.html?
http://www.cnblogs.com/coser/archive/2013/04/10/3013044.html


PCA和MDA的推導過程都是手推,本來想拍照發上來,但前幾次‘作’過之后實在提不起興趣,還好有小伙伴(妹子)總結的很好:?


http://blog.csdn.net/totodum/article/details/51049165?
http://blog.csdn.net/totodum/article/details/51097329


來看我們這次課的任務:

?數據Cat4D3Groups是4維觀察數據,?
?請先采用MDS方法降維到3D,形成Cat3D3Groups數據,顯示并觀察。?
?對Cat3D3Groups數據采用線性PCA方法降維到2D,形成Cat2D3Groups數據,顯示并觀察。?


?對Cat2D3Groups數據采用K-Mean方法對數據進行分類并最終確定K,顯示分類結果。?
?對Cat2D3Groups數據采用Hierarchical分類法對數據進行分類,并顯示分類結果。

理論一旦推導完成,代碼寫起來就很輕松:

Part 1:降維處理?
MDA:

# -*- coding:gb2312 -*-from pylab import *import numpy as npfrom mpl_toolkits.mplot3d import Axes3Ddef print_D(data):
 ? ?N = np.shape(data)[0]
 ? ?d = np.zeros((N, N)) ? ?for i in range(N):
 ? ? ? ?c = data[i, :] ? ? ? ?for j in range(N):
 ? ? ? ? ? ?e = data[j, :]
 ? ? ? ? ? ?d[i, j] = np.sqrt(np.sum(np.power(c - e, 2))) ? ?return ddef MDS(D, K):
 ? ?N = np.shape(D)[0]
 ? ?D2 = D ** 2
 ? ?H = np.eye(N) - 1.0/N
 ? ?T = -0.5 * np.dot(np.dot(H, D2), H)
 ? ?eigVal, eigVec = np.linalg.eig(T)
 ? ?indices = np.argsort(eigVal) # 返回從小到大的索引值
 ? ?indices = indices[::-1] # 反轉

 ? ?eigVal = eigVal[indices] # 特征值從大到小排列
 ? ?eigVec = eigVec[:, indices] # 排列對應特征向量

 ? ?m = eigVec[:, :K]
 ? ?n = np.diag(np.sqrt(eigVal[:K]))
 ? ?X = np.dot(m, n) ? ?return X# test'''
data = genfromtxt("CAT4D3GROUPS.txt")
D = print_D(data)
# print D

# 4D 轉 3D
CAT3D3GROUPS = MDS(D, 3)
# print CAT3D3GROUPS
# D_3D = print_D(CAT3D3GROUPS)
# print D_3D
figure(1)
ax = subplot(111,projection='3d')
ax.scatter(CAT3D3GROUPS[:, 0], CAT3D3GROUPS[:, 1], CAT3D3GROUPS[:, 2], c = 'b')
ax.set_zlabel('Z') #坐標軸
ax.set_ylabel('Y')
ax.set_xlabel('X')
title('MDS_4to3')

# 4D 轉 2D
CAT2D3GROUPS = MDS(D, 2)
# print CAT2D3GROUPS
# D_2D = print_D(CAT2D3GROUPS)
# print D_2D
figure(2)
plot(CAT2D3GROUPS[:, 0], CAT2D3GROUPS[:, 1], 'b.')
xlabel('x')
ylabel('y')
title('MDS_4to2')
'''


PDA:

# -*- coding:gb2312 -*-from pylab import *from numpy import *from mpl_toolkits.mplot3d import Axes3Ddef PCA(data, K):
 ? ?# 數據標準化
 ? ?m = mean(data, axis=0) # 每列均值
 ? ?data -= m ? ?# 協方差矩陣
 ? ?C = cov(transpose(data)) ? ?# 計算特征值特征向量,按降序排序
 ? ?evals, evecs = linalg.eig(C)
 ? ?indices = argsort(evals) # 返回從小到大的索引值
 ? ?indices = indices[::-1] # 反轉

 ? ?evals = evals[indices] # 特征值從大到小排列
 ? ?evecs = evecs[:, indices] # 排列對應特征向量
 ? ?evecs_K_max = evecs[:, :K] # 取最大的前K個特征值對應的特征向量

 ? ?# 產生新的數據矩陣
 ? ?finaldata = dot(data, evecs_K_max) ? ?return finaldata# test'''
data = genfromtxt("CAT4D3GROUPS.txt")

# 4D 轉 3D
data_PCA = PCA(data, 3)
# print data_PCA
figure(1)
ax = subplot(111, projection='3d')
ax.scatter(data_PCA[:, 0], data_PCA[:, 1], data_PCA[:, 2], c='b')
ax.set_zlabel('Z') #坐標軸
ax.set_ylabel('Y')
ax.set_xlabel('X')
title('PCA_4to3')

# 4D 轉 2D
data_PCA = PCA(data, 2)
print data_PCA
figure(2)
plot(data_PCA[:, 0], data_PCA[:, 1], 'b.')
xlabel('x')
ylabel('y')
title('PCA_4to2')
'''

代碼里的注釋啰啰嗦嗦應該解釋的很清楚,這里不再贅述,看結果:

1、用MDS方法降維到3D,形成Cat3D3Groups數據:?
共兩個函數,輔助函數用來生成歐氏距離矩陣,MDS函數用于降維。?


通過輸出的距離矩陣可以看出,降維前后歐氏距離誤差小于10^-4,證明算法有效。同時旋轉3D圖像也可以明顯找出2D平面圖的視角


2、用PCA方法降維到2D,形成Cat2D3Groups數據:?


用PCA直接對4D數據降維后的結果與MDS等價,證明算法有效。同時旋轉3D圖像也可以明顯找出2D平面圖的視角。

3、總結分析:?
先用MDS算法將4D數據降到3D,再用PCA降到2D。?

與MDS降維生成的2D圖像及數據對比,誤差忽略不計,證明算法有效,同時證明MDS和PCA算法在進行小批量數據降維處理上效果類似。


Part 2:聚類分析:?
數據用前面降維之后的二維數據。K-means聚類分析的程序主要參考《Machine Learning in Action》- Peter Harrington這本書第十章,我自己添加了選擇最優K值的功能:

# -*- coding:gb2312 -*-import numpy as npfrom pylab import *from numpy import *# 求歐氏距離def euclDistance(vector1, vector2):
 ? ?return np.sqrt(np.sum(np.power(vector2 - vector1, 2)))# 返回某個值在列表中的全部索引值def myfind(x, y):
 ? ?return [ a for a in range(len(y)) if y[a] == x]# 初始化聚類點def initCentroids(data, k):
 ? ?numSamples, dim = data.shape
 ? ?centroids = np.zeros((k, dim)) ? ?for i in range(k):
 ? ? ? ?index = int(np.random.uniform(0, numSamples))
 ? ? ? ?centroids[i, :] = data[index, :] ? ?return centroids# K-mean 聚類def K_mean(data, k):
 ? ?## step 1: 初始化聚點
 ? ?centroides = initCentroids(data, k)
 ? ?numSamples = data.shape[0]
 ? ?clusterAssment = np.zeros((numSamples, 2)) # 保存每個樣本點的簇索引值和誤差
 ? ?clusterChanged = True

 ? ?while clusterChanged:
 ? ? ? ?clusterChanged = False
 ? ? ? ?global sum
 ? ? ? ?sum = [] ? ? ? ?# 對每一個樣本點
 ? ? ? ?for i in xrange(numSamples):
 ? ? ? ? ? ?minDist = np.inf # 記錄最近距離
 ? ? ? ? ? ?minIndex = 0 # 記錄聚點

 ? ? ? ? ? ?## step 2: 找到距離最近的聚點
 ? ? ? ? ? ?for j in range(k):
 ? ? ? ? ? ? ? ?distance = euclDistance(centroides[j, :], data[i, :]) ? ? ? ? ? ? ? ?if distance < minDist:
 ? ? ? ? ? ? ? ? ? ?minDist ?= distance
 ? ? ? ? ? ? ? ? ? ?minIndex = j ? ? ? ? ? ?## step 3: 將該樣本歸到該簇
 ? ? ? ? ? ?if clusterAssment[i, 0] != minIndex:
 ? ? ? ? ? ? ? ?clusterChanged = True # 前后分類相同時停止循環
 ? ? ? ? ? ? ? ?clusterAssment[i, :] = minIndex, minDist ** 2 # 記錄簇索引值和誤差


 ? ? ? ?## step 4: 更新聚點
 ? ? ? ?for j in range(k):
 ? ? ? ? ? ?index = myfind(j, clusterAssment[:, 0])
 ? ? ? ? ? ?pointsInCluster = data[index, :] # 返回屬于j簇的data中非零樣本的目錄值,并取出樣本

 ? ? ? ? ? ?centroides[j, :] = np.mean(pointsInCluster, axis=0) # 求列平均
 ? ? ? ? ? ?# 返回cost funktion值
 ? ? ? ? ? ?suml = 0
 ? ? ? ? ? ?lenght = pointsInCluster.shape[0] ? ? ? ? ? ?for l in range(lenght):
 ? ? ? ? ? ? ? ?dis = euclDistance(centroides[j, :], pointsInCluster[l, :])
 ? ? ? ? ? ? ? ?suml += dis ** 2 / lenght
 ? ? ? ? ? ?sum.append(suml)
 ? ?cost = np.sum(sum) / k ? ?print cost ? ?return centroides, clusterAssment, cost# 畫出分類前后結果def showCluster(data, k, centroides, clusterAssment):
 ? ?numSamples, dim = data.shape
 ? ?mark = ['r.', 'b.', 'g.', 'k.', '^r', '+r', 'sr', 'dr', '<r', 'pr'] ? ?# draw all samples
 ? ?for i in xrange(numSamples):
 ? ? ? ?markIndex = int(clusterAssment[i, 0])
 ? ? ? ?figure(2)
 ? ? ? ?plt.plot(data[i, 0], data[i, 1], mark[markIndex])
 ? ? ? ?plt.title('K-means')
 ? ?mark = ['Dr', 'Db', 'Dg', 'Dk', '^b', '+b', 'sb', 'db', '<b', 'pb'] ? ?# draw the centroids
 ? ?for i in range(k):
 ? ? ? ?figure(2)
 ? ? ? ?plt.plot(centroides[i, 0], centroides[i, 1], mark[i], markersize = 12)

其中三個輔助函數用于求歐氏距離,返回矩陣索引值和畫圖,k-mean函數用于聚類,當所有樣本點到其所屬聚類中心距離不變時,輸出聚類結果,并返回cost function的值。


Cost function計算方法:對每個簇,求所有點到所屬聚類中心的歐氏距離,平方后取均值E。聚類結束后,所有簇E值求和取平均得到cost function的值。

不同K值下的分類結果如下(標明聚類中心):?




主觀判斷,k = 4時聚類結果最優。用Elbow方法選擇K值結果如下:?

發現在K = 2時cost function值下降最為明顯,與之前判斷的結果不符。思考后發現,K=1時聚類沒有意義,所以上圖并不能有效選擇K值,調整后結果如下:?

明顯看出,k = 4時cost function下降極為明顯。與主觀判斷結果相符。

Hierarchical分類,參考網上代碼,出處不記得了:

# -*- coding:gb2312 -*-import numpy as npimport matplotlib.pyplot as pltimport MDSimport PCAdef yezi(clust):
 ? ?if clust.left == None and clust.right == None : ? ? ? ?return [clust.id] ? ?return yezi(clust.left) + yezi(clust.right)def Euclidean_distance(vector1,vector2):
 ? ?length = len(vector1)
 ? ?TSum = sum([pow((vector1[i] - vector2[i]),2) for i in range(len(vector1))])
 ? ?SSum = np.sqrt(TSum) ? ?return SSumdef loadDataSet(fileName):
 ? ?a = [] ? ?with open(fileName, 'r') as f:
 ? ? ? ?data = f.readlines() ?#txt中所有字符串讀入data
 ? ? ? ?for line in data:
 ? ? ? ? ? ?odom = line.split() ? ? ? ?#將單個數據分隔開存好
 ? ? ? ? ? ?numbers_float = map(float, odom) #轉化為浮點數
 ? ? ? ? ? ?a.append(numbers_float) #print numbers_float
 ? ?a = np.array(a) ? ?return aclass bicluster:
 ? ?def __init__(self, vec, left=None,right=None,distance=0.0,id=None):
 ? ? ? ?self.left = left
 ? ? ? ?self.right = right ?#每次聚類都是一對數據,left保存其中一個數據,right保存另一個
 ? ? ? ?self.vec = vec ? ? ?#保存兩個數據聚類后形成新的中心
 ? ? ? ?self.id = id
 ? ? ? ?self.distance = distancedef list_array(wd, clo):
 ? ?len_=len(wd)
 ? ?xc=np.zeros([len_, clo]) ? ?for i in range(len_):
 ? ? ? ?ad = wd[i]
 ? ? ? ?xc[i, :] = ad ? ?return xcdef hcluster(data, n) :
 ? ?[row,column] = data.shape
 ? ?data = list_array(data, column)

 ? ?biclusters = [bicluster(vec = data[i], id = i) for i in range(len(data))]
 ? ?distances = {}
 ? ?flag = None
 ? ?currentclusted = -1
 ? ?while(len(biclusters) > n) : #假設聚成n個類
 ? ? ? ?min_val = np.inf ? #Python的無窮大
 ? ? ? ?biclusters_len = len(biclusters) ? ? ? ?for i in range(biclusters_len-1) : ? ? ? ? ? ?for j in range(i + 1, biclusters_len): ? ? ? ? ? ? ? ?#print biclusters[i].vec

 ? ? ? ? ? ? ? ?if distances.get((biclusters[i].id,biclusters[j].id)) == None: ? ? ? ? ? ? ? ? ? ?#print biclusters[i].vec
 ? ? ? ? ? ? ? ? ? ?distances[(biclusters[i].id,biclusters[j].id)] = Euclidean_distance(biclusters[i].vec,biclusters[j].vec)
 ? ? ? ? ? ? ? ?d = distances[(biclusters[i].id,biclusters[j].id)] ? ? ? ? ? ? ? ?if d < min_val:
 ? ? ? ? ? ? ? ? ? ?min_val = d
 ? ? ? ? ? ? ? ? ? ?flag = (i,j)
 ? ? ? ?bic1,bic2 = flag #解包bic1 = i , bic2 = j
 ? ? ? ?newvec = [(biclusters[bic1].vec[i] + biclusters[bic2].vec[i])/2 for i in range(len(biclusters[bic1].vec))] #形成新的類中心,平均
 ? ? ? ?newbic = bicluster(newvec, left=biclusters[bic1], right=biclusters[bic2], distance=min_val, id = currentclusted) #二合一
 ? ? ? ?currentclusted -= 1
 ? ? ? ?del biclusters[bic2] #刪除聚成一起的兩個數據,由于這兩個數據要聚成一起
 ? ? ? ?del biclusters[bic1]
 ? ? ? ?biclusters.append(newbic)#補回新聚類中心
 ? ? ? ?clusters = [yezi(biclusters[i]) for i in range(len(biclusters))] #深度優先搜索葉子節點,用于輸出顯示
 ? ?return biclusters,clustersdef showCluster(dataSet, k, num_mark):
 ? ?numSamples, dim = dataSet.shape
 ? ?mark = ['r.', 'b.', 'g.', 'k.', '^r', '+r', 'sr', 'dr', '<r', 'pr'] ? ?# draw all samples
 ? ?for i in xrange(numSamples):
 ? ? ? ?plt.plot(dataSet[i, 0], dataSet[i, 1], mark[num_mark])

 ? ?plt.xlabel('X')
 ? ?plt.ylabel('Y')
 ? ?plt.title('Hierarchical')if __name__ == "__main__": ? ?# 加載數據
 ? ?dataMat =np.genfromtxt('CAT4D3GROUPS.txt') ? ?#400*4
 ? ?dataSet = PCA.PCA(dataMat, 2)

 ? ?k,l = hcluster(dataSet, 10) ?# ?l返回了聚類的索引
 ? ?# 選取規模最大的k個簇,其他簇歸為噪音點
 ? ?for j in range(len(l)):
 ? ? ? ?m = [] ? ? ? ?for ii in range(len(l[j])):
 ? ? ? ? ? ?m.append(l[j][ii])
 ? ? ? ?m = np.array(m)
 ? ? ? ?a = dataSet[m]
 ? ? ? ?showCluster(a,len(l),j)
 ? ?plt.show()


當一個類集合中包含多個樣本點時,類與類之間的距離取Group Average:把兩個集合中的點兩兩的歐氏距離全部放在一起求平均值,分類結果如下:?



重復運行后分類結果并未有太大變化。主觀判斷,從分成3類及4類的結果看,Hierarchical分類方法效果不如K-mean聚類效果好。


數據挖掘入門與實戰

搜索添加微信公眾號:datadw


教你機器學習,教你數據挖掘


長按圖片,識別二維碼,點關注



公眾號推薦: weic2c? ?


數據分析入門與實戰

從哪里做起學習數據分析?

如何培養數據分析的能力?



長按圖片,識別二維碼,點關注


發表評論
用戶反饋
客戶端