python和spss的主成分分析结果不一样 python主成分分析实例

文章目录

• 前置知识
• PCA

• 均值和协方差

• 均值
• 协方差矩阵

• 主成分分析

• 变量代换
• 多维数据的降维
• 算法步骤

• 代码实现

• 特征值分解
• 奇异值分解

前置知识

• 二次型
• 对角化
• 奇异值分解
• 施密特正交化

PCA

主成分分析是一种从原始数据中消除冗余信息，只需要一个或者几个属性合成的属性就可以提供大部分信息的方法。大概过程是找出一个特殊的线性组合，给要分析的元素赋予一个权值，综合得到一个新的数据。

假设每一个样本有p个属性，一共有n个样本，那么这些向量可以构成一个$p\times{n}$的矩阵，称为观测矩阵。

主成分分析实际上相当于变换坐标轴，并降低维度。

如一个二维数据集如下图：

经过PCA进行降维之后可以使用一个维度进行表示，如下图

均值和协方差

均值

令$[\boldsymbol{X}_1\cdots \boldsymbol{X}_N ]$为一个$p\times{N}$观测矩阵，则其样本平均值M为：

$\boldsymbol{M}=\frac{1}{N}(\boldsymbol{X}_1+\cdots+\boldsymbol{X}_N)$

令

$\boldsymbol{\hat{X}}_k=\boldsymbol{X}_k-\boldsymbol{M},k=1,2,\cdots,N$

则新矩阵

$B=[\boldsymbol{\hat{X}}_1\cdots\boldsymbol{\hat{X}}_N]$

具有零样本均值，这种矩阵

$B$被称为

平均偏差形式

协方差矩阵

协方差矩阵是一个$p\times{p}$的矩阵$S$,其定义为

$S=\frac{1}{N-1}BB^T$

由于任何具有$BB^T$形式的矩阵都是半正定的，所以$S$也是半正定的。

PS： 虽然协方差矩阵是半正定矩阵，但是由于计算误差，算出来的特征值可能是负数。

其中的对角线元素$s_{jj}$称为$x_j$的方差。方差用来度量$x_j$的分散性。

数据的总方差被称为矩阵的迹，即$S$的对角线元素之和，记做$tr(S)$

$S$中的元素$s_{ij},i\not ={j}$称为$x_i$和$x_j$的协方差。若协方差为0，则称$x_i$和$x_j$是无关的。当协方差矩阵是对角阵或几乎是对角阵时，$\boldsymbol{X}_1 \cdots\boldsymbol{X}_N$中多变量分析可以简化。

主成分分析

变量代换

设矩阵$[\boldsymbol{X}_1\cdots \boldsymbol{X}_N ]$已经称为平均偏差形式，主成分分析的目标是找到一个$p\times{p}$的正交矩阵$P=[\boldsymbol{u}_1 \cdots \boldsymbol{u}_p]$,确定一个变量代换$\boldsymbol{X}=P\boldsymbol{Y}$,或
$\left[\begin{matrix} x_1 \\ x_2 \\ \vdots \\ x_p \end{matrix}\right] =[\boldsymbol{u_1} \cdots \boldsymbol{u_p}]\left[\begin{matrix} y_1 \\ y_2 \\ \vdots \\ y_p \end{matrix}\right]$
满足新的变量$y_1,\cdots,y_p$两两无关，并保证整理后的方差具有递减顺序。
变量正交变换$X=PY$说明，每一个观测向量$\boldsymbol{X}_k$得到一个新”名称“$\boldsymbol{Y}_k$,即$\boldsymbol{Y}_k$是$\boldsymbol{X}_k$以$P$为基的坐标。
易得对任何正交矩阵$P$，$\boldsymbol{Y}_1,\cdots,\boldsymbol{Y_N}$的协方差是$P^TSP$，于是可以构造符合要求的P。设$D$是对角矩阵且$S$的特征值在$D$的对角线上，重新整理使得$\lambda_1\geq\lambda_2\geq\cdots\geq\lambda_p\geq 0$即
$D = \left[\begin{matrix} \lambda_1 & 0 & \cdots & 0 \\ 0 & \lambda_2 & \cdots & 0 \\ \vdots & \enspace & \ddots & \vdots \\ 0 & \cdots &\cdots & \lambda_p \end{matrix}\right]$找到正交矩阵$P$，它的列对应单位特征向量$\boldsymbol{u}_1,\cdots\boldsymbol{u}_p$,且有$S=PDP^T$。
协方差矩阵$S$的单位特征向量$\boldsymbol{u}_1,\cdots\boldsymbol{u}_p$称为(观测矩阵中)数据的主成分。第一主成分是$S$中最大特征值对应的特征向量，第二主成分以此类推。
设$c1,c2,\cdots,c_p$为$\boldsymbol{u}_1$中的元素，第一主成分$\boldsymbol{u}_1$可以通过公式$y_1=\boldsymbol{u}_1^T\boldsymbol{X} =c_1x_1+c_2x_2+\cdots+c_px_p$来确定新变量$y_1$，可以看出$y_1$是原变量$x_1,\cdots,x_p$的线性组合。通过同样的方式，可以确定$y_2$等。
通过这种方式，可以将变量进行代换，但是优点并不明显。但是如果在对角矩阵$D$中，前$p$个方差和比其他方差大得多，那么我们就可以近似的把数据当成$p$维数据而不是$p$维的。

多维数据的降维

选取一个重构阈值$t$，选取使得下式成立的最小$p$值

$\frac{\sum^{p$

PCA仅需保留

$P$和样本均值

$M$即可将新样本投影到低维空间。转换过程不可避免的造成数据的损失，但是舍弃这部分信息往往是必要的。一方面降维之后可以使样本的采样密度增大，这也是降维的主要动机，其次可以在一定程度上起到降噪的效果，因为最小的特征值对应的特征向量往往与噪声有关。

算法步骤

input

样本集$[\boldsymbol{X}_1\cdots \boldsymbol{X}_N]$，重构阈值$t$

process

1. 进行去中心化,得到矩阵$B$
2. 计算样本的协方差矩阵$S=BB^T$(常用对输入矩阵的奇异值分解代替特征值分解)
3. 对协方差矩阵进行特征值分解
4. 取出最大的$p$个特征值对应的特征向量 $\boldsymbol{u}_1,\cdots\boldsymbol{u}_{p$

Output

投影矩阵$P$ = [$\boldsymbol{u}_1\cdots\boldsymbol{u}_{p$]

代码实现

特征值分解

from numpy import *

def PCA( data, t ):
    m, n = data.shape
    M = data.sum(1) / n
    M = tile(M,(1,n))
    B = data - M
    D, P = linalg.eig( B*B.T /(n - 1) )
    index = argsort(-D)
    D = D[index]; P = P[index];D[D<0] = 0
    denominator = sum(D); numerator = 0
    for i in range(len(D)):
        numerator += D[i]
        if numerator / denominator >= t:
            return diag(D[:i + 1]),P[:,:i + 1], M
    return diag(D), P, M

奇异值分解

from numpy import *

def PCA( data, k ):
    m, n = data.shape
    M = data.sum(1) / n
    M = tile(M,(1,n))
    B = data - M
    U, S, V = linalg.svd( B*B.T /(n - 1) )
    return U[:,:k+1], S[:,:k+1]