2017年2月23日 星期四

機器學習(6)--主成分分析(Principal component analysis,PCA)



    這一篇介紹主成分分析(Principal component analysis,PCA),是一種特徵提取的技術,利用特徵降維來避免因維度災難所造成的過度適合(Overfitting)現象。就是說假如訓練資料集裡有很多個特徵數據,我們可以利用主成分分析(Principal component analysis,PCA)來選取最有影響的幾個特徵來做分類器模型的訓練而不需要使用所有的特徵數據來做訓練。
    在這個範例裡,我們將使用葡萄酒數據集來做示範,把原本有的13個特徵,利用PCA分析後只找出兩個主要的特徵值作為分類器模型分類的特徵依據。

    PCA著眼在高維度的數據,"最大化變異數"並投影到與原數據及相同維數或較低維數的"新特徵子空間"。新特徵子空間的正交軸即是"主成分"。

<圖一>x1,x2軸維原本的,PC1,PC2為主成分







    使用PCA來做降維,需要一個dxk維的轉換矩陣W。如下列原本d維的數據x,對應到新的k維
特徵子空間z,這個新的特徵子空間z會比原本d維的特徵空間小,也就是k<d。

    從d維x特徵空間,轉換成新k維的z特徵子空間,將把多變量(multivariate)轉換成好幾個"主成分"(Principal component),並依變異數的大小排序,擁有最大變異數的"主成分"會排在最前面。另外,主成分彼此之間不相關,相關係數為0。
   在轉換為新特徵子空間前,要先做"特徵標準化",原因是PCA非常敏感於數據縮放,因此,如果原特徵是以不同尺度規範量度的話,那就一定要做"特徵標準化",讓所有的特徵都有相
同的尺規。


PCA演算法步驟:
1.標準化d維原"數據集"
2.建立共變異數矩陣(covariance matrix)
3.分解共變異數矩陣(covariance matrix)為特徵向量(eigenvector)與特徵值(eigenvalues)。
4.選取k個最大特徵值(eigenvalues)相對應k個的特徵向量(eigenvector),其中k即為新特徵子空間的維數。
5.使用排序最上面的k個的特徵向量(eigenvector),建立投影矩陣(project matrix)W。
6.使用投影矩陣(project matrix)W轉換原本d維的原數據至新的k維特徵子空間


葡萄酒數據集可從以下網址下載:
https://archive.ics.uci.edu/ml/datasets/Wine

1.標準化d維原"數據集"
    直接利用利用sci-kit learn StandardScaler 函數便可標準化數據
from sklearn.preprocessing import StandardScaler
stdsc = StandardScaler()
X_train_std = stdsc.fit_transform(X_train)
X_test_std = stdsc.fit_transform(X_test)

2. 共變異數矩陣
    在對稱dxd維共變異數矩陣中,儲存著不同特徵兩兩之間的變異數,
    兩個特徵xj與xk之間的共變異數,可用以下方式計算:


uj,uk分別為特徵j與k的平均值,這個平均值會等於0,因為我們已經過標準化處理。
而三個特徵的共變異數矩陣可用以下方式表示:Σ為sigma,勿與加總符號混淆:

    共變異數矩陣中的特徵向量代表主成分最大變異數的方向,而特徵值即代表其幅度(magnitude)。在葡萄酒數據集中,我們會從13x13維(d=13)的共變異數矩陣中,找到13個
特徵向量與特徵值。
         
3.分解共變異數矩陣(covariance matrix)為特徵向量(eigenvector)與特徵值(eigenvalues)。
    我們可以直接使用Numpy 中的linalg.eig函數來計算共變異數矩陣中的特徵對:

cov_mat = np.cov(X_train_std.T)                              #cov_mat 共變異數矩陣
eigen_vals, eigen_vecs = np.linalg.eig(cov_mat)    #eigen_vals為特徵值,eigen_vecs為 特徵向量

接著我們來看哪一主成分的"解釋變異數比率"比較高
特徵值λj的解釋變異數比率,可定義為特徵值λj對全特徵值總和的比例:

由下圖可以看出第一個主成分佔40%變異數,第二個主成分約佔20%,合併60%。





4,5 選定k=2,並建立投影矩陣(project matrix)W
    接著,以特徵值為依據以遞減的方式來排序,並選定特徵向量(eigenvector)建立投影矩陣(project matrix)W。在此選定13x2維的投影矩陣(project matrix)W。





6.使用投影矩陣(project matrix)W轉換原本d維的原數據至新的k維特徵子空間

     使用公式x'=xW,計算矩陣內積,便可將原葡萄酒數據集124x13維訓練數據集(在此範例分70%訓練數據集,30%為測試數據集),轉換到兩個主成分。將散布圖畫出如下:可以看出利用線性分類器,應該就可以分類的很成功。




    最後我們將此轉換後的訓練數據集以邏輯斯迴歸分類器來做分類,底下即為決策區域圖:
可以看出只有兩個分類錯誤。



接著對測試數據集做分類,可以看出只有一個分類錯誤。由於測試數據集也正確分類,所以並無過度適合(overfitting)的現象。


    整體簡單的說就是原本葡萄酒數據集有13個特徵值,經過PCA分析後我們只選定兩個主成分,將投影過後的數據利用邏輯斯迴歸便可分類適當,也不用擔心因使用太多特徵值維數災難造成的過度適合(Overfitting)。


<PCA範例程式:>

from sklearn.ensemble import RandomForestClassifier
import pandas as pd
import numpy as np
import warnings,math
import matplotlib.pyplot as plt

##讀取wine 資料+
df_wine = pd.read_csv('wine_data.csv',
       header=None)      
       
df_wine.columns = ['Class label', 'Alcohol', 'Malic acid', 'Ash',
       'Alcalinity of ash', 'Magnesium', 'Total phenols',
       'Flavanoids', 'Nonflavanoid phenols', 'Proanthocyanins',
       'Color intensity', 'Hue', 'OD280/OD315 of diluted wines',
       'Proline']
from sklearn.model_selection import train_test_split
X, y = df_wine.iloc[:, 1:].values, df_wine.iloc[:, 0].values #0:d為 class  1:之後為13個特徵行
X_train, X_test, y_train, y_test = \
 train_test_split(X, y, test_size=0.3, random_state=0)

from sklearn.preprocessing import StandardScaler
stdsc = StandardScaler()
X_train_std = stdsc.fit_transform(X_train)
X_test_std = stdsc.fit_transform(X_test)
##讀取wine 資料-

#求共變異係數矩陣
cov_mat = np.cov(X_train_std.T)
print("共變異係數矩陣.shape=",cov_mat.shape)
print("共變異係數矩陣=",cov_mat)

#求共變異係數矩陣 的特徵向量及特徵值
eigen_vals, eigen_vecs = np.linalg.eig(cov_mat)
print("特徵向量.shape=",eigen_vecs.shape)
print("特徵向量=",eigen_vecs)
print("特徵值=",eigen_vals)


#計算解釋變異數比率 各特徵值/特徵值總和
tot = sum(eigen_vals)
var_exp = [(i / tot) for i in sorted(eigen_vals, reverse=True)]
cum_var_exp = np.cumsum(var_exp)
print("各特徵值變異數比率:",var_exp)
print("特徵值變異數比率累加:",cum_var_exp)

#畫圖 :解釋變異數比率 ,各特徵值/特徵值總和
plt.bar(range(1, 14), var_exp, alpha=0.5, align='center',
        label='individual explained variance')
plt.step(range(1, 14), cum_var_exp, where='mid',
         label='cumulative explained variance')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal components')
plt.legend(loc='best')
plt.tight_layout()
plt.show()

# Make a list of (eigenvalue, eigenvector) tuples
eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:, i])
               for i in range(len(eigen_vals))]
print("特徵值,特徵向量length:",len(eigen_pairs))
print("特徵值,特徵向量:",eigen_pairs)
# Sort the (eigenvalue, eigenvector) tuples from high to low
eigen_pairs.sort(key=lambda k: k[0], reverse=True)
print("特徵值,特徵向量排序:",eigen_pairs.sort(key=lambda k: k[0], reverse=True))
#保留兩個最具影響力的特徵向量組成13x2 的投影矩陣W
w = np.hstack((eigen_pairs[0][1][:, np.newaxis],
               eigen_pairs[1][1][:, np.newaxis]))
print('Matrix W:\n', w)

#畫出轉換後的數據集 散點圖
#print("X_train_std[0].dot(w)=",X_train_std[0].dot(w))
X_train_pca = X_train_std.dot(w)
colors = ['r', 'b', 'g']
markers = ['s', 'x', 'o']
for l, c, m in zip(np.unique(y_train), colors, markers):
    plt.scatter(X_train_pca[y_train == l, 0], 
                X_train_pca[y_train == l, 1], 
                c=c, label=l, marker=m)

plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.tight_layout()
plt.show()

##劃出決策分布圖
from matplotlib.colors import ListedColormap
def plot_decision_regions(X, y, classifier, resolution=0.02):
    # setup marker generator and color map
    markers = ('s', 'x', 'o', '^', 'v')
    colors = ('red', 'blue', 'lightgreen', 'gray', 'cyan')
    cmap = ListedColormap(colors[:len(np.unique(y))])
    # plot the decision surface
    x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution),
                           np.arange(x2_min, x2_max, resolution))
    Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
    Z = Z.reshape(xx1.shape)
    plt.contourf(xx1, xx2, Z, alpha=0.4, cmap=cmap)
    plt.xlim(xx1.min(), xx1.max())
    plt.ylim(xx2.min(), xx2.max())
    # plot class samples
    for idx, cl in enumerate(np.unique(y)):
        plt.scatter(x=X[y == cl, 0], y=X[y == cl, 1],
                    alpha=0.8, c=cmap(idx),
                    marker=markers[idx], label=cl)
     
#使用LogisticRegression 並用兩個PCA 主成分做訓練分類 
from sklearn.linear_model import LogisticRegression
lr = LogisticRegression()
lr = lr.fit(X_train_pca, y_train)
plot_decision_regions(X_train_pca, y_train, classifier=lr)
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower right')
plt.tight_layout()
plt.show()






<參考資料>書名:Python機器學習,作者:Sebastian Raschka

加入阿布拉機的3D列印與機器人的FB專頁
https://www.facebook.com/arbu00/


<其他相關文章>
人工神經網路(1)--使用Python實作perceptron(感知器)
人工神經網路(2)--使用Python實作後向傳遞神經網路演算法(Backprogation artificial neature network)
深度學習(1)-如何在windows安裝Theano +Keras +Tensorflow並使用GPU加速訓練神經網路
機器學習(1)--使用OPENCV KNN實作手寫辨識
機器學習(2)--使用OPENCV SVM實作手寫辨識
演算法(1)--蒙地卡羅法求圓周率及橢圓面積(Monte carlo)
機器學習(3)--適應線性神經元與梯度下降法(Adaline neuron and Gradient descent)
機器學習(4)--資料標準常態化與隨機梯度下降法( standardization & Stochastic Gradient descent)
機器學習(5)--邏輯斯迴歸,過度適合與正規化( Logistic regression,overfitting and regularization)