@sprt
*寫在開頭：博主在開始學習機器學習和Python之前從未有過任何程式設計經驗，這個系列寫在學習這個領域一個月之後，完全從一個入門級菜鳥的角度記錄我的學習歷程，程式碼未經優化，僅供參考。有錯誤之處歡迎大家指正。
系統：win7-CPU;
程式設計環境：Anaconda2-Python2.7，IDE：pycharm5;
參考書籍：
《Neural Networks and Learning Machines(Third Edition)》- Simon Haykin;
《Machine Learning in Action》- Peter Harrington;
《Building Machine Learning Systems with Python》- Wili Richert;
C站裡都有資源，也有中文譯本。
我很慶幸能跟隨老師從最基礎的東西學起，進入機器學習的世界。*

來看我們這次課的任務：

•資料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 np
from mpl_toolkits.mplot3d import Axes3D

def 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 d

def 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 Axes3D

def 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 np
from 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 np
import matplotlib.pyplot as plt
import MDS
import PCA

def 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 SSum

def 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 a

class 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 = distance

def list_array(wd, clo):
    len_=len(wd)
    xc=np.zeros([len_, clo])
    for i in range(len_):
        ad = wd[i]
        xc[i, :] = ad
    return xc

def 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,clusters

def 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聚類效果好。