线性判别分析之python代码分析

前几天主要更新了一下机器学习的相关理论,主要介绍了感知机,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

最终得到的分类结果图如下:

线性判别分析之python代码分析

上一篇:好用的终端SSH工具SecureCRT Mac破解版


下一篇:支持SSH的终端模拟程序