机器学习算法（二）：1 逻辑回归的从零实现（普通实现+多项式特征实现非线性分类+正则化实现三个版本）

---

前言

今天我们开始介绍逻辑回归的从零开始实现代码了，其中内容会包括普通实现、多项式特征实现非线性分类、正则化实现三个版本。相信看完底层实现你对逻辑回归的理解也会上升一个层次。

一、普通实现

1 数据集准备

在训练的初始阶段，我们将要构建一个逻辑回归模型来预测，某个学生是否被大学录取。设想你是大学相关部分的管理者，想通过申请学生两次测试的评分，来决定他们是否被录取。现在你拥有之前申请学生的可以用于训练逻辑回归的训练样本集。对于每一个训练样本，你有他们两次测试的评分和最后是被录取的结果。为了完成这个预测任务，我们准备构建一个可以基于两次测试评分来评估录取可能性的分类模型。
让我们从检查数据开始。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

path = 'ex2data1.txt'
data = pd.read_csv(path, header=None, names=['Exam 1', 'Exam 2', 'Admitted'])
data.head()

输出：

# 可视化一下该二分类数据
fig, ax = plt.subplots(1,1,figsize=(4,3))
ax.scatter(data[data['Admitted']==1]['Exam 1'], data[data['Admitted']==1]['Exam 2'], color = 'lightgreen', marker='o', label='Admitted')
ax.scatter(data[data['Admitted']==0]['Exam 1'], data[data['Admitted']==0]['Exam 2'], color = 'red', marker='x', label='Not Admitted')

plt.xlabel('Exam 1 Score')
plt.ylabel('Exam 2 Score')
plt.legend(loc='upper right')
plt.grid(True)
plt.show()

输出：

看起来在两类间，有一个清晰的决策边界。现在我们需要实现逻辑回归，那样就可以训练一个模型来预测结果。

# 数据准备
X_train = data.iloc[:,0:2].values   # X_train是一个(m,n)的矩阵，m是样本数，n是特征数
y_train = data.iloc[:,2].values     # y_train是一个(m,)的向量
print(f"X_train: {
     
     X_train}")
print(f"y_train: {
     
     y_train}")

输出：

2 逻辑回归模型

$$f_{\mathbf{w},b}(x)=g(\mathbf{w}\cdot \mathbf{x}+b)$$

$$g(z)=\frac{1}{1+e^{-z}}$$

回忆逻辑回归模型，最外层是一个sigmoid函数，因此我们需要先实现sigmoid函数。

def sigmoid(z):
    return 1 / (1 + np.exp(-z))

# 可视化一下sigmoid函数
nums = np.arange(-10, 10, step=1)
fig, ax = plt.subplots(1,1,figsize=(4,3))
ax.plot(nums, sigmoid(nums), color='lightgreen')
plt.grid(True)
plt.show()

输出：

模型实现了，接下来我们需要实现损失函数，以及梯度下降算法。

3 损失函数

$$loss(f_{\mathbf{w},b}({\mathbf{x}}^{(i)}),y^{(i)})=(-y^{(i)}\log \left(f_{\mathbf{w},b}\left({\mathbf{x}}^{(i)}\right)\right)-\left(1-y^{(i)}\right)\log \left(1-f_{\mathbf{w},b}\left({\mathbf{x}}^{(i)}\right)\right)\text{\hspace{1em}(2)}$$

• $f_{\mathbf{w},b}({\mathbf{x}}^{(i)})$ is the model’s prediction, while $y^{(i)}$, which is the actual label

• $f_{\mathbf{w},b}({\mathbf{x}}^{(i)})=g(\mathbf{w}\cdot {\mathbf{x}}^{(\mathbf{i})}+b)$ where function $g$ is the sigmoid function.

  • It might be helpful to first calculate an intermediate variable $z_{\mathbf{w},b}({\mathbf{x}}^{(i)})=\mathbf{w}\cdot {\mathbf{x}}^{(\mathbf{i})}+b=w_0{x}_0^{(i)}+...+w_{n-1}{x}_{n-1}^{(i)}+b$ where $n$ is the number of features, before calculating $f_{\mathbf{w},b}({\mathbf{x}}^{(i)})=g(z_{\mathbf{w},b}({\mathbf{x}}^{(i)}))$
    *
    $J\left(\theta \right)=\frac{1}{m}\sum _{i=1}^m[-y^{(i)}\log \left(h_{\theta }\left(x^{(i)}\right)\right)-\left(1-y^{(i)}\right)\log \left(1-h_{\theta }\left(x^{(i)}\right)\right)]$

def compute_cost_logistic(X, y, w, b):
    """
    Computes cost for logistic regression with numerical stability.
    
    Args:
      X (ndarray (m,n)): Data, m examples with n features
      y (ndarray (m,)) : target values
      w (ndarray (n,)) : model parameters  
      b (scalar)       : model parameter
      
    Returns:
      cost (scalar): cost
    """
    
    m = X.shape[0]
    cost = 0.0
    epsilon = 1e-15  # Small value to ensure numerical stability

    for i in range(m):
        z_i = np.dot(X[i], w) + b              
        f_wb_i = sigmoid(z_i)
        
        # Clip f_wb_i to avoid log(0) and log(1)
        f_wb_i = np.clip(f_wb_i, epsilon, 1 - epsilon)       # 利用clip限制f_wb_i的范围为(epsilon, 1-epsilon)间，避免log(0)和log(1)的情况，这样可以保证数值稳定性 epsilon=1e-15是一个很小的数
        
        cost += -y[i] * np.log(f_wb_i) - (1 - y[i]) * np.log(1 - f_wb_i)
    
    cost = cost / m
    return cost

补充：上面计算损失函数有一个数值不稳定的问题，cost += -y[i]*np.log(f_wb_i) - (1-y[i])*np.log(1-f_wb_i)中对数在训练中出现了 对数里面为0的情况，这样会导致计算出来的cost为nan，在使用逻辑回归计算损失时，log 函数会在 f_wb_i（预测概率）非常接近0或1时引发 divide by zero 警告。这是因为对数函数在接近0时趋向于负无穷大。在计算损失时，这种数值不稳定性会导致非常大的损失值。

因此我们需要对这个进行处理，这里我们可以使用数值稳定性修正的方式来避免，即：
数值稳定性修正：在计算 log 时，将输入值限制在一个非常小的范围之内，以避免数值问题。具体来说，常见的做法是对 f_wb_i 进行裁剪，使其不小于一个很小的数值 ε（如 1e-15），也不大于 1 - ε。这种方法可以有效防止 log(0) 和 log(1) 这样的问题。

接下来就是进行梯度下降了，我们需要先计算损失函数的梯度。

4 计算损失函数的梯度

Recall the gradient descent algorithm utilizes the gradient calculation: