主成分分析:Principal Component Analysis.PCA

PCA就是数据降维，最重要的能力就是数据可视化（高维数据）。

一个二维数据，我们可以在这个图中找到一根直线，使得所有的数据点在这根直线上的投影差异化最大（其实就是方差最大），很显然，在A，B，C三条线中，B是最符合的

PCA的本质就是找一些相互正交的投影方向，使得数据在这些投影方向上的方差最大（找新的一组正交基）。计算原始数据在这些正交基上投影的方差，方差越大，就说明在对应正交基上包含了更多的信息量。
原始数据协方差矩阵的特征值越大，对应的方差越大，在对应的特征向量上投影的信息量就越大，就是主成分。反之特征值小，数据在这些特征向量上投影的信息量很小，可以将小特征值对应方向的数据删除，从而达到了降维的目的。
我们可以选择最大的n个特征值，将数据降到n维，做个归一化然后看这几个特征值的权重，可以是超过80%或者95%，具体大小看业务的要求。

PCA的主要逻辑

* 去除平均值
* 计算协方差矩阵
* 计算协方差矩阵的特征值和特征向量
* 将特征值从大到小排序
* 保留最大的N个特征向量
* 将数据转换到上述N个特征向量构建的新空间中

教材实例（一）：某工厂半导体特征数据

《ML in Action》p246，

数据源

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

#格式化数据函数
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) #协方差矩阵
    # rowvar=0很重要，表明每一行是一个样本，否则numpy认为整个数据是一个样本
    eigVals,eigVects = linalg.eig(mat(covMat))
    # numpy内建函数，一次性求出特征值和特征向量
    eigValInd = argsort(eigVals)
    eigValInd = eigValInd[:-(topNfeat+1):-1]
    redEigVects = eigVects[:,eigValInd]
    # 对特征值从小到大排序
    lowDDataMat = meanRemoved * redEigVects
    # 将数据集投影到新的空间，结果是一个低维数据
    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打印出来以后（由于结果太冗长就不再打印出来），会发现前30个特征值就占了99%以上的权重，如果这符合业务要求的话，那么我们就成功地将一个590维数据降到了30维。

也可以写成权重百分比版本：

def pca(dataMat,percentage=0.99):
    newData,meanVal=zeroMean(dataMat)
    covMat=np.cov(newData,rowvar=0)    #求协方差矩阵,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个特征值对应的特征向量
    lowDDataMat=newData*n_eigVect               #低维特征空间的数据
    reconMat=(lowDDataMat*n_eigVect.T)+meanVal  #重构数据
    return lowDDataMat,reconMat

教材实例（二）：鸢尾花数据 Iris

PCA经典数据，数据源
利用sklearn中现有的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说明

参数：

• n_components: 主成分个数。
  赋值int，比如n_components=1，将把原始数据降到一个维度。
  赋值string，比如n_components='mle'，将自动选取特征个数n，使得满足所要求的方差百分比。
  默认为None，只求特征值，不做降维处理。
• copy:bool，True或者False，缺省时默认为True。表示是否在运行算法时，将原始训练数据复制一份。若为True，在原始数据的副本上进行运算；若为False，则运行PCA算法后原始训练数据的值会改。
• whiten: bool，缺省时默认为False白化，使得每个特征具有相同的方差。白化

属性：

• components_ ：返回具有最大方差的成分。
• explained_variance_ratio_：返回 所保留的n个成分各自的方差百分比。
• n_components_：返回所保留的成分个数n。
• mean_：
• noise_variance_：

方法：

• fit(X,y=None)
  fit()可以说是scikit-learn中通用的方法，每个需要训练的算法都会有fit()方法，它其实就是算法中的“训练”这一步骤。因为PCA是无监督学习算法，此处y自然等于None。fit(X)，表示用数据X来训练PCA模型。
  函数返回值：调用fit方法的对象本身。比如pca.fit(X)，表示用X对pca这个对象进行训练。
• fit_transform(X)
  用X来训练PCA模型，同时返回降维后的数据。
  newX=pca.fit_transform(X)，newX就是降维后的数据。
• inverse_transform()
  将降维后的数据转换成原始数据，X=pca.inverse_transform(newX)
• transform(X)
  将数据X转换成降维后的数据。当模型训练好后，对于新输入的数据，都可以用transform方法来降维。
• 此外，还有get_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'鸢尾花数据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转换成ndarray
        x1_label, x2_label = columns[idx]
        title = u'鸢尾花数据特征选择'

    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 '最优参数：', 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)                    # 生成网格采样点
    x_show = np.stack((x1.flat, x2.flat), axis=1)   # 测试点
    y_hat = model.predict(x_show)  # 预测值
    y_hat = y_hat.reshape(x1.shape)  # 使之与输入的形状相同
    plt.figure(facecolor='w')
    plt.pcolormesh(x1, x2, y_hat, cmap=cm_light)  # 预测值的显示
    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同独立成分分析等方法一样，是不假设样本处于何种具体分布的。