机器学习（二）—逻辑回归理论与代码详解

1. 概要

逻辑回归是一个分类模型，其与线性回归的区别在于最后加了一个sigmoid激活函数将输出转化为概率，从而完成分类的功能

2. 理论

逻辑回归函数是：

$$y=\frac{1}{1+e^{-({\mathbf{w}}^{\mathrm{T}}\mathbf{x}+b)}}=p\text{\hspace{1em}(1)}$$

这里的输出值指的是样本属于类别的概率用$p$表示.

公式1可以变换为（由2推导1还是很直接的）：

$$\ln \frac{y}{1-y}={\mathbf{w}}^{\mathrm{T}}\mathbf{x}+b\text{\hspace{1em}(2)}$$

将 $y$ 视为样本为正例的可能性，$1-y$ 是其反例的可能性，再取对数构成对数几率，通过线性回归模型的预测结果逼近真实标记的对数几率，因此称为”对数几率回归”. 使得每个样本属于其真实标记的概率越大越好。将$\mathbf{w}$与$b$统一表示成$\theta$，则逻辑回归的似然函数：

$$L(\theta )=\prod _{i=1}^{n}P(y_i|{\mathbf{x}}_i;\theta )=\prod _{i=1}^n(h_{\theta }({\mathbf{x}}_i){)}^{y_i}(1-h_{\theta }({\mathbf{x}}_i{)}^{1-y_i})\text{\hspace{1em}(3)}$$

这里的$h_{\theta }(\cdot )$对应着概率输出函数，就是样本属于对应标签的概率值，$n$代表样本数量，这里指的是二分类。

对数似然函数：

$$l(\theta )=\log L(\theta )=\sum _{i=1}^n(y_i\log h_{\theta }({\mathbf{x}}_i)+(1-y_i)\log (1-h_{\theta }({\mathbf{x}}_i)))\text{\hspace{1em}(4)}$$

此问题是一个最大化问题，引入$J(\theta )=-\frac{1}{n}l(\theta )$从而转换为最小化问题。

这里的损失函数是凸函数，可以对一个样本的损失函数求二阶导数判断恒大于等于0，因此可以判断存在全局最小值，可以通过梯度下降法、牛顿法求解参数。参考链接：逻辑回归LR(Logistic Regression)原理以及代码实践 - 简书 (jianshu.com)

损失函数对$\theta$进行求导：

$$\frac{\partial J(\theta )}{\partial \theta }=\frac{1}{n}\sum _{i=1}^n(h_{\theta }({\mathbf{x}}_i)-y_i){\mathbf{x}}_i\text{\hspace{1em}(5)}$$

这里的${\mathbf{x}}_i$组成带有b的数据矩阵$X$形状是[n x (m+1)]，输出与标签的差值形状是[n, 1]，因此这里用矩阵乘法为$\frac{1}{n}{X}^{T}E$，这一步在写代码时特别重要

3. 代码实现

算法流程：

1. 计算逻辑回归输出
2. 计算损失函数对参数的梯度(公式5)
3. 通过梯度下降算法进行参数的更新

注意首先需要将数据矩阵进行填充，偏置项纳入权重向量中，这样就不用单独更新偏置，数据矩阵[n x m], 在第一列和最后一列进行添加，那么权重中的第一项或者最后一项代表的是偏置，这样写代码量减少也更加清晰，输出要明确矩阵的形状，这样进行矩阵乘法时不易出错

import matplotlib.pyplot as plt
import numpy as np

def LogisticRegression_self(X, y, iterations, learning_rate):

    # X:[n, m], y:[n, 1]
    
    X = np.concatenate([X, np.ones((X.shape[0], 1))], axis=1)   # [n, m+1]
    y = y.reshape(len(y), 1)                                    # [n, 1]
    weights = np.zeros((X.shape[1], 1))                         # [m+1, 1]
    errors_cache = []                                           # loss存储列表

    for i in range(iterations): 
				# 1. 计算输出
        out = 1/(1+np.exp(-(X.dot(weights))))                   # X[n, m+1] x Weights[m+1, 1] -> [n, 1] 
        loss = - np.sum(y * np.log(out) + (1 - y) * np.log(1 - out))  # 误差
        errors_cache.append(loss)

        # 2. 计算梯度
				weights_grad = 1/len(X)*(X.T).dot(out - y)      # 梯度

        # 3. 通过梯度进行更新权重
        weights = weights - learning_rate*(weights_grad)  # 更新权重

    return np.array(weights), errors_cache

# 数据准备
X = np.array([[1, 1], [2, 2], [3, 3], [4, 4], [5, 5], [2, 1], [3, 2], [4, 3], [5, 3], [6, 5]])
Y = np.array([1, 1, 1, 1, 1, 0, 0, 0, 0, 0])

# 逻辑回归进行二分类
weights, errors = LogisticRegression_self(X, Y, iterations=3000, learning_rate=0.1)

x1 = np.linspace(0, 7, 3000)
y1 = -(weights[0]*x1 + weights[2]) / weights[1]
plt.plot(x1, y1)                                 # 可视化超平面
plt.plot(np.linspace(0, 7, 3000), errors)        # 可视化损失
plt.scatter(X[:, 0], X[:, 1], c=Y)               # 可视化样本
plt.show()

注意到代码部分，可视化超平面，这里相当于$ax+by+c=0$，移项后$y=-\frac{-ax+c}{b}$，从这个层面理解超平面表达式更容易理解。

实验结果：