代数与逻辑：作业二 主成分分析法

代数与逻辑：作业二 主成分分析法

一、作业要求

1. 给出PCA的数学推导过程
2. 编程实现PCA
3. 选择公开数据集，编程实现PCA并进行降维。（降低到累计贡献率大于0.9）

二、PCA的数学推导过程

我们的目标是找到传播最大的方向，并将数据点投射到该方向上。

让我们尝试找到一条最大化投影点到原点的距离的线，即最大化投影距离的方差。

这种方法称为方差最大化方法。还有另一种方法可以构造 PCA 的优化函数。

我们考虑 PCA 的另一种方式是，它适合通过我们的数据的最佳线，旨在最小化每个点的投影误差“d”。这种方法称为距离最小化方法。

这两个优化问题虽然看起来不同，但都是相同的。由于 (x^T *x) 项与 u 无关，因此为了最小化函数，我们必须最大化 (u^Tu)²，这与我们的第一个优化问题相同。

假设 X 是我们的数据矩阵，具有 n 个观察值和 d 个特征。

1、协方差矩阵

协方差矩阵（或方差协方差矩阵）是一个方形矩阵，其形状=特征数，其对角线元素是每个特征的方差，非对角线元素是特征之间的协方差。

2、特征值和特征向量

每个特征值对应一个特征向量。所有的特征向量都是正交的。

3、具体步骤

假设我们的数据矩阵 X 有 d 个特征，我们想将它们减少到 k 个特征。

1. 列标准化我们的数据。
2. 找到协方差矩阵。
3. 找出协方差矩阵的所有特征值和特征向量。
4. 那么最大特征值λ1对应的v1就是方差最大的方向，λ2对应的v2就是方差第二大的方向，以此类推。
5. 要获得 k 个特征，我们需要将原始数据矩阵与对应于前 k 个最大特征值的特征向量矩阵相乘。

得到的矩阵是特征减少的矩阵。

三、编程实现PCA

1、导入库

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
1
2
3

2、导入数据集

导入数据集，将数据集分发到 X 和 y 分量进行数据分析。

dataset = pd.read_csv('wine.csv')

X = dataset.iloc[:, 0:13].values
y = dataset.iloc[:, 13].values
1
2
3
4

3、将数据集拆分为训练集和测试集

from sklearn.model_selection import train_test_split
  
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.2, random_state = 0)
1
2
3

4、特征缩放

对训练和测试集进行预处理部分，例如拟合标准量表。

from sklearn.preprocessing import StandardScaler
sc = StandardScaler()
  
X_train = sc.fit_transform(X_train)
X_test = sc.transform(X_test)
1
2
3
4
5

5、应用PCA函数

将 PCA 函数应用到训练和测试集中进行分析。

from sklearn.decomposition import PCA
  
pca = PCA(n_components = 2)
  
X_train = pca.fit_transform(X_train)
X_test = pca.transform(X_test)
  
explained_variance = pca.explained_variance_ratio_
1
2
3
4
5
6
7
8

6、将 Logistic 回归拟合到训练集

from sklearn.linear_model import LogisticRegression  
  
classifier = LogisticRegression(random_state = 0)
classifier.fit(X_train, y_train)
1
2
3
4

7、预测测试集结果

y_pred = classifier.predict(X_test)
1

8、制作混淆矩阵

# Y 和预测值的测试集。
from sklearn.metrics import confusion_matrix
  
cm = confusion_matrix(y_test, y_pred)
1
2
3
4

9、预测训练集结果

# 预测训练集
# 通过散点图得出结果
from matplotlib.colors import ListedColormap
  
X_set, y_set = X_train, y_train
X1, X2 = np.meshgrid(np.arange(start = X_set[:, 0].min() - 1,
                     stop = X_set[:, 0].max() + 1, step = 0.01),
                     np.arange(start = X_set[:, 1].min() - 1,
                     stop = X_set[:, 1].max() + 1, step = 0.01))
  
plt.contourf(X1, X2, classifier.predict(np.array([X1.ravel(),
             X2.ravel()]).T).reshape(X1.shape), alpha = 0.75,
             cmap = ListedColormap(('yellow', 'white', 'aquamarine')))
  
plt.xlim(X1.min(), X1.max())
plt.ylim(X2.min(), X2.max())
  
for i, j in enumerate(np.unique(y_set)):
    plt.scatter(X_set[y_set == j, 0], X_set[y_set == j, 1],
                c = ListedColormap(('red', 'green', 'blue'))(i), label = j)
  
plt.title('Logistic Regression (Training set)')
plt.xlabel('PC1') 
plt.ylabel('PC2') 
plt.legend()
  
# 显示散点图
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

10、可视化测试集结果

# 通过散点图可视化测试集结果
from matplotlib.colors import ListedColormap
  
X_set, y_set = X_test, y_test
  
X1, X2 = np.meshgrid(np.arange(start = X_set[:, 0].min() - 1,
                     stop = X_set[:, 0].max() + 1, step = 0.01),
                     np.arange(start = X_set[:, 1].min() - 1,
                     stop = X_set[:, 1].max() + 1, step = 0.01))
  
plt.contourf(X1, X2, classifier.predict(np.array([X1.ravel(),
             X2.ravel()]).T).reshape(X1.shape), alpha = 0.75,
             cmap = ListedColormap(('yellow', 'white', 'aquamarine'))) 
  
plt.xlim(X1.min(), X1.max())
plt.ylim(X2.min(), X2.max())
  
for i, j in enumerate(np.unique(y_set)):
    plt.scatter(X_set[y_set == j, 0], X_set[y_set == j, 1],
                c = ListedColormap(('red', 'green', 'blue'))(i), label = j)
  
# 散点图的标题
plt.title('Logistic Regression (Test set)') 
plt.xlabel('PC1') 
plt.ylabel('PC2') 
plt.legend()
  
# 显示散点图
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

四、选择公开数据集，编程实现PCA并进行降维

import pandas as pd
import numpy as np
import seaborn as sns

# 导入数据
data= pd.read_csv("./Iris.csv")
data.head()

X = data.iloc[:, 0:4]
y = data.Species

# 第 1 步：缩放数据
from sklearn.preprocessing import StandardScaler
X_scaled = StandardScaler().fit_transform(X)

# 第 2 步：找到协方差矩阵
covar_matrix = (1 / X_scaled.shape[0]) * np.matmul(X_scaled.T,X_scaled)
print (covar_matrix.shape)


# 第 3 步：找到特征值和特征向量
from scipy.linalg import eigh

# eigh 函数按升序返回特征值
# 我们指定前 2 个特征值（在 0、1、2、3 中）
values, vectors = eigh(covar_matrix, eigvals=(2, 3))
print (vectors.shape)
print (vectors)


# 第 4 步：将原始数据投影到特征向量上
pca_components = np.matmul(X_scaled, vectors)
pca_data = pd.DataFrame(np.hstack((pca_components, y.to_numpy().reshape(-1, 1))), columns=["Component 1", "Component 2", "Species"])
pca_data.head()

# 计算累计贡献率
print (values.sum() / eigh(covar_matrix)[0].sum() * 100)

# 绘制主成分
sns.scatterplot(x= "Component 1", y= "Component 2", hue= "Species", data= pca_data)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

运行的结果是：

(4, 4)
(4, 2)
[[-0.00902118 -0.55318314]
 [-0.42093567 -0.51774664]
 [-0.90471285  0.28847469]
 [-0.06501105 -0.58541369]]
93.18211437893089
1
2
3
4
5
6
7