主成分分析（Principal Component Analysis，PCA）， 是一种统计方法。通过正交变换将一组可能存在相关性的变量转换为一组线性不相关的变量，转换后的这组变量叫主成分。

在实际课题中，为了全面分析问题，往往提出很多与此有关的变量（或因素），因为每个变量都在不同程度上反映这个课题的某些信息。

主成分分析首先是由K.皮尔森（Karl Pearson）对非随机变量引入的，尔后H.霍特林将此方法推广到随机向量的情形。信息的大小通常用离差平方和或方差来衡量。

# -*- coding: utf-8 -*-
import numpy as np
from matplotlib import pyplot as plt
from scipy import io as spio
from sklearn.decomposition import pca

'''
主成分分析_2维数据降维1维演示函数
'''


def PCA_2D():
    data_2d = spio.loadmat("data.mat")
    X = data_2d['X']
    m = X.shape[0]
    plt = plot_data_2d(X, 'bo')  # 显示二维的数据
    plt.show()

    X_copy = X.copy()
    X_norm, mu, sigma = featureNormalize(X_copy)  # 归一化数据
    # plot_data_2d(X_norm)    # 显示归一化后的数据
    # plt.show()

    Sigma = np.dot(np.transpose(X_norm), X_norm) / m  # 求Sigma
    U, S, V = np.linalg.svd(Sigma)  # 求Sigma的奇异值分解

    plt = plot_data_2d(X, 'bo')  # 显示原本数据
    drawline(plt, mu, mu + S[0] * (U[:, 0]), 'r-')  # 线，为投影的方向

    plt.axis('square')
    plt.show()

    K = 1  # 定义降维多少维（本来是2维的，这里降维1维）
    '''投影之后数据（降维之后）'''
    Z = projectData(X_norm, U, K)  # 投影
    '''恢复数据'''
    X_rec = recoverData(Z, U, K)  # 恢复
    '''作图-----原数据与恢复的数据'''
    plt = plot_data_2d(X_norm, 'bo')
    plot_data_2d(X_rec, 'ro')
    for i in range(X_norm.shape[0]):
        drawline(plt, X_norm[i, :], X_rec[i, :], '--k')
    plt.axis('square')
    plt.show()


'''主成分分析_PCA图像数据降维'''


def PCA_faceImage():
    print(u'加载图像数据.....')
    data_image = spio.loadmat('data_faces.mat')
    X = data_image['X']
    display_imageData(X[0:100, :])
    m = X.shape[0]  # 数据条数

    print(u'运行PCA....')
    X_norm, mu, sigma = featureNormalize(X)  # 归一化

    Sigma = np.dot(np.transpose(X_norm), X_norm) / m  # 求Sigma
    U, S, V = np.linalg.svd(Sigma)  # 奇异值分解
    display_imageData(np.transpose(U[:, 0:36]))  # 显示U的数据

    print(u'对face数据降维.....')
    K = 100  # 降维100维(原先是32*32=1024维的)
    Z = projectData(X_norm, U, K)
    print(u'投影之后Z向量的大小：%d %d' % Z.shape)

    print(u'显示降维之后的数据......')
    X_rec = recoverData(Z, U, K)  # 恢复数据
    display_imageData(X_rec[0:100, :])


# 可视化二维数据
def plot_data_2d(X, marker):
    plt.plot(X[:, 0], X[:, 1], marker)
    return plt


# 归一化数据
def featureNormalize(X):
    '''（每一个数据-当前列的均值）/当前列的标准差'''
    n = X.shape[1]
    mu = np.zeros((1, n));
    sigma = np.zeros((1, n))

    mu = np.mean(X, axis=0)  # axis=0表示列
    sigma = np.std(X, axis=0)
    for i in range(n):
        X[:, i] = (X[:, i] - mu[i]) / sigma[i]
    return X, mu, sigma


# 映射数据
def projectData(X_norm, U, K):
    Z = np.zeros((X_norm.shape[0], K))

    U_reduce = U[:, 0:K]  # 取前K个
    Z = np.dot(X_norm, U_reduce)
    return Z


# 画一条线
def drawline(plt, p1, p2, line_type):
    plt.plot(np.array([p1[0], p2[0]]), np.array([p1[1], p2[1]]), line_type)


# 恢复数据
def recoverData(Z, U, K):
    X_rec = np.zeros((Z.shape[0], U.shape[0]))
    U_recude = U[:, 0:K]
    X_rec = np.dot(Z, np.transpose(U_recude))  # 还原数据（近似）
    return X_rec


# 显示图片
def display_imageData(imgData):
    sum = 0
    '''
    显示100个数（若是一个一个绘制将会非常慢，可以将要画的图片整理好，放到一个矩阵中，显示这个矩阵即可）
    - 初始化一个二维数组
    - 将每行的数据调整成图像的矩阵，放进二维数组
    - 显示即可
    '''
    m, n = imgData.shape
    width = np.int32(np.round(np.sqrt(n)))
    height = np.int32(n / width);
    rows_count = np.int32(np.floor(np.sqrt(m)))
    cols_count = np.int32(np.ceil(m / rows_count))
    pad = 1
    display_array = -np.ones((pad + rows_count * (height + pad), pad + cols_count * (width + pad)))
    for i in range(rows_count):
        for j in range(cols_count):
            max_val = np.max(np.abs(imgData[sum, :]))
            display_array[pad + i * (height + pad):pad + i * (height + pad) + height,
            pad + j * (width + pad):pad + j * (width + pad) + width] = imgData[sum, :].reshape(height, width,
                                                                                               order="F") / max_val  # order=F指定以列优先，在matlab中是这样的，python中需要指定，默认以行
            sum += 1

    plt.imshow(display_array, cmap='gray')  # 显示灰度图像
    plt.axis('off')
    plt.show()


if __name__ == "__main__":
    PCA_2D()
    PCA_faceImage()


作者：WangB

在上一篇文章中，我汇总了大佬们写的关于代谢组学的一些基础知识，让大家对代谢组学有了一个整体性的认识。这次给大家分享的是一些如何对数据进行可视化的文章。
Principal Component Analysis （主成分分析）
很多同学在网上搜的时候，会发现大多数关于主成分分析的文章都是给你从数学理论方面推导来说明什么是PCA，而像我这种数学小白，表示数学公式都看不懂，直到看到一个不存在的网站上关于PCA的视频。自己一个人看肯定不够爽，然后就征得原作者同意，把系列视频搬运到了B站，关于这个，大家可以看一下这篇文章：
久久琼殷：StatQuest with Josh Starmer 来到中国啦​zhuanlan.zhihu.com
以下是关于PCA系列的视频，一定要多看几遍，不懂的时候再回来多看几遍：
StatQuest - PCA中的主要概念​www.bilibili.comStatQuest -主成分分析（Principal Component Analysis)​www.bilibili.comStatQuest - 主成分分析（更新版）​www.bilibili.comStatQuest - R实现主成分分析 （PCA in R）​www.bilibili.com
如果想利用 R 来绘制PCA的话，更多的细节可以仔细阅读下面的文章（也可百度一下）：
PCA - Principal Component Analysis Essentials​www.sthda.com
一文看懂主成分分析​mp.weixin.qq.com
如果你有一些数据分析软件的话，可以参考以下的文章（需要对应软件）：
【干货】SIMCA14.1操作教程-杨树胁迫案例（GC-MS）​mp.weixin.qq.com
JMP——Example of PCA​www.jmp.com
医小咖：SPSS超详细教程：主成分分析​zhuanlan.zhihu.com
当然，你可以用在线的网站进行分析：
MetaboAnalyst​www.metaboanalyst.ca
Multivariate data analysis​scienceinside.shinyapps.io

主成分分析
主成分分析（Principal Component Analysis，PCA，通过正交变换将一组可能存在相关性的变量转换为一组线性不相关的变量，转换后的这组变量叫主成分。设法将原来变量重新组合成一组新的互相无关的几个综合变量，同时根据实际需要从中可以取出几个较少的综合变量尽可能多地反映原来变量的信息的统计方法叫做主成分分析或称主分量分析，也是数学上用来降维的一种方法。
1.准备数据
在这里我们使用sklearn提供的鸢尾花数据集。
import numpy as np
from sklearn import datasets

 # 加载鸢尾花数据集中的特征作为原始数据集
    def loadIris(self):
        data = datasets.load_iris()["data"]
        return data

    # 标准化数据
    def Standard(self,data):
        # axis=0按列取均值
        mean_vector=np.mean(data,axis=0)
        return mean_vector,data - mean_vector

2.计算得到协方差矩阵
    def getCovMatrix(self,newData):
        # rowvar=0表示数据的每一列代表一个feature
        return np.cov(newData,rowvar=0)

    # 计算协方差矩阵的特征值和特征向量
    def getFValueAndFVector(self,covMatrix):
        fValue,fVector = np.linalg.eig(covMatrix)
        return fValue,fVector

    # 得到特征向量矩阵
    def getVectorMatrix(self,fValue,fVector,k):
        fValueSort = np.argsort(fValue)
        fValueTopN = fValueSort[:-(k + 1):-1]
        return fVector[:,fValueTopN]

    # 得到降维后的数据
    def getResult(self,data,vectorMatrix):
        return np.dot(data,vectorMatrix)

3.计算特征值和特征向量
	def getFValueAndFVector(self,covMatrix):
        fValue,fVector = np.linalg.eig(covMatrix)
        return fValue,fVector
 	# 得到特征向量矩阵
   	def getVectorMatrix(self,fValue,fVector,k):
        fValueSort = np.argsort(fValue)
        fValueTopN = fValueSort[:-(k + 1):-1]
        return fVector[:,fValueTopN]

4.得到降维后的数据
	def getResult(self,data,vectorMatrix):
        return np.dot(data,vectorMatrix)

主成分分析是一种降维方法，主要用于数据压缩，数据可视化以及特征提取等方面。 现实中我们经常可以遇到维数很高的数据，如一张28*28的图片，可以看作维度为784。类似图片这样的高维数据，实际上各个维度之间具有高度的关联性，即维度之间并非完全独立的。通过进行主成分分析,可以将数据的主要特征提取出来,忽略非 重要特征,从而对数据进行压缩。这里的主成分代表一组规范正交基,每个基用表示,并且满足。假设X表示我们的观测数据矩阵,大小为N*D,
 我们希望将其压缩到 N*M的矩阵Z。主成分分析有两种直观上的理解,并且可以 证明这两种直观理解是等价的。

第一种方式,我们希望找到这样一组正交基,使得映射后的数据方差最大。这样的目标函数符合我们的直观认识,因为方差表示了数据的差异性，方差大表示数据的差异性大, 即信息量大。如果在某个方向上的方差太小，或者等于 0，那么该方向的信息量小，因此在要求不高的情况下，可以忽略这样的方向。主成分分析正是基于这样的直观认识进行处理的。为了书写简单，我们假设预测数据的均值是0，如果不是0，让每个数据减去均值。一个数据 映射到方向u后的坐标可以记为。
 因此映射后在该方向上的方差为


    (1)


其中

              （2）

为了使映射方差最大,使用拉格朗日乘子法最大化以下目标函数





            （3）





对u进行求导,然后令其等于零,可以得到如下解,

                   (4)

因此,u是数据协方差矩阵S的特征向量,进一步得到

                 (5)


即λ是一个极值点。因此u对应的是最大特征值的 特征向量,并且方差的最大值是S的最大特征值。换句话

说, S的各个特征值是数据在各个方向上映射方差的极值点。 另外计算一个矩阵A在某个方向u上的方差可

以如下计算

  
   (6)

即它表示u在各个正交基上的方差的加权和,权重就是相应方向的 特征值，如果u恰好是一个特征向量,那

么该结果就变成了。 这说明了特征值和方差之间的内在关系。




主成分分析的另一种解释是找到一组正交基使得映射后的误差 最小。同样我们用表示我们的正次基，并且假设我们的主成分空间大小为M，数据空间大小为D，那么对于每个点我们可以有如下近似，

                (7)


上面的近似理解为，对 于前M个向量，我们为每个数据每给出精确的映射坐标，而对于非主成分，所有数据的值都用一个常数表示。于是我们得到如下的损失函数，

             (8)


对损失函数求导后，得到以下解，


                   (9)





         (10)




      
   (11)





然后将式(11)代入式(8)得到如下




  
   (12)





从上面的形式可以看出,它与最大化方差的形式相同,即的形式, 只不过现在是最小化。通过简单求导即可导出,当u是特征向量,可以取得 极值λ, 因此为了使得上式最小可以选择特征值最小的特征向量作为非主成分,相应的特征值最大的特征向量作为主成分。

上面给出了主成分分析的两种直观上的理解，一种是最大化方差，一种是最小化损失函数，无论用哪种方法，得到的结果是相同的。通过主成分分析我们发现，一个矩阵的特征值实际上是该矩阵在各个方向上方差的极值，并且在特征向量处取得。




1. pattern recognition and machine learning              Christopher M.Bishop   p559-565