当所要分析的样本特征过多时,我们可以采用主成分分析即PCA(principal component analysis)对数据进行降维和可视化。代码引自《python机器学习》

PCA算法及其实现

PCA算法的步骤如下:
1)对原始 d d d维数据集做标准化处理。
2)构造样本的协方差矩阵。
3)计算协方差矩阵的特征值和相应的特征向量。
4)选择与前 k k k个最大特征值对应的特征向量,其中 k k k为新特征空间的维度 ( k ≤ d ) (k\le d) (kd)
5)通过前 k k k个特征向量构建映射矩阵 W \bm{W} W
6)通过映射矩阵 W \bm{W} W d d d维的输入数据集 X \bm{X} X转换到新的 k k k维特征子空间。
加载数据集并将原始数据集划分为训练集和测试集并对数据进行标准化:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
df_wine = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data', header=None)
X, y = df_wine.iloc[:, 1:].values, df_wine.iloc[:, 0].values
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0)
sc = StandardScaler()
X_train_std = sc.fit_transform(X_train)
X_test_std = sc.fit_transform(X_test)

两个特征间的协方差为:
σ j k = 1 n ∑ i = 1 n ( x j ( i ) − μ j ) ( x k ( i ) − μ k ) \sigma_{jk}=\frac{1}{n}\sum_{i=1}^{n}(x_j^{(i)}-\mu_j)(x_k^{(i)}-\mu_k) σjk=n1i=1n(xj(i)μj)(xk(i)μk)
一个包含3个特征的协方差矩阵如下:
Σ = ( σ 11 σ 12 σ 13 σ 21 σ 22 σ 23 σ 31 σ 32 σ 33 ) \Sigma=\left( \begin{array}{ccc} \sigma_{11} &\sigma_{12} & \sigma_{13} \\ \sigma_{21} & \sigma_{22} & \sigma_{23} \\ \sigma_{31}&\sigma_{32}& \sigma_{33} \end{array} \right) Σ=σ11σ21σ31σ12σ22σ32σ13σ23σ33
协方差矩阵的特征向量代表主成分(最大方差方向),其对应特征值大小就决定了特征向量的重要性。
特征值满足:
Σ ν = λ ν \Sigma\bm{\nu}=\lambda\bm{\nu} Σν=λν
可通过如下方法计算协方差矩阵及其特征对:

cov_mat = np.cov(X_train_std.T)# 计算协方差矩阵
eigen_vals, eigen_vecs = np.linalg.eig(cov_mat)# 计算特征值与特征向量
print('\n Eigenvalues \n%s' % eigen_vals)

输出结果为:
Eigenvalues [4.8923083 2.46635032 1.42809973 1.01233462 0.84906459 0.60181514 0.52251546 0.08414846 0.33051429 0.29595018 0.16831254 0.21432212 0.2399553 ]
要实现数据的降维,我们只选择包含最多信息(方差最大)的特征向量组成子集。特征值的大小决定了特征向量的重要性,故需将特征值按降序排列并取出排在前面的 k k k个特征向量。
定义方差贡献率为:
λ j ∑ j = 1 d λ j \frac{\lambda_j}{\sum_{j=1}^d\lambda_j} j=1dλjλj
我们也可以通过以下代码绘制特征值的方差贡献率图像:

cov_mat = np.cov(X_train_std.T)# 计算协方差矩阵
eigen_vals, eigen_vecs = np.linalg.eig(cov_mat)# 计算特征值与特征向量
print('\n Eigenvalues \n%s' % eigen_vals)

tot = sum(eigen_vals)
var_exp = [(i / tot) for i in sorted(eigen_vals, reverse=True)]
cum_var_exp = np.cumsum(var_exp)# 未在参数中提供轴,则将数组展平,并为结果数组计算累计和即为累计方差

plt.bar(range(1, 14), var_exp, alpha=0.5, align='center', label='individual explained variance')
plt.step(range(1, 14), cum_var_exp, where='mid', label='cumulative explained variance')
plt.ylabel('Explained variance ratio')# 方差贡献率
plt.xlabel('Principal components')# 主成分数
plt.legend(loc='best')
plt.show()

在这里插入图片描述
可看出前两个主成分占总方差的近60%,故选取这两个主成分。
通过如下代码,可以获得降维后的葡萄酒特征与种类图:

eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:, i]) for i in range(len(eigen_vals))]
eigen_pairs.sort(reverse=True)# 按特征值降序排列特征对
'''
选取两个对应特征值最大的特征向量
'''
w = np.hstack((eigen_pairs[0][1][:, np.newaxis], eigen_pairs[1][1][:, np.newaxis]))#得到13*2维的映射矩阵W
print('Matrix W:\n', w)
X_train_pca = X_train_std.dot(w) #将124*13维的训练集数据转换到包含两个主成分的空间上
'''
压缩到二维的数据进行可视化展示
'''
colors = ['r', 'b', 'g']
markers = ['s', 'x', 'o']
for l, c, m in zip(np.unique(y_train), colors, markers):
    plt.scatter(X_train_pca[y_train == l, 0], X_train_pca[y_train == l, 1], c=c, label=l, marker=m)
plt.xlabel('PC 1')
plt.ylabel('PC 2')
plt.legend(loc='lower left')
plt.show()

绘制出的图片如下:
在这里插入图片描述
可以使用线性分类器对这三种葡萄酒进行分类。

使用scikit-learn对葡萄酒进行PCA降维分类

scikit-learn提供了PCA类,我们可以使用PCA类和逻辑斯谛回归的方法对葡萄酒进行较方便的分类,代码如下:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.decomposition import PCA


def plot_decision_regions(X, y, classifier, test_idx=None, resolution=0.02):
    # 设置marker generator和color map
    markers = ('s', 'x', 'o', '^', 'v')
    colors = ('red', 'blue', 'lightgreen', 'gray', 'cyan')
    cmap = ListedColormap(colors[:len(np.unique(y))])
    # 画出决策区域
    x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution),
                           np.arange(x2_min, x2_max, resolution))
    Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
    Z = Z.reshape(xx1.shape)
    plt.contourf(xx1, xx2, Z, alpha=0.4, cmap=cmap)
    plt.xlim(xx1.min(), xx1.max())
    plt.ylim = (xx2.min(), xx2.max())
    # 画出所有样本
    X_test, y_test = X[test_idx, :], y[test_idx]
    for idx, cl in enumerate(np.unique(y)):  # np.unique:去除重复数据
        plt.scatter(x=X[y == cl, 0], y=X[y == cl, 1], alpha=0.8, c=cmap(idx), marker=markers[idx], label=cl)
    # 标出需测试的样本
    if test_idx:
        X_test, y_test = X[test_idx, :], y[test_idx]
        plt.scatter(X_test[:, 0], X_test[:, 1], c='black', alpha=0.8, linewidths=1, marker='o', s=10, label='test set')


df_wine = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/wine/wine.data', header=None)
X, y = df_wine.iloc[:, 1:].values, df_wine.iloc[:, 0].values
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0)
sc = StandardScaler()
X_train_std = sc.fit_transform(X_train)
X_test_std = sc.fit_transform(X_test)

pca = PCA(n_components=2)
lr = LogisticRegression()
X_train_pca = pca.fit_transform(X_train_std)
X_test_pca = pca.fit_transform(X_test_std)
lr.fit(X_train_pca, y_train)
plot_decision_regions(X_train_pca, y_train, classifier=lr)
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.legend(loc='lower left')
plt.show()

分类结果如下图,分类效果较好:
在这里插入图片描述

Logo

为开发者提供学习成长、分享交流、生态实践、资源工具等服务,帮助开发者快速成长。

更多推荐