跳转至

多元 Logistic 回归的推导与 Python 实现

多元 Logistic 回归的目标函数推导,并用 Python 实现 OvO、OvR 和直接构造多元 Logistic 模型的方法。

image-20220928222830418

数学推导

\(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)

image-20220928222618524

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)

image-20220928222704381

直接构造多元逻辑回归模型

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)

image-20220928222738001

比较三种方法的误差率

Python
# 输出各方法的误差率
for method in ['OvO', 'OvR', 'multiclass']:
    print(method+'的分类误差率{:.2%}'.format(globals()['error_rate_'+method]), sep='')

image-20220928222830418

比较三种方法的运行速度

  • OvO 需要构造\(\frac{K\times (K-1)}{2}\)个二元分类器,因此耗时较长(本例耗时 16 秒),但结果也更准确。

  • OvR 需要构造\(K\)个二元分类器,因此耗时较短(本例耗时 3 秒),但结果不如 OvO 准确。

评论