多元 Logistic 回归的推导与 Python 实现¶
多元 Logistic 回归的目标函数推导，并用 Python 实现 OvO、OvR 和直接构造多元 Logistic 模型的方法。
数学推导¶
记 \(x=\left[\begin{array}{l} 1 \\\ x \end{array}\right]\)，即下文的推导中，\(x\) 的含义是在原始特征的上方添加元素“ \(1\) ”。
对数几率模型是线性模型¶
由上述关于 \(x\) 的定义，在多元 Logistic 回归下，有
\[
\log \frac{P(Y=k \mid X=x)}{P(Y=K \mid X=x)}=\beta_k^{\top} x,\quad k=1,2, \cdots, K-1
\]
上式中，\(x\) 和 \(\beta_k\) 均是 \(n+1\) 维列向量。
待求解的 Logistic 模型¶
多元 Logistic 回归的模型参数为
\[
\theta=\lbrace\beta_1, \beta_2, \ldots, \beta_{k-1}\rbrace
\]
当\(k=1,2, \cdots, K-1\)时，
\[
\begin{aligned}
P\left(y_i=k \mid x\right)&=\frac{P\left(y_i=k \mid x\right)}{1}\\\
&=\frac{P\left(y_i=k \mid x\right)}{P\left(y_i=K \mid x\right)+\sum\limits_{l=1}^{K-1} P\left(y_i=l \mid x\right)}\\\
&=\frac{\frac{P\left(y_i=k \mid x\right)}{P\left(y_i=K \mid x\right)}}
{\frac{P\left(y_i=K \mid x\right)}{P\left(y_i=K \mid x\right)}+\frac{\sum\limits_{l=1}^{K-1} P\left(y_i=l \mid x\right)}{P\left(y_i=K \mid x\right)}}\\\
&=\frac{e^{\beta_{k} ^{\top} x}}{1+\sum\limits_{l=1}^{K-1} e^{\beta_{l} ^{\top} x}}
\end{aligned}
\]
当\(k=K\)时，
\[
\begin{aligned}
P\left(y_i=K \mid x\right)=\frac{1}{1+\sum_{l=1}^{K-1} e^{\beta_l ^{\top} x}}
\end{aligned}
\]
似然函数¶
\[
L(\theta)=\prod_{i=1}^n \prod_{k=1}^K P\left(y_i=k \mid x ; \theta\right)^{I\left(y_i=k\right)}
\]
这里的\(I\left(y_i=k\right)\)作为示性函数放在指数部分，方便识别样本的观测值，也方便后续求导，因此是一个十分巧妙的办法。
目标函数（损失函数）¶
寻找最优参数\(\theta\)，使得似然函数最大，即对数似然函数的相反数最小。
\[
\begin{aligned}
l(\theta)=&-\log L(\theta) \\\
=&-\sum_{i=1}^n \sum_{k=1}^K I\left(y_i=k\right) \log P\left(y_{i}=k \mid x_{i; \theta}\right) \\\
=&-\sum_{i=1} ^n \left(\sum_{k=1}^{K-1} I\left(y_{i}=k\right)\left(\beta_k ^{\top} x_i-\log \left(1+\sum_{l=1}^{K-1} e^{\beta_l ^{\top} x_i}\right)\right)\right.\\\
&\left.+I\left(y_i=K\right)\left(-\log \left(1+\sum_{l=1}^{K-1} e^{\beta_l ^{\top} x_i}\right)\right)\right) \\\
=&-\sum_{i=1} ^n \left(\sum_{k=1}^{K-1} I\left(y_i=k\right) \beta_k ^{\top} x_i-\log \left(1+\sum_{l=1}^{K-1} e^{\beta_l ^{\top} x_i}\right) \sum_{k=1}^K I\left(y_i=k\right)\right) \\\
=& \sum_{i=1} ^n \left[-\sum_{k=1}^{K-1} I\left(y_i=k\right) \beta_k ^{\top} x_i+\log \left(1+\sum_{l=1}^{K-1} e^{\beta_l ^{\top} x_i}\right)\right]
\end{aligned}
\]
求梯度下降法的迭代公式¶
先将损失函数对\(\theta\)求导，对于\(k=1,2, \cdots, K-1\)，有
\[
\begin{aligned}
\frac{\partial l(\theta)}{\partial \beta_k} &=\sum_{i=1} ^n \left[-I\left(y_i=k\right) x_i+\frac{e^{\beta_k ^{\top} x_i} x_i}{1+\sum\limits_{l=1}^{K-1} e^{\beta_l ^{\top} x_i}}\right] \\\
&=\sum_{i=1} ^n \left(P\left(y_i=k \mid x\right)-I\left(y_i=k\right)\right) x_i
\end{aligned}
\]
因此，梯度下降法的迭代公式为
\[
\beta_{k}^{'}=\beta_k-\sum_{i=1}^n\left(P\left(y_i=k \mid x\right)-I\left(y_i=k\right)\right) x_i
\]
Python 实现多元逻辑回归¶
导入相关包¶
Python
# 导入数据处理相关的包
import pandas as pd
import numpy as np
from scipy import stats
# 导入逻辑回归的包
from sklearn.linear_model import LogisticRegression as LR
from sklearn.preprocessing import LabelEncoder
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score
# 导入绘制循环进度条的包
from tqdm import tqdm
# 忽略警告
import warnings
warnings.filterwarnings("ignore")
导入数据¶
Python
# 导入数据
data = pd.read_csv("letter_recognition.csv")
# 划分特征和标签
X = data.iloc[:, 1:]
y = data.iloc[:, :1]
# 为字母进行编码
y = LabelEncoder().fit_transform(y.values.ravel())
# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.4, shuffle=False)
# 提取唯一的标签
labels = np.unique(y_train)
One vs. One 算法¶
Python
# OvO
count = 0
# 将预测值初始化为 0
y_pred = np.zeros((len(y_test), 1))
pbar = tqdm(range(len(labels)))
for i in pbar:
for j in range(i + 1, len(labels)):
# 输出进度条
pbar.set_description("正在构建第{}_{}分类器".format(labels[i] + 1, labels[j] + 1))
count += 1
# 只使用属于第 i 类和第 j 类的样本训练模型
y_train_ = np.where(
np.logical_or(y_train == labels[i], y_train == labels[j]), y_train, np.NAN
)
X_train_ = X_train.iloc[~np.isnan(y_train_)]
# 二元逻辑回归
globals()["LR_{}_{}".format(labels[i], labels[j])] = LR().fit(
X_train_, y_train_[~np.isnan(y_train_)]
)
# 测试集预测
globals()["y_pred_{}_{}".format(labels[i], labels[j])] = globals()[
"LR_{}_{}".format(labels[i], labels[j])
].predict(X_test)
# 合并预测结果
y_pred = np.append(
y_pred,
globals()["y_pred_{}_{}".format(labels[i], labels[j])].reshape(-1, 1),
axis=1,
)
print("一共构建了{}个分类器".format(count))
# 删除第一列的初始零值
y_pred = np.delete(y_pred, 0, axis=1)
# 以所有分类器的分类结果的众数作为最终的预测结果
y_pred = np.squeeze(stats.mode(y_pred, axis=1)[0])
error_rate_OvO = 1 - accuracy_score(y_test, y_pred)
# 查看预测结果
pd.DataFrame(
np.append(y_pred.reshape(-1, 1), y_test.reshape(-1, 1), axis=1),
columns=["预测值", "真实值"],
dtype=int,
)
One vs. Rest 算法¶
Python
# OvR
count = 0
# 将预测值初始化为 0
y_pred = np.zeros((len(y_test), 1))
pbar = tqdm(range(len(labels)))
for i in pbar:
# 输出进度条
pbar.set_description("正在构建第{}分类器".format(labels[i] + 1))
count += 1
# 只使用属于第 i 类和第 j 类的样本训练模型
y_train_ = np.where(y_train == labels[i], y_train, -1)
# 二元逻辑回归
globals()["LR_{}".format(labels[i])] = LR().fit(X_train, y_train_)
# 测试集预测
globals()["y_pred_{}".format(labels[i])] = globals()[
"LR_{}".format(labels[i])
].predict_proba(X_test)[:, 1]
# 合并预测结果
y_pred = np.append(
y_pred, globals()["y_pred_{}".format(labels[i])].reshape(-1, 1), axis=1
)
print("一共构建了{}个分类器".format(count))
# 删除第一列的初始零值
y_pred = np.delete(y_pred, 0, axis=1)
# 以所有分类器的分类概率的最大值对应的值作为最终的预测结果
y_pred = np.argmax(y_pred, axis=1)
error_rate_OvR = 1 - accuracy_score(y_test, y_pred)
# 查看预测结果
pd.DataFrame(
np.append(y_pred.reshape(-1, 1), y_test.reshape(-1, 1), axis=1),
columns=["预测值", "真实值"],
dtype=int,
)
直接构造多元逻辑回归模型¶
Python
# 直接构造多元逻辑回归模型
# 多元逻辑回归
LR_multiclass = LR(multi_class="multinomial").fit(X_train, y_train)
# 测试集预测
y_pred_multiclass = LR_multiclass.predict(X_test)
error_rate_multiclass = 1 - accuracy_score(y_test, y_pred_multiclass)
# 查看预测结果
pd.DataFrame(
np.append(y_pred_multiclass.reshape(-1, 1), y_test.reshape(-1, 1), axis=1),
columns=["预测值", "真实值"],
dtype=int,
)
比较三种方法的误差率¶
Python
# 输出各方法的误差率
for method in ["OvO", "OvR", "multiclass"]:
print(method + "的分类误差率{:.2%}".format(globals()["error_rate_" + method]), sep="")
比较三种方法的运行速度¶
-
OvO 需要构造\(\frac{K\times (K-1)}{2}\)个二元分类器，因此耗时较长（本例耗时 16 秒），但结果也更准确。
-
OvR 需要构造\(K\)个二元分类器，因此耗时较短（本例耗时 3 秒），但结果不如 OvO 准确。