主成分分析:Principal Component Analysis.PCA

PCA就是數(shù)據(jù)降維，最重要的能力就是數(shù)據(jù)可視化（高維數(shù)據(jù)）。

一個二維數(shù)據(jù)，我們可以在這個圖中找到一根直線，使得所有的數(shù)據(jù)點在這根直線上的投影差異化最大（其實就是方差最大），很顯然，在A，B，C三條線中，B是最符合的

PCA的本質(zhì)就是找一些相互正交的投影方向，使得數(shù)據(jù)在這些投影方向上的方差最大（找新的一組正交基）。計算原始數(shù)據(jù)在這些正交基上投影的方差，方差越大，就說明在對應(yīng)正交基上包含了更多的信息量。
原始數(shù)據(jù)協(xié)方差矩陣的特征值越大，對應(yīng)的方差越大，在對應(yīng)的特征向量上投影的信息量就越大，就是主成分。反之特征值小，數(shù)據(jù)在這些特征向量上投影的信息量很小，可以將小特征值對應(yīng)方向的數(shù)據(jù)刪除，從而達到了降維的目的。
我們可以選擇最大的n個特征值，將數(shù)據(jù)降到n維，做個歸一化然后看這幾個特征值的權(quán)重，可以是超過80%或者95%，具體大小看業(yè)務(wù)的要求。

PCA的主要邏輯

* 去除平均值
* 計算協(xié)方差矩陣
* 計算協(xié)方差矩陣的特征值和特征向量
* 將特征值從大到小排序
* 保留最大的N個特征向量
* 將數(shù)據(jù)轉(zhuǎn)換到上述N個特征向量構(gòu)建的新空間中

教材實例（一）：某工廠半導體特征數(shù)據(jù)

《ML in Action》p246，

數(shù)據(jù)源

# -*- coding:utf-8 -*-
from numpy import *

#格式化數(shù)據(jù)函數(shù)
def loadDataSet(filename, delim='\t'):
    fr = open(filename)
    stringArr = [line.strip().split(delim) for line in fr.readlines()]
    datArr = [map(float, line) for line in stringArr]
    return mat(datArr)

def pca(dataMat, topNfeat=999999):
    meanVals = mean(dataMat)
    # meanVals = mean(dataMat, axis=0)
    # axis =0 說明按列求均值，否則是按行
    meanRemoved = dataMat - meanVals #去均值
    covMat = cov(meanRemoved, rowvar=0) #協(xié)方差矩陣
    # rowvar=0很重要，表明每一行是一個樣本，否則numpy認為整個數(shù)據(jù)是一個樣本
    eigVals,eigVects = linalg.eig(mat(covMat))
    # numpy內(nèi)建函數(shù)，一次性求出特征值和特征向量
    eigValInd = argsort(eigVals)
    eigValInd = eigValInd[:-(topNfeat+1):-1]
    redEigVects = eigVects[:,eigValInd]
    # 對特征值從小到大排序
    lowDDataMat = meanRemoved * redEigVects
    # 將數(shù)據(jù)集投影到新的空間，結(jié)果是一個低維數(shù)據(jù)
    reconMat = (lowDDataMat * redEigVects.T) + meanVals
    return lowDDataMat, reconMat

# 缺失值處理
def replaceNanWithMean():
    datMat = loadDataSet('secom.data', ' ')
    numFeat = shape(datMat)[1]
    for i in range(numFeat):
        meanVal = mean(datMat[nonzero(~isnan(datMat[:,i].A))[0],i])
        datMat[nonzero(isnan(datMat[:,i].A))[0],i] = meanVal
    return datMat

if __name__ == '__main__':
    dataMat = replaceNanWithMean()
    meanVals = mean(dataMat, axis=0)
    meanRemoved = dataMat - meanVals
    covMat = cov(meanRemoved, rowvar=0)
    eigVals, eigVects = linalg.eig(mat(covMat))

當我們試著把eigVals打印出來以后（由于結(jié)果太冗長就不再打印出來），會發(fā)現(xiàn)前30個特征值就占了99%以上的權(quán)重，如果這符合業(yè)務(wù)要求的話，那么我們就成功地將一個590維數(shù)據(jù)降到了30維。

也可以寫成權(quán)重百分比版本：

def pca(dataMat,percentage=0.99):
    newData,meanVal=zeroMean(dataMat)
    covMat=np.cov(newData,rowvar=0)    #求協(xié)方差矩陣,return ndarray；若rowvar非0，一列代表一個樣本，為0，一行代表一個樣本
    eigVals,eigVects=np.linalg.eig(np.mat(covMat))#求特征值和特征向量,特征向量是按列放的，即一列代表一個特征向量
    n=percentage2n(eigVals,percentage)                 #要達到percent的方差百分比，需要前n個特征向量
    eigValIndice=np.argsort(eigVals)            #對特征值從小到大排序
    n_eigValIndice=eigValIndice[-1:-(n+1):-1]   #最大的n個特征值的下標
    n_eigVect=eigVects[:,n_eigValIndice]        #最大的n個特征值對應(yīng)的特征向量
    lowDDataMat=newData*n_eigVect               #低維特征空間的數(shù)據(jù)
    reconMat=(lowDDataMat*n_eigVect.T)+meanVal  #重構(gòu)數(shù)據(jù)
    return lowDDataMat,reconMat

教材實例（二）：鳶尾花數(shù)據(jù) Iris

PCA經(jīng)典數(shù)據(jù)，數(shù)據(jù)源
利用sklearn中現(xiàn)有的PCA工具。

Data Preview:
5.1,3.5,1.4,0.2,Iris-setosa
4.9,3.0,1.4,0.2,Iris-setosa
5.5,2.3,4.0,1.3,Iris-versicolor
...
Attribute Information:
   1. sepal length in cm
   2. sepal width in cm
   3. petal length in cm
   4. petal width in cm
   5. class: 
      -- Iris Setosa
      -- Iris Versicolour
      -- Iris Virginica

代碼如下：

# -*- coding:utf-8 -*-

import pandas as pd
import numpy as np
from sklearn.decomposition import PCA
from sklearn.feature_selection import SelectKBest, SelectPercentile, chi2
from sklearn.linear_model import LogisticRegressionCV
from sklearn import metrics
from sklearn.model_selection import train_test_split
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures


def extend(a, b):
    return 1.05*a-0.05*b, 1.05*b-0.05*a

if __name__ == '__main__':
    pca = True
    pd.set_option('display.width', 200)
    data = pd.read_csv('iris.data', header=None)
    # columns = np.array(['sepal_length', 'sepal_width', 'petal_length', 'petal_width', 'type'])
    columns = np.array([u'花萼長度', u'花萼寬度', u'花瓣長度', u'花瓣寬度', u'類型'])
    data.rename(columns=dict(zip(np.arange(5), columns)), inplace=True)

    data[u'類型'] = pd.Categorical(data[u'類型']).codes 
    #pd的類型變量功能
    print data.head(5)

    x = data[columns[:-1]]
    y = data[columns[-1]]

sklearn.PCA說明

參數(shù)：

• n_components: 主成分個數(shù)。
  賦值int，比如n_components=1，將把原始數(shù)據(jù)降到一個維度。
  賦值string，比如n_components='mle'，將自動選取特征個數(shù)n，使得滿足所要求的方差百分比。
  默認為None，只求特征值，不做降維處理。
• copy:bool，True或者False，缺省時默認為True。表示是否在運行算法時，將原始訓練數(shù)據(jù)復(fù)制一份。若為True，在原始數(shù)據(jù)的副本上進行運算；若為False，則運行PCA算法后原始訓練數(shù)據(jù)的值會改。
• whiten: bool，缺省時默認為False白化，使得每個特征具有相同的方差。白化

屬性：

• components_ ：返回具有最大方差的成分。
• explained_variance_ratio_：返回 所保留的n個成分各自的方差百分比。
• n_components_：返回所保留的成分個數(shù)n。
• mean_：
• noise_variance_：

方法：

• fit(X,y=None)
  fit()可以說是scikit-learn中通用的方法，每個需要訓練的算法都會有fit()方法，它其實就是算法中的“訓練”這一步驟。因為PCA是無監(jiān)督學習算法，此處y自然等于None。fit(X)，表示用數(shù)據(jù)X來訓練PCA模型。
  函數(shù)返回值：調(diào)用fit方法的對象本身。比如pca.fit(X)，表示用X對pca這個對象進行訓練。
• fit_transform(X)
  用X來訓練PCA模型，同時返回降維后的數(shù)據(jù)。
  newX=pca.fit_transform(X)，newX就是降維后的數(shù)據(jù)。
• inverse_transform()
  將降維后的數(shù)據(jù)轉(zhuǎn)換成原始數(shù)據(jù)，X=pca.inverse_transform(newX)
• transform(X)
  將數(shù)據(jù)X轉(zhuǎn)換成降維后的數(shù)據(jù)。當模型訓練好后，對于新輸入的數(shù)據(jù)，都可以用transform方法來降維。
• 此外，還有g(shù)et_covariance()、get_precision()、get_params(deep=True)、score(X, y=None)等方法。

    if pca:
        pca = PCA(n_components=2, whiten=True, random_state=0)
        x = pca.fit_transform(x)
        print '各方向方差：', pca.explained_variance_
        print '方差所占比例：', pca.explained_variance_ratio_
        x1_label, x2_label = u'組分1', u'組分2'
        title = u'鳶尾花數(shù)據(jù)PCA降維'
    else:
        fs = SelectKBest(chi2, k=2)
        # fs = SelectPercentile(chi2, percentile=60)
        fs.fit(x, y)
        idx = fs.get_support(indices=True)
        print 'fs.get_support() = ', idx
        x = x[idx]
        x = x.values    # 為下面使用方便，DataFrame轉(zhuǎn)換成ndarray
        x1_label, x2_label = columns[idx]
        title = u'鳶尾花數(shù)據(jù)特征選擇'

    print x[:5]
    cm_light = mpl.colors.ListedColormap(['#77E0A0', '#FF8080', '#A0A0FF'])
    cm_dark = mpl.colors.ListedColormap(['g', 'r', 'b'])
    mpl.rcParams['font.sans-serif'] = u'SimHei'
    mpl.rcParams['axes.unicode_minus'] = False
    plt.figure(facecolor='w')
    plt.scatter(x[:, 0], x[:, 1], s=30, c=y, marker='o', cmap=cm_dark)
    plt.grid(b=True, ls=':')
    plt.xlabel(x1_label, fontsize=14)
    plt.ylabel(x2_label, fontsize=14)
    plt.title(title, fontsize=18)
    # plt.savefig('1.png')
    plt.show()

    x, x_test, y, y_test = train_test_split(x, y, train_size=0.7)
    model = Pipeline([
        ('poly', PolynomialFeatures(degree=2, include_bias=True)),
        ('lr', LogisticRegressionCV(Cs=np.logspace(-3, 4, 8), cv=5, fit_intercept=False))
    ])
    model.fit(x, y)
    print '最優(yōu)參數(shù)：', model.get_params('lr')['lr'].C_
    y_hat = model.predict(x)
    print '訓練集精確度：', metrics.accuracy_score(y, y_hat)
    y_test_hat = model.predict(x_test)
    print '測試集精確度：', metrics.accuracy_score(y_test, y_test_hat)

    N, M = 500, 500     # 橫縱各采樣多少個值
    x1_min, x1_max = extend(x[:, 0].min(), x[:, 0].max())   # 第0列的范圍
    x2_min, x2_max = extend(x[:, 1].min(), x[:, 1].max())   # 第1列的范圍
    t1 = np.linspace(x1_min, x1_max, N)
    t2 = np.linspace(x2_min, x2_max, M)
    x1, x2 = np.meshgrid(t1, t2)                    # 生成網(wǎng)格采樣點
    x_show = np.stack((x1.flat, x2.flat), axis=1)   # 測試點
    y_hat = model.predict(x_show)  # 預(yù)測值
    y_hat = y_hat.reshape(x1.shape)  # 使之與輸入的形狀相同
    plt.figure(facecolor='w')
    plt.pcolormesh(x1, x2, y_hat, cmap=cm_light)  # 預(yù)測值的顯示
    plt.scatter(x[:, 0], x[:, 1], s=30, c=y, edgecolors='k', cmap=cm_dark)  # 樣本的顯示
    plt.xlabel(x1_label, fontsize=14)
    plt.ylabel(x2_label, fontsize=14)
    plt.xlim(x1_min, x1_max)
    plt.ylim(x2_min, x2_max)
    plt.grid(b=True, ls=':')
    patchs = [mpatches.Patch(color='#77E0A0', label='Iris-setosa'),
              mpatches.Patch(color='#FF8080', label='Iris-versicolor'),
              mpatches.Patch(color='#A0A0FF', label='Iris-virginica')]
    plt.legend(handles=patchs, fancybox=True, framealpha=0.8, loc='lower right')
    plt.title(u'鳶尾花Logistic回歸分類效果', fontsize=17)
    # plt.savefig('2.png')
    plt.show()

最后要注意的是，PCA同獨立成分分析等方法一樣，是不假設(shè)樣本處于何種具體分布的。