一个简单的PCA实现
具体的计算过程包括以下步骤:
- 数据标准化: 将原始数据进行标准化,使得每个变量具有相同的尺度
- 计算
协方差矩阵
: 计算标准化后的数据的协方差矩阵 - 求解
特征值
和特征向量
: 对协方差矩阵进行特征值分解,得到特征值和对应的特征向量。 - 选择主成分: 按照特征值的大小选择主成分,特征值越大,表示包含的信息越多
- 计算
主成分得分
: 将原始数据投影到选定的主成分上,得到主成分得分 - 计算
贡献率
: 计算每个主成分的贡献率
主成分得分:通过将原始数据点投影到已确定的主成分上得到。投影是通过将数据点乘以主成分的权重向量来完成的
贡献率:每个主成分的特征值除以所有主成分的特征值之和
下面具体看代码
import numpy as np
def pca(data, num_components): # num_components表示需要保留的主成分数量
# 1. 数据标准化
standardized_data = (data - np.mean(data, axis=0)) / np.std(data, axis=0)
# 2. 计算标准化后的数据的协方差矩阵
covariance_matrix = np.cov(standardized_data, rowvar=False)
# 3. 求特征值和特征向量
eigenvalues, eigenvectors = np.linalg.eig(covariance_matrix)
# 4. 选择主成分
sorted_indices = np.argsort(eigenvalues)[::-1]
selected_indices = sorted_indices[:num_components]
sorted_eigenvalues = np.sort(eigenvalues)[::-1]
selected_eigenvalues = sorted_eigenvalues[:num_components]
print("排序后的前{}个特征值:".format(num_components))
print(selected_eigenvalues)
selected_eigenvectors = eigenvectors[:, selected_indices]
print("排序后的前{}个特征向量:".format(num_components))
print(selected_eigenvectors)
# 5. 计算主成分得分
principal_component_scores = np.dot(standardized_data, principal_components)
# 6. 计算贡献率(解释方差)
contribution_rate = eigenvalues[selected_indices] / np.sum(eigenvalues)
return principal_component_scores, principal_components, contribution_rate
在选择主成分时,首先使用np.argsort函数对特征值进行排序,得到按降序排列的特征值的索引;
接着取出前num_components的特征值然后输出,一般取3个;最后选择特征值对应的特征向量,输出的结果每一列对应一个特征向量
举个例子
# 示例数据
data = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
# 设置要保留的主成分数量
num_components = 2
# 进行PCA
scores, components, contribution_rate = pca(data, num_components)
# 打印结果
print("主成分得分:")
print(scores)
print("\n贡献率:")
print(contribution_rate)
输出结果:
再看一个详细的例子:
(1)对这个数据集进行主成分分析
我们将数据存入data之后调用上面的函数,输出的结果如下:
原样本中有8项指标,从结果我们看到,取两个主成分的累计贡献率为75%,虽然这个贡献率不小,但若要求损失信息不超过15%,我们应该取三个主成分,此时的累计贡献率为86.7%
(2)根据前三个主成分的特征向量对其进行解释:
- 第一个主成分:
特征向量[-0.47665026, -0.47280848, -0.42384532, 0.21289287, 0.38845966, 0.35242704, -0.21483468, -0.05503433]
选择绝对值最大的元素,占的权重较大,即前三个特征权值较大,所以可能与年末固定资产净值、职工人数、工业总产值相关
- 第二个主成分:
特征向量[-0.29599107, -0.27789355, -0.37795129, -0.45140781, -0.33094476, -0.40273683, 0.37741499, -0.27273579]
也是前三个特征权值较大,所以可能与年末固定资产净值、职工人数、工业总产值相关
- 第三个主成分:
特征向量[-0.10418975, -0.16298327, -0.1562551, 0.00854434, -0.32113305, -0.14514414, -0.14045919, 0.8911623]
第三、第八个个特征权值较大,所以可能与工业总产值、能源利用效果相关
(3)接下来根据主成分得分对13个行业进行排序和分类
我们可以计算每个样本在所有主成分得分的总和,对总和进行排序
# 计算每个样本在所有主成分上得分的总和
total_scores = np.sum(scores, axis=1)
# 排序样本
sorted_indices = np.argsort(total_scores)[::-1]
# 打印排序结果
print("样本排序结果:")
for i, index in enumerate(sorted_indices):
print("样本 {}: 总得分 = {}".format(index+1, total_scores[index]))
结果:
根据该得分,我们可以将2、3、11、13样本作为一类,6、7、10、12作为一类,剩下的1、4、5、8、9作为一类