前几天主要更新了一下机器学习的相关理论，主要介绍了感知机，SVM以及线性判别分析。现在用代码来实现一下其中的模型，一方面对存粹理论的理解，另一方面也提升一下代码的能力。本文就先从线性判别分析开始讲起，不熟悉的可以先移步至线性判别分析（Linear Discriminant Analysis, LDA） - ZhiboZhao - 博客园 (cnblogs.com)对基础知识做一个大概的了解。在代码分析过程中，本文重点从应用入手，只讲API中最常用的参数，能够完成任务即可。
本文代码参考链接：https://github.com/han1057578619/MachineLearning_Zhouzhihua_ProblemSets
一、数据准备
数据集部分我采用周志华《机器学习》书中的 watermelon数据集，数据集前5行如下：
| 编号 | 色泽 | 根蒂 | 敲声 | 纹理 | 脐部 | 触感 | 密度 | 含糖率 | 好瓜 |
|---|---|---|---|---|---|---|---|---|---|
| 1 | 青绿 | 蜷缩 | 浊响 | 清晰 | 凹陷 | 硬滑 | 0.697 | 0.46 | 是 |
| 2 | 乌黑 | 蜷缩 | 沉闷 | 清晰 | 凹陷 | 硬滑 | 0.774 | 0.376 | 是 |
| 3 | 乌黑 | 蜷缩 | 浊响 | 清晰 | 凹陷 | 硬滑 | 0.634 | 0.264 | 是 |
| 4 | 青绿 | 蜷缩 | 沉闷 | 清晰 | 凹陷 | 硬滑 | 0.608 | 0.318 | 是 |
| 5 | 浅白 | 蜷缩 | 浊响 | 清晰 | 凹陷 | 硬滑 | 0.556 | 0.215 | 是 |
1.1 读取数据：
import pandas as pd
data_path = './watermelon3_0_ch.csv'
data = pd.read_csv(data_path).values # 读取数据并转为np.array类型
这里主要运用 pd.read_csv() 进行 .csv 文件的读取，该模块主要用到的参数如下：
pd.read_csv(file_path, sep, header)
其中：file_path 是目标文件的路径；sep 是目标文件中的分隔符，默认 .csv 文件以 ‘,’ 分隔；header 是整数类型的，它的数值决定了读取 .csv 文件时从第几行开始。举个例子：
# header = 0, 默认第0行为表头，从表头往下开始读取
head_0 = pd.read_csv(data_path, header = 0)
# header = 1, 默认第1行为表头，从表头往下开始读取
head_0 = pd.read_csv(data_path, header = 1)
header_0的结果为：
| 编号 | 色泽 | 根蒂 | 敲声 | 纹理 | 脐部 | 触感 | 密度 | 含糖率 | 好瓜 |
|---|---|---|---|---|---|---|---|---|---|
| 1 | 青绿 | 蜷缩 | 浊响 | 清晰 | 凹陷 | 硬滑 | 0.697 | 0.46 | 是 |
| 2 | 乌黑 | 蜷缩 | 沉闷 | 清晰 | 凹陷 | 硬滑 | 0.774 | 0.376 | 是 |
header_1的结果为：
| 1 | 青绿 | 蜷缩 | 浊响 | 清晰 | 凹陷 | 硬滑 | 0.697 | 0.46 | 是 |
|---|---|---|---|---|---|---|---|---|---|
| 2 | 乌黑 | 蜷缩 | 沉闷 | 清晰 | 凹陷 | 硬滑 | 0.774 | 0.376 | 是 |
| 3 | 乌黑 | 蜷缩 | 浊响 | 清晰 | 凹陷 | 硬滑 | 0.634 | 0.264 | 是 |
1.2 对数据进行 "one-hot" 编码
我们以二维线性判别分析为例，只根据 "密度" 和 "含糖量" 来确定是否是好瓜
X = data[:, 7:9].astype(float) # 提取密度和含糖量的数据作为输入特征
y = data[:, 9] # 提取最后一列作为判别类型
y[y == '是'] = 1 # 需要进行one-hot编码，将瓜分类
y[y == '否'] = 0
y = y.astype(int)
'''
以好瓜/坏瓜 来对样本进行分类
'''
pos = y == 1, neg = y == 0 # 分别找到正负样本的位置
X0 = X[neg]， X1 = X[pos] # 以提取正负样本的输入特征
二、线性判别分析
2.1 根据对应模型进行求解
从上一讲中我们得到，线性分类判别模型的最优解为：
\[w = S_{w}^{-1}(u_{0}-u_{1}) \]其中，
\[u_{0} = \dfrac{1}{m} \sum_{i=1}^{m}x_{i},\quad u_{1} = \dfrac{1}{n} \sum_{i=1}^{n}x_{i}\\ S_{w} = \dfrac{1}{m} \sum_{i=1}^{m}(x_{i}-u_{0})(x_{i}-u_{0})^{T} +\dfrac{1}{n} \sum_{i=1}^{n}(x_{i}-u_{1})(x_{i}-u_{1})^{T}\\ \]这里面注意一点，为了更符合人的理解习惯，我们在公式 (3) 中，定义的 \(S_w\) 是单个向量相乘之后求和；但是矩阵形式则更方便被计算机描述，设 $ X_{0} = {x_{1},x_{2},...,x_{m} }^{T}, X_{1} = {x_{1},x_{2},...,x_{n} }^{T}$，由于 \(x_{i} \in R^{p \times 1}\)，因此\(X_{0}, X_{1} \in R^{m \times p}\)，改写成矩阵形式：
\[S_{w} = \dfrac{1}{m} (X_{0}-u_{0})^{T}(X_{0}-u_{0}) + \dfrac{1}{n}(X_{1}-u_{1})^{T}(X_{1}-u_{1}) \]于是，对应代码为：
u0 = X0.mean(0, keepdims=True) # (1, p)
u1 = X1.mean(0, keepdims=True)
sw = np.dot((X0 - u0).T, X0 - u0) + np.dot((X1 - u1).T, X1 - u1)
w = np.dot(np.linalg.inv(sw), (u0 - u1).T).reshape(1, -1) # (1, p)
说明：
mean() 函数在指定维度上求均值，由于 \(X_{0} \in R^{m \times p}\)，所有指定维度为0之后相当于对所有 \(m\) 个样本进行求平均，得到 \(u_{0} \in R^{1\times p}\)
2.2 模型可视化
这一部分代码主要是绘图的一些格式，本文就不多做解释了。
fig, ax = plt.subplots()
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
ax.spines['left'].set_position(('data', 0))
ax.spines['bottom'].set_position(('data', 0))
plt.scatter(X1[:, 0], X1[:, 1], c='k', marker='o', label='good')
plt.scatter(X0[:, 0], X0[:, 1], c='r', marker='x', label='bad')
plt.xlabel('密度', labelpad=1)
plt.ylabel('含糖量')
plt.legend(loc='upper right')
x_tmp = np.linspace(-0.05, 0.15)
y_tmp = x_tmp * w[0, 1] / w[0, 0]
plt.plot(x_tmp, y_tmp, '#808080', linewidth=1)
wu = w / np.linalg.norm(w)
# 正负样板店
X0_project = np.dot(X0, np.dot(wu.T, wu))
plt.scatter(X0_project[:, 0], X0_project[:, 1], c='r', s=15)
for i in range(X0.shape[0]):
plt.plot([X0[i, 0], X0_project[i, 0]], [X0[i, 1], X0_project[i, 1]], '--r', linewidth=1)
X1_project = np.dot(X1, np.dot(wu.T, wu))
plt.scatter(X1_project[:, 0], X1_project[:, 1], c='k', s=15)
for i in range(X1.shape[0]):
plt.plot([X1[i, 0], X1_project[i, 0]], [X1[i, 1], X1_project[i, 1]], '--k', linewidth=1)
# 中心点的投影
u0_project = np.dot(u0, np.dot(wu.T, wu))
plt.scatter(u0_project[:, 0], u0_project[:, 1], c='#FF4500', s=60)
u1_project = np.dot(u1, np.dot(wu.T, wu))
plt.scatter(u1_project[:, 0], u1_project[:, 1], c='#696969', s=60)
ax.annotate(r'u0 投影点',
xy=(u0_project[:, 0], u0_project[:, 1]),
xytext=(u0_project[:, 0] - 0.2, u0_project[:, 1] - 0.1),
size=13,
va="center", ha="left",
arrowprops=dict(arrowstyle="->",
color="k",
)
)
ax.annotate(r'u1 投影点',
xy=(u1_project[:, 0], u1_project[:, 1]),
xytext=(u1_project[:, 0] - 0.1, u1_project[:, 1] + 0.1),
size=13,
va="center", ha="left",
arrowprops=dict(arrowstyle="->",
color="k",
)
)
plt.axis("equal") # 两坐标轴的单位刻度长度保存一致
plt.show()
self.w = w
self.u0 = u0
self.u1 = u1
return self
最终得到的分类结果图如下：