Python 实现主成分分析

2018年5月17日

在统计学里,主成分分析(Principal Component Analysis,简称 PCA)是一种用于分析、简化数据集的技术。在实际问题中,为了能够全面、系统地分析问题,往往会考虑许多影响的因素。比如分析一个地区的发展情况,我们可能会去查询生产总值、人均生产总值、人口数、人均收入、资源环境、景区数目等等众多因素,而这些因素中可能存在一定的相关性,从而存在信息的重叠。人们自然希望通过克服相关性、重叠性,用较少的变量来代替原来就多的变量,且这种代替可以反映原来多个变量的大部分信息,这就是一种“降维”的思想。通过主成分分析可以降低数据的“维数”,又保留数据的大部分信息。

主成分分析的步骤

  1. 对数据进行标准化处理
  2. 计算标准化后的协方差矩阵
  3. 计算协方差矩阵的特征值和特征向量
  4. 对特征值进行排序,确定各成分贡献率,并按大小排列
  5. 依据所需主成份贡献率的大小,选取前 m 个主成份,使这 m 个主成份累积贡献率大于 p

问题实例

假设有一份数据如下,存在一个 Excel 表格中,我们使用 Python 来实现主成份分析,对该数据进行分析。要求选取的主成份涵盖信息量的 85% 以上。

步骤实现

  1. 引入相关包
    我们需要利用 numpy 和 pandas:
import numpy as np
import pandas as pd
  1. 读取数据
filename = 'data.xls'
df = pd.DataFrame(pd.read_excel(filename)) # 打开文件
data_mat = df.drop(df.columns[0], axis=1).values # 删除第一列,即“地区”那一列
  1. 对数据进行标准化处理

在实际应用时,指标的量纲往往不同,所以在主成分计算之前应先消除量纲的影响。消除数据的量纲有很多方法,常用方法是将原始数据标准化,即做如下数据变换:

$$x_{ij}^*=\frac{x_{ij}-\bar{x}_j}{s_j}$$
其中 $$ \bar{x}_j=\frac{1}{n}\sum_{i=1}^{n}x_{ij}, s_j^2=\frac{1}{n-1}\sum_{i=1}^{n}(x_{ij}-\bar{x}_j)^2, j=1,2,…,p $$

这里运用到的标准化处理的方法为 z-score 标准化(zero-mean normalization),这也是 SPSS 中默认的标准化方法,经过处理的数据符合标准正态分布,即均值为0,标准差为1,数据无量纲。

mean_values = np.mean(data_mat, axis=0)
std_mat = (data_mat - mean_values) / np.std(data_mat, axis=0, ddof=1)
  1. 计算协方差矩阵
cov_mat = np.cov(std_mat, rowvar=0)

注:原始数据经标准化后变量的协方差矩阵(Covariance Matrix)$ \sum = (s_{ij})_{p \times p}$,就是原始数据的相关系数矩阵(Correlation Matrix)$R=(r_{ij})_{p \times p} $:

$$s_{ij}=\frac{1}{n-1}\sum_{k=1}^{n} x_{ki}^* x_{kj}^* =\frac{1}{n-1}\sum_{k=1}^{n} \frac{x_{ki} - \bar{x}_i}{\sqrt{\frac{\sum_{t=1}^{n} (x_{ti} - \bar{x}_i)^2 }{n-1}}} \frac{x_{kj} - \bar{x}_j}{\sqrt{\frac{\sum_{t=1}^{n} (x_{tj} - \bar{x}_j)^2 }{n-1}}}=\frac{\sum_{k=1}^{n}(x_{ki}-\bar{x}_i)(x_{kj}-\bar{x}_j)}{\sqrt{\sum_{t=1}^{n}(x_{ti}-\bar{x}_i)^2}\sqrt{\sum_{t=1}^{n}(x_{tj}-\bar{x}_j)^2}}=r_{ij}$$

$$i,j=1,2,…,p$$

我们还可以直接用 numpy.corrcoef() 函数来计算相关系数矩阵。

  1. 计算协方差矩阵的特征值和特征向量
eig_values, eig_vectors = np.linalg.eig(np.mat(cov_mat))

我们求出协方差矩阵 $sum$ 的特征值 $\lambda_1 \geq \lambda_2 \geq … \geq \lambda_q > 0$ 及相应的正交化单位特征向量:
$$
\alpha_1=\begin{pmatrix}
a_{11} \\
a_{21} \\
… \\
a_{p1}
\end{pmatrix}
,
\alpha_2=\begin{pmatrix}
a_{12} \\
a_{22} \\
… \\
a_{p2}
\end{pmatrix}
,

,
\alpha_p=\begin{pmatrix}
a_{1p} \\
a_{2p} \\
… \\
a_{pp}
\end{pmatrix}
$$
则 $X$ 的第 $i$ 个主成份为 $F_i = \alpha_{i}^* X^*, i=1,2,…,p$。其中 $X^*$ 为标准化后的数据。

  1. 根据贡献率 p 选取 m 个主成份

在已确定的全部 $p$ 个主成分中合理选择 $m$ 个来实现最终的评价分析。一般用方差贡献率:
$$ \alpha_i = \frac{\lambda_i}{\sum_{i=1}^{p}\lambda_k}$$
解释主成份 $F_i$ 所反映的信息量的大小, $m$ 的确定以累计贡献率:
$$G(m)=\frac{\sum_{i=1}^{m}\lambda_i}{\sum_{k=1}^{p}\lambda_k}$$
达到足够大(一般在85%以上)为原则。

indices = np.argsort(eig_values)   # 返回特征值从小到大的索引值
indices = indices[::-1]   # 索引值设为从大到小
eig_values_sort = eig_values[indices]
eig_vectors_sort = eig_vectors[:, indices]
rate = eig_values_sort / sum(eig_values_sort)
rate_cumulative = np.cumsum(rate)
  1. 输出
    以上,我们已经实现了 PCA 的算法,下面为方便观察,我们将数据输出到一个 Excel 文件中。详见完整代码。

运行结果

用 SPSS 软件来做一遍,SPSS 输出结果如下:

可以看到,用该算法得到的结果与 SPSS 运行得到的结果一致。

完整源代码

import numpy as np
import pandas as pd


filename = 'data.xls' #要分析的文件名
tofile = 'pca.xls' #保存所产生数据的文件名

df = pd.DataFrame(pd.read_excel(filename)) # 打开文件
data_mat = df.drop(df.columns[0], axis=1).values # 删除第一列,即“地区”那一列

# PCA Algorithm begin
mean_values = np.mean(data_mat, axis=0)
std_mat = (data_mat - mean_values) / np.std(data_mat, axis=0, ddof=1)

cov_mat = np.cov(std_mat, rowvar=0)
eig_values, eig_vectors = np.linalg.eig(np.mat(cov_mat))

indices = np.argsort(eig_values)[::-1]
eig_values = eig_values[indices]
eig_vectors = eig_vectors[:, indices]

m = np.dot(std_mat, eig_vectors)
explained_variance_ratio = eig_values / sum(eig_values) # 计算贡献率
explained_variance_ratio_cumulative = np.cumsum(explained_variance_ratio) # 计算累计贡献率
# PCA Algorithm end


writer = pd.ExcelWriter(tofile)

# 写入相关系数矩阵
data_df = pd.DataFrame(cov_mat)
data_df.columns = df.columns[1:]
data_df.index = df.columns[1:]
data_df.to_excel(writer, 'Correlation Matrix', float_format='%.5f')
writer.save()

eigenvalues = {'component':[], 'value':[], 'difference':[], 'proportion':[], 'cumulative':[]}
for i in range(len(explained_variance_ratio)):
    eigenvalues['component'].append(i + 1)
    eigenvalues['value'].append(eig_values[i])
    if i != len(explained_variance_ratio) - 1:
        eigenvalues['difference'].append(eig_values[i] - eig_values[i + 1])
    else:
        eigenvalues['difference'].append(0)
    eigenvalues['proportion'].append(explained_variance_ratio[i])
    eigenvalues['cumulative'].append(explained_variance_ratio_cumulative[i])

data_df = pd.DataFrame(eigenvalues, columns=eigenvalues.keys())
data_df.to_excel(writer, 'Eigenvalues', float_format='%.5f')
writer.save()
打赏 赞(0)
微信
支付宝
微信二维码图片

微信扫描二维码打赏

支付宝二维码图片

支付宝扫描二维码打赏

One thought on “Python 实现主成分分析

  1. 头像

    您好,我看到您最后算得的主成分是有九个,原数据也是有九个特征,这是说算出的是原数据每一个特征的贡献率吗?不了解这个,有点看不懂。。。谢谢您解答一下疑问好吗?

    Reply

发表评论

电子邮件地址不会被公开。 必填项已用*标注