s03 逻辑回归 — demo.py 代码详解

运行方式

bash

cd s03_logistic_regression/code
python demo.py

代码逐段详解

第1步：导入库 — 每个库是做什么的

python

import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression as SklearnLR

os：操作文件路径，创建 images/ 目录
numpy：数值计算核心。关键用法：np.exp() 指数函数（Sigmoid 和 Softmax 的核心），np.clip() 裁剪数值范围防止溢出，np.log() 对数函数（交叉熵损失），np.argmax() 取最大值索引，np.unique() 统计类别数
matplotlib：绘图，包括函数曲线、散点图、决策边界热力图、等高线图
load_iris（sklearn）：加载经典的 Iris（鸢尾花）数据集——150 个样本，4 个特征，3 个类别
train_test_split（sklearn）：按比例（80/20）随机划分训练集和测试集
sklearn.linear_model.LogisticRegression：sklearn 的实现，用作基准对比

第2步：Sigmoid 函数 — 从实数到概率

python

def sigmoid(z: np.ndarray) -> np.ndarray:
    z_clipped = np.clip(z, -500, 500)
    return 1.0 / (1.0 + np.exp(-z_clipped))

Sigmoid 函数的数学形式：

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

Sigmoid 的核心性质：

值域： $σ (z) \in (0, 1)$ ，天然适合解释为概率。当 $z \to + \infty$ 时 $σ \to 1$ ，当 $z \to - \infty$ 时 $σ \to 0$
对称性： $σ (- z) = 1 - σ (z)$ ，得分取反，概率互补
导数： $σ^{'} (z) = σ (z) (1 - σ (z))$ ，可以用自身表达——这让反向传播计算异常简便
中心点： $σ (0) = 0.5$ ，对应决策边界 $w \cdot x + b = 0$

数值稳定技巧：np.clip(z, -500, 500) 将 $z$ 限制在 $[- 500, 500]$ 内。当 $z > 500$ 时 $e^{- z} \approx 0$ （下溢出），当 $z < - 500$ 时 $e^{- z} \approx \infty$ （上溢出）。裁剪确保计算稳定，而 $[- 500, 500]$ 已经足够覆盖所有实际应用场景（ $σ (500) \approx 1$ ， $σ (- 500) \approx 0$ ）。

第3步：Softmax 函数 — 从得分到概率分布

python

def softmax(z: np.ndarray) -> np.ndarray:
    z_stable = z - np.max(z, axis=1, keepdims=True)
    exp_z = np.exp(z_stable)
    return exp_z / np.sum(exp_z, axis=1, keepdims=True)

Softmax 将 $K$ 个原始得分 $[z_{1}, \dots, z_{K}]$ 转换为概率分布：

softmax (z_{k}) = \frac{e^{z_{k}}}{\sum_{j = 1}^{K} e^{z_{j}}}

三个关键性质：

非负性：每个输出 $\in (0, 1)$
归一性：所有输出之和 $= 1$ （构成一个合法的概率分布）
保序性：如果 $z_{i} > z_{j}$ ，则 $softmax (z_{i}) > softmax (z_{j})$

数值稳定技巧：z - np.max(z, axis=1, keepdims=True) 在每行（每个样本）减去该行的最大值。这等价于分子分母同除 $e^{max (z)}$ ，不改变概率值但避免了 $e^{z}$ 上溢出。axis=1 表示沿列方向（对每个样本的所有类别做 max），keepdims=True 保持维度以便广播。

第4步：二分类逻辑回归模型

4.1 模型定义

逻辑回归 = 线性模型 + Sigmoid：

P (y = 1 ∣ x) = \hat{y} = σ (w^{T} x + b) = \frac{1}{1 + e^{- (w^{T} x + b)}}

这个值被解释为"给定输入 $x$ ，样本属于正类的概率"。逻辑回归本质上就是一个没有隐藏层的神经网络。

预测概率 _predict_proba()：

python

def _predict_proba(self, X):
    z = X @ self.w + self.b    # 线性组合 (n,)
    return sigmoid(z)           # 通过 Sigmoid 得到概率 (n,)

X @ self.w 是矩阵乘法 $X_{n \times d} \cdot w_{d}$ ，结果是 $n$ 个得分值。然后逐元素通过 Sigmoid 转化为概率。

类别预测 predict()：

python

def predict(self, X, threshold=0.5):
    proba = self._predict_proba(X)
    return (proba >= threshold).astype(int)

当预测概率 $\geq 0.5$ （即 $w^{T} x + b \geq 0$ ）时判为正类（1），否则判为负类（0）。

4.2 损失函数：二元交叉熵

为什么不继续用 MSE？两个致命问题：

Sigmoid 的饱和特性导致梯度消失（当预测值接近 0 或 1 且标签相反时，MSE 的梯度极小）
MSE + Sigmoid 在参数空间中是非凸的，存在多个局部最小值

解决方案是二元交叉熵（Binary Cross-Entropy）：

J (w, b) = - \frac{1}{n} \sum_{i = 1}^{n} [y_{i} \log ({\hat{y}}_{i}) + (1 - y_{i}) \log (1 - {\hat{y}}_{i})]

这个公式设计得非常精妙：当 $y = 1$ 时只有第一项起作用（ $- \log (\hat{y})$ ），预测概率越高损失越小；当 $y = 0$ 时只有第二项起作用（ $- \log (1 - \hat{y})$ ），预测概率越低损失越小。

代码实现：

python

def _compute_loss(self, X, y):
    y_pred = self._predict_proba(X)
    n = len(y)
    eps = 1e-15
    y_pred = np.clip(y_pred, eps, 1 - eps)  # 防止 log(0)
    loss = -(1.0 / n) * np.sum(
        y * np.log(y_pred) + (1 - y) * np.log(1 - y_pred)
    )
    return loss

np.clip(y_pred, eps, 1 - eps) 将预测概率限制在 $[10^{- 15}, 1 - 10^{- 15}]$ 之间，防止 $\log (0) = - \infty$ 导致数值崩溃。 $10^{- 15}$ 已经足够小，对实际损失值几乎没有影响。

4.3 梯度：数学之美的巅峰

Sigmoid + 交叉熵组合的梯度有一个惊人简洁的形式。经过数学推导（链式法则 + Sigmoid 导数性质）：

\frac{\partial J}{\partial z} = \hat{y} - y

这被称为"黄金组合"——损失对原始得分 $z$ 的梯度就等于预测误差！利用链式法则传播到参数：

\frac{\partial J}{\partial w} = \frac{1}{n} X^{T} (\hat{y} - y)

\frac{\partial J}{\partial b} = \frac{1}{n} \sum ({\hat{y}}_{i} - y_{i})

代码实现：

python

def _compute_gradients(self, X, y):
    y_pred = self._predict_proba(X)
    n = len(y)
    errors = y_pred - y                           # (n,) 预测误差
    dw = (1.0 / n) * (X.T @ errors)               # (d,) 权重梯度
    db = (1.0 / n) * np.sum(errors)               # 标量 偏置梯度
    return dw, db

X.T @ errors 是矩阵乘法： $X_{(d \times n)}^{T} \cdot {errors}_{(n)}$ ，每行是将所有样本的误差按对应特征加权求和。这个梯度形式与线性回归 + MSE 的梯度完全一致——唯一的区别在于 $\hat{y}$ 的计算方式不同（逻辑回归中多了 Sigmoid）。

4.4 训练循环

python

def fit(self, X, y, verbose=True):
    self.w = np.random.randn(n_features) * 0.01
    self.b = 0.0
    for epoch in range(self.max_epochs):
        loss = self._compute_loss(X, y)       # 前向 + 损失
        dw, db = self._compute_gradients(X, y) # 反向传播
        self.w -= self.learning_rate * dw       # 参数更新
        self.b -= self.learning_rate * db

与线性回归的训练循环结构完全相同——这揭示了机器学习的一个统一框架：

前向传播：计算 $\hat{y} = f_{θ} (x)$
计算损失： $J (\hat{y}, y)$
反向传播：求 $\nabla_{θ} J$
更新参数： $θ \leftarrow θ - η \nabla_{θ} J$

不同模型只是换用了不同的 $f_{θ}$ 和 $J$ ，但四步循环的结构不变。

第5步：多分类 Softmax 回归

二分类 $\to$ 多分类的推广：

权重：从向量 $w_{(d,)}$ 变为矩阵 $W_{(d \times K)}$ ，每个类别有自己的权重向量
偏置：从标量 $b$ 变为向量 $b_{(K,)}$ ，每个类别有自己的偏置
激活：从 Sigmoid 变为 Softmax
损失：从二元交叉熵变为多分类交叉熵

多分类交叉熵损失：

J = - \frac{1}{n} \sum_{i = 1}^{n} \sum_{k = 1}^{K} y_{i k} \log ({\hat{y}}_{i k})

由于真实标签是 one-hot 编码（只有正确类别的位置为 1，其余为 0），实际计算时只需取正确类别位置的 $\log$ 概率：

python

def _compute_loss(self, X, y):
    proba = self._predict_proba(X)          # (n, K)
    proba = np.clip(proba, eps, 1 - eps)
    loss = -(1.0 / n) * np.sum(np.log(proba[np.arange(n), y]))
    return loss

proba[np.arange(n), y] 是高级索引（fancy indexing）：对每个样本 $i$ ，取第 $i$ 行、第 $y_{i}$ 列的概率值——即模型对该样本正确类别的预测概率。

One-hot 编码在梯度计算中扮演关键角色：

python

y_onehot = np.zeros((n, self.n_classes))
y_onehot[np.arange(n), y] = 1
errors = proba - y_onehot

y_onehot 中，每行只有一个位置为 1（真实类别），其余为 0。errors = proba - y_onehot 的结果是：正确类别位置为 ${\hat{y}}_{k} - 1$ （负值），其他位置为 ${\hat{y}}_{k}$ （正值）。梯度方向就是"增大正确类别的概率，减小错误类别的概率"。

第6步：可视化

6.1 Sigmoid 函数曲线

在 $[- 8, 8]$ 区间绘制 Sigmoid 曲线，关键标注：

红色虚线标注 $z = 0$ 和 $σ = 0.5$ ——这是决策边界的位置
灰色虚线标注渐近线 $y = 0$ 和 $y = 1$ ——Sigmoid 的值域下界和上界
绿色区域标注 $z > 0$ 的"正类区域"，红色区域标注 $z < 0$ 的"负类区域"

6.2 决策边界与概率热力图

这是逻辑回归最直观的可视化。代码使用 np.meshgrid() 在二维平面上创建 $300 \times 300$ 的网格，对每个网格点计算 $P (y = 1 ∣ x)$ ，然后用 contourf() 填充颜色：

蓝色区域：模型输出低概率，倾向于预测负类
红色区域：模型输出高概率，倾向于预测正类
绿色轮廓线： $σ = 0.5$ 的等概率线，即决策边界 $w^{T} x + b = 0$

颜色渐变展示了模型的"置信度"——远离决策边界的区域颜色更深（更确信），决策边界附近的区域颜色较浅（不确定）。

第7步：模型评估 — 混淆矩阵

对于分类问题，仅看准确率可能不够。代码计算了**混淆矩阵（Confusion Matrix）**的四个元素：

	预测正类 (1)	预测负类 (0)
真实正类 (1)	TP (真正例)	FN (假负例)
真实负类 (0)	FP (假正例)	TN (真负例)

准确率： $\frac{T P + T N}{T P + T N + F P + F N}$ ，所有预测中正确的比例
从混淆矩阵可以进一步计算精确率 $P = \frac{T P}{T P + F P}$ 、召回率 $R = \frac{T P}{T P + F N}$ 、F1 分数等更细致的指标

关键概念速查表

概念	数学形式	代码位置	关键说明
Sigmoid	$σ (z) = \frac{1}{1 + e^{- z}}$	`sigmoid()`	实数→(0,1)概率映射
Softmax	$\frac{e^{z_{k}}}{\sum e^{z_{j}}}$	`softmax()`	得分→概率分布
交叉熵 (二元)	$- [y \log \hat{y} + (1 - y) \log (1 - \hat{y})]$	`_compute_loss()`	分类问题的标准损失
黄金梯度	$\partial J / \partial z = \hat{y} - y$	`_compute_gradients()`	Sigmoid导数被约掉
决策边界	$w^{T} x + b = 0$	可视化	$σ = 0.5$ 的等概率线
One-hot 编码	正确类别=1，其余=0	`_compute_gradients()`	多分类梯度的关键
混淆矩阵	TP/TN/FP/FN	`main()`	分类问题精细化评估

完整代码

# -*- coding: utf-8 -*-
"""
===============================================================================
s03_logistic_regression/code/demo.py — 逻辑回归从零实现
===============================================================================
本演示从零实现逻辑回归，涵盖二分类和多分类（Softmax），使用 Iris 数据集。
内容包括 Sigmoid 函数、交叉熵损失、梯度下降、决策边界可视化等。

通过本演示，你将理解：
  1. Sigmoid 函数如何将实数映射为 (0,1) 的概率
  2. 交叉熵损失为何是分类问题的标准选择
  3. 逻辑回归的梯度为何如此简洁：∂L/∂z = ŷ - y
  4. 决策边界的几何含义
  5. Softmax 如何将二分类推广到多分类
  6. 概率热力图——展示模型的「置信度」分布

作者：learn-ai 项目
日期：2025
===============================================================================
"""

import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from sklearn.datasets import load_iris  # 加载经典的 Iris 数据集
from sklearn.model_selection import train_test_split  # 划分训练/测试集
from sklearn.linear_model import LogisticRegression as SklearnLR  # sklearn 实现（对比用）
matplotlib.rcParams['axes.unicode_minus'] = False

# 图片保存目录：固定为本章节的 images/ 目录（相对于本脚本的 ../images/）
_SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
_IMAGES_DIR = os.path.join(_SCRIPT_DIR, '..', 'images')
os.makedirs(_IMAGES_DIR, exist_ok=True)


# ============================================================================
# 第一部分：核心函数
# ============================================================================

def sigmoid(z: np.ndarray) -> np.ndarray:
    """
    Sigmoid 激活函数。

    σ(z) = 1 / (1 + e^(-z))
    将任意实数 z 映射到 (0, 1) 区间，解释为概率。

    参数:
        z: np.ndarray, 输入值（可以是标量、向量或矩阵）

    返回:
        np.ndarray, Sigmoid 输出，形状与 z 相同
    """
    # np.clip 防止数值溢出：当 z 很大时 e^(-z) ≈ 0，z 很小时 e^(-z) ≈ inf
    z_clipped = np.clip(z, -500, 500)  # 将 z 限制在 [-500, 500] 内
    return 1.0 / (1.0 + np.exp(-z_clipped))  # 计算 Sigmoid 函数


def softmax(z: np.ndarray) -> np.ndarray:
    """
    Softmax 函数：将 K 个原始得分转化为概率分布。

    softmax(z_k) = e^{z_k} / Σ_j e^{z_j}

    技巧：先减去最大值 (z - max(z)) 以提高数值稳定性，
    因为 e^{z_i - max(z)} 的最大值为 1，不会溢出。

    参数:
        z: np.ndarray, 形状 (n_samples, n_classes)，原始得分矩阵

    返回:
        np.ndarray, 形状 (n_samples, n_classes)，概率分布矩阵
    """
    # 数值稳定技巧：减去每行的最大值
    z_stable = z - np.max(z, axis=1, keepdims=True)  # 保持形状 (n, 1) 以便广播
    exp_z = np.exp(z_stable)  # 计算指数
    return exp_z / np.sum(exp_z, axis=1, keepdims=True)  # 归一化为概率


# ============================================================================
# 第二部分：二分类逻辑回归
# ============================================================================

class LogisticRegression:
    """
    二分类逻辑回归（使用梯度下降法）。

    模型: P(y=1|x) = σ(w^T x + b)，其中 σ 是 Sigmoid 函数
    损失: 二元交叉熵 J = -(1/n) Σ [y_i log(ŷ_i) + (1-y_i) log(1-ŷ_i)]
    梯度: ∂J/∂w = (1/n) X^T (ŷ - y), ∂J/∂b = (1/n) Σ (ŷ_i - y_i)

    注意梯度的简洁性：∂L/∂z = ŷ - y。这正是 Sigmoid + 交叉熵「黄金组合」
    的数学之美——Sigmoid 的导数项在链式法则中被约掉。

    属性:
        w: np.ndarray, 权重向量，形状 (n_features,)
        b: float, 偏置
        loss_history: list, 训练过程中的损失记录
    """

    def __init__(self, learning_rate: float = 0.1, max_epochs: int = 1000,
                 tolerance: float = 1e-6):
        """
        初始化逻辑回归模型。

        参数:
            learning_rate: float, 学习率 η
            max_epochs: int, 最大训练轮数
            tolerance: float, 收敛容忍度
        """
        self.learning_rate = learning_rate
        self.max_epochs = max_epochs
        self.tolerance = tolerance
        self.w = None  # 权重向量
        self.b = None  # 偏置
        self.loss_history = []  # 损失记录

    def _predict_proba(self, X: np.ndarray) -> np.ndarray:
        """
        计算样本属于正类的概率 P(y=1|x)。

        ŷ = σ(X @ w + b)

        参数:
            X: np.ndarray, 形状 (n_samples, n_features)

        返回:
            np.ndarray, 形状 (n_samples,)，每个样本属于正类的概率
        """
        z = X @ self.w + self.b  # 线性组合，形状 (n_samples,)
        return sigmoid(z)  # 通过 Sigmoid 得到概率

    def predict(self, X: np.ndarray, threshold: float = 0.5) -> np.ndarray:
        """
        对样本进行类别预测。

        ŷ >= threshold 预测为正类 (1)，否则为负类 (0)。

        参数:
            X: np.ndarray, 形状 (n_samples, n_features)
            threshold: float, 分类阈值，默认 0.5

        返回:
            np.ndarray, 形状 (n_samples,)，预测的类别标签 (0 或 1)
        """
        proba = self._predict_proba(X)  # 计算概率
        return (proba >= threshold).astype(int)  # 转换为 0/1 标签

    def _compute_loss(self, X: np.ndarray, y: np.ndarray) -> float:
        """
        计算二元交叉熵损失。

        J = -(1/n) Σ [y_i log(ŷ_i) + (1 - y_i) log(1 - ŷ_i)]

        为了防止 log(0)，对 ŷ 加一个极小值 eps。

        参数:
            X: np.ndarray, 特征矩阵
            y: np.ndarray, 真实标签 (0 或 1)

        返回:
            float, 交叉熵损失
        """
        y_pred = self._predict_proba(X)  # 预测概率
        n = len(y)  # 样本数
        eps = 1e-15  # 小常数，防止 log(0)
        # 限制概率在 [eps, 1-eps] 之间以保证 log 的数值稳定
        y_pred = np.clip(y_pred, eps, 1 - eps)
        # 交叉熵公式
        loss = -(1.0 / n) * np.sum(
            y * np.log(y_pred) + (1 - y) * np.log(1 - y_pred)
        )
        return loss

    def _compute_gradients(self, X: np.ndarray, y: np.ndarray):
        """
        计算损失对 w 和 b 的梯度。

        ∂J/∂w = (1/n) X^T (ŷ - y)
        ∂J/∂b = (1/n) Σ (ŷ_i - y_i)

        这是 Sigmoid + 交叉熵组合的优美结果——梯度等于「预测误差」的加权和。

        参数:
            X: np.ndarray, 特征矩阵 (n_samples, n_features)
            y: np.ndarray, 真实标签 (n_samples,)

        返回:
            dw: np.ndarray, ∂J/∂w, 形状 (n_features,)
            db: float, ∂J/∂b
        """
        y_pred = self._predict_proba(X)  # 预测概率
        n = len(y)  # 样本数
        errors = y_pred - y  # 预测误差 (n_samples,)

        # ∂J/∂w = (1/n) X^T @ errors
        dw = (1.0 / n) * (X.T @ errors)  # 形状 (n_features,)
        # ∂J/∂b = (1/n) Σ errors
        db = (1.0 / n) * np.sum(errors)  # 标量

        return dw, db

    def fit(self, X: np.ndarray, y: np.ndarray, verbose: bool = True):
        """
        使用梯度下降法训练逻辑回归模型。

        算法步骤:
            1. 随机初始化 w 和 b
            2. 重复:
                a. 前向计算: ŷ = σ(Xw + b)
                b. 计算损失: J = cross_entropy(ŷ, y)
                c. 反向传播: 计算 ∂J/∂w 和 ∂J/∂b
                d. 参数更新: w ← w - η·∂J/∂w, b ← b - η·∂J/∂b
            3. 直到收敛或达到最大轮数

        参数:
            X: np.ndarray, 特征矩阵 (n_samples, n_features)
            y: np.ndarray, 真实标签 (n_samples,)，取值为 0 或 1
            verbose: bool, 是否打印训练日志
        """
        n_samples, n_features = X.shape

        # 初始化参数：使用小随机数，偏置初始化为 0
        self.w = np.random.randn(n_features) * 0.01  # 权重：小随机数
        self.b = 0.0  # 偏置：初始化为 0
        self.loss_history = []

        for epoch in range(self.max_epochs):
            # 前向计算 + 损失
            loss = self._compute_loss(X, y)
            self.loss_history.append(loss)

            # 计算梯度
            dw, db = self._compute_gradients(X, y)

            # 梯度下降更新参数
            self.w -= self.learning_rate * dw  # w ← w - η·∂J/∂w
            self.b -= self.learning_rate * db  # b ← b - η·∂J/∂b

            # 检查收敛
            if len(self.loss_history) > 1:
                if abs(self.loss_history[-2] - loss) < self.tolerance:
                    if verbose:
                        print(f"第 {epoch + 1} 轮收敛！损失变化 < {self.tolerance}")
                    break

            if verbose and (epoch + 1) % 100 == 0:
                print(f"Epoch {epoch + 1:4d}: loss={loss:.6f}")

        if verbose:
            print(f"\n训练完成，共 {epoch + 1} 轮，最终损失: {self.loss_history[-1]:.6f}")


# ============================================================================
# 第三部分：多分类逻辑回归（Softmax 回归）
# ============================================================================

class SoftmaxRegression:
    """
    多分类逻辑回归（Softmax 回归）。

    模型: P(y=k|x) = softmax(X @ W + b)_k
    损失: 多分类交叉熵 J = -(1/n) Σ Σ y_{ik} log(ŷ_{ik})
    梯度: ∂J/∂W = (1/n) X^T (ŷ - y_onehot)

    其中 y_onehot 是真实标签的 one-hot 编码。

    属性:
        W: np.ndarray, 权重矩阵，形状 (n_features, n_classes)
        b: np.ndarray, 偏置向量，形状 (n_classes,)
        loss_history: list, 损失记录
    """

    def __init__(self, learning_rate: float = 0.1, max_epochs: int = 2000,
                 tolerance: float = 1e-6):
        """初始化多分类逻辑回归。"""
        self.learning_rate = learning_rate
        self.max_epochs = max_epochs
        self.tolerance = tolerance
        self.W = None  # 权重矩阵 (n_features, n_classes)
        self.b = None  # 偏置向量 (n_classes,)
        self.loss_history = []
        self.n_classes = None

    def _predict_proba(self, X: np.ndarray) -> np.ndarray:
        """
        计算每个样本属于每个类别的概率。

        ŷ = softmax(X @ W + b)

        参数:
            X: np.ndarray, (n_samples, n_features)

        返回:
            np.ndarray, (n_samples, n_classes)，每行是一个概率分布
        """
        z = X @ self.W + self.b  # 线性组合，形状 (n_samples, n_classes)
        return softmax(z)  # Softmax 归一化

    def predict(self, X: np.ndarray) -> np.ndarray:
        """
        预测类别（取概率最大的类别）。

        参数:
            X: np.ndarray, (n_samples, n_features)

        返回:
            np.ndarray, (n_samples,)，类别标签
        """
        proba = self._predict_proba(X)  # 计算概率
        return np.argmax(proba, axis=1)  # 取每行最大概率的索引

    def _compute_loss(self, X: np.ndarray, y: np.ndarray) -> float:
        """
        计算多分类交叉熵损失。

        J = -(1/n) Σ_i Σ_k y_{ik} log(ŷ_{ik})

        参数:
            X: np.ndarray, 特征矩阵
            y: np.ndarray, 类别标签 (0, 1, ..., K-1)

        返回:
            float, 交叉熵损失
        """
        proba = self._predict_proba(X)  # 预测概率矩阵
        n = len(y)  # 样本数
        eps = 1e-15  # 防止 log(0)
        proba = np.clip(proba, eps, 1 - eps)  # 数值稳定

        # 交叉熵：只取真实类别位置的对数概率
        # proba[np.arange(n), y] 取出每个样本正确类别的预测概率
        loss = -(1.0 / n) * np.sum(np.log(proba[np.arange(n), y]))
        return loss

    def _compute_gradients(self, X: np.ndarray, y: np.ndarray):
        """
        计算多分类交叉熵对 W 和 b 的梯度。

        ∂J/∂W = (1/n) X^T (ŷ - y_onehot)
        ∂J/∂b = (1/n) Σ (ŷ_i - y_onehot_i)

        参数:
            X: np.ndarray, (n_samples, n_features)
            y: np.ndarray, (n_samples,)

        返回:
            dW: np.ndarray, (n_features, n_classes)
            db: np.ndarray, (n_classes,)
        """
        n = len(y)
        proba = self._predict_proba(X)  # (n, K)

        # 构建 one-hot 编码的真实标签
        y_onehot = np.zeros((n, self.n_classes))  # (n, K) 全零矩阵
        y_onehot[np.arange(n), y] = 1  # 在真实类别位置置 1

        errors = proba - y_onehot  # 预测误差矩阵 (n, K)

        dW = (1.0 / n) * (X.T @ errors)  # (d, K)
        db = (1.0 / n) * np.sum(errors, axis=0)  # (K,)

        return dW, db

    def fit(self, X: np.ndarray, y: np.ndarray, verbose: bool = True):
        """
        训练 Softmax 回归模型。

        参数:
            X: np.ndarray, (n_samples, n_features)
            y: np.ndarray, (n_samples,)，类别标签为 0, 1, ..., K-1
            verbose: bool, 是否打印日志
        """
        n_samples, n_features = X.shape
        self.n_classes = len(np.unique(y))  # 类别数

        # 初始化参数
        self.W = np.random.randn(n_features, self.n_classes) * 0.01  # (d, K)
        self.b = np.zeros(self.n_classes)  # (K,)
        self.loss_history = []

        for epoch in range(self.max_epochs):
            loss = self._compute_loss(X, y)
            self.loss_history.append(loss)

            dW, db = self._compute_gradients(X, y)

            self.W -= self.learning_rate * dW  # 矩阵更新
            self.b -= self.learning_rate * db  # 向量更新

            if len(self.loss_history) > 1:
                if abs(self.loss_history[-2] - loss) < self.tolerance:
                    if verbose:
                        print(f"第 {epoch + 1} 轮收敛！")
                    break

            if verbose and (epoch + 1) % 200 == 0:
                print(f"Epoch {epoch + 1:4d}: loss={loss:.6f}")

        if verbose:
            print(f"\nSoftmax 训练完成，共 {epoch + 1} 轮，最终损失: {self.loss_history[-1]:.6f}")


# ============================================================================
# 第四部分：可视化
# ============================================================================

def plot_sigmoid_curve():
    """绘制 Sigmoid 函数曲线，帮助理解其形态。"""
    z = np.linspace(-8, 8, 500)  # z 从 -8 到 8
    s = sigmoid(z)  # 对应的 σ(z) 值

    fig, ax = plt.subplots(figsize=(8, 5))
    ax.plot(z, s, 'b-', linewidth=2.5, label=r'$\sigma(z) = 1/(1+e^{-z})$')

    # 关键点标记
    ax.scatter([0], [0.5], c='red', s=100, zorder=5)  # z=0 时 σ=0.5
    ax.axhline(y=0.5, color='red', linestyle='--', alpha=0.5, linewidth=1)  # 0.5 水平线
    ax.axvline(x=0, color='red', linestyle='--', alpha=0.5, linewidth=1)  # z=0 竖直线
    ax.axhline(y=0, color='gray', linestyle=':', linewidth=0.8)  # 渐近线 y=0
    ax.axhline(y=1, color='gray', linestyle=':', linewidth=0.8)  # 渐近线 y=1

    # 区域标注
    ax.annotate('Positive Region\n(z > 0 -> sigma > 0.5)', xy=(3, 0.95), fontsize=11,
                ha='center', bbox=dict(boxstyle='round', facecolor='lightgreen', alpha=0.7))
    ax.annotate('Negative Region\n(z < 0 -> sigma < 0.5)', xy=(-3, 0.05), fontsize=11,
                ha='center', bbox=dict(boxstyle='round', facecolor='lightcoral', alpha=0.7))
    ax.annotate('Decision Boundary\n(z = 0, sigma = 0.5)', xy=(0, 0.5), fontsize=10,
                xytext=(1.5, 0.35), arrowprops=dict(arrowstyle='->', color='red'),
                ha='center', color='red')

    ax.set_xlabel('z = w*x + b', fontsize=13)
    ax.set_ylabel('sigma(z)', fontsize=13)
    ax.set_title('Sigmoid Function: Mapping Reals to [0, 1]', fontsize=14)
    ax.legend(fontsize=11)
    ax.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.savefig(os.path.join(_IMAGES_DIR, 'sigmoid_curve.png'), dpi=150, bbox_inches='tight')
    plt.show()
    print(f"Sigmoid 曲线已保存为 {os.path.join(_IMAGES_DIR, 'sigmoid_curve.png')}")


def plot_decision_boundary(model, X, y, title='Logistic Regression Decision Boundary'):
    """
    绘制二维特征空间中的决策边界和概率热力图。

    对于每个 (x1, x2) 坐标点，计算模型输出的概率 P(y=1|x)，
    然后用颜色表示概率大小，形成热力图。

    参数:
        model: 训练好的 LogisticRegression 模型
        X: np.ndarray, 形状 (n, 2) 的特征矩阵（仅用前两个特征）
        y: np.ndarray, 真实标签
        title: str, 图表标题
    """
    # 创建网格
    x_min, x_max = X[:, 0].min() - 0.5, X[:, 0].max() + 0.5
    y_min, y_max = X[:, 1].min() - 0.5, X[:, 1].max() + 0.5
    xx, yy = np.meshgrid(np.linspace(x_min, x_max, 300),
                         np.linspace(y_min, y_max, 300))

    # 计算网格上每个点的预测概率
    grid_points = np.c_[xx.ravel(), yy.ravel()]  # 展开为 (N, 2)
    Z = model._predict_proba(grid_points)  # 概率值
    Z = Z.reshape(xx.shape)  # 恢复为网格形状

    fig, ax = plt.subplots(figsize=(9, 7))

    # 绘制概率热力图（蓝色=低概率，红色=高概率）
    contour = ax.contourf(xx, yy, Z, levels=20, cmap='RdBu', alpha=0.6)

    # 绘制决策边界线（σ = 0.5 的等高线）
    ax.contour(xx, yy, Z, levels=[0.5], colors='green', linewidths=2.5,
               linestyles='-')

    # 绘制数据点
    ax.scatter(X[y == 0, 0], X[y == 0, 1], c='blue', marker='o',
               edgecolors='k', s=60, label='Negative Class (y=0)', alpha=0.8)
    ax.scatter(X[y == 1, 0], X[y == 1, 1], c='red', marker='^',
               edgecolors='k', s=60, label='Positive Class (y=1)', alpha=0.8)

    # 添加颜色条
    cbar = plt.colorbar(contour, ax=ax)
    cbar.set_label('P(y=1|x)', fontsize=11)

    ax.set_xlabel('Feature x1', fontsize=13)
    ax.set_ylabel('Feature x2', fontsize=13)
    ax.set_title(title, fontsize=14)
    ax.legend(fontsize=11, loc='upper left')
    ax.grid(True, alpha=0.2)

    plt.tight_layout()
    plt.savefig(os.path.join(_IMAGES_DIR, 'logistic_regression_boundary.png'), dpi=150, bbox_inches='tight')
    plt.show()
    print(f"决策边界图已保存为 {os.path.join(_IMAGES_DIR, 'logistic_regression_boundary.png')}")


def plot_loss_curve(loss_history, title='Training Loss Curve'):
    """绘制训练过程中的损失变化。"""
    fig, ax = plt.subplots(figsize=(8, 5))
    ax.plot(range(1, len(loss_history) + 1), loss_history, 'b-', linewidth=1.5)
    ax.set_xlabel('Epoch', fontsize=12)
    ax.set_ylabel('Cross-Entropy Loss', fontsize=12)
    ax.set_title(title, fontsize=14)
    ax.grid(True, alpha=0.3)
    plt.tight_layout()
    plt.savefig(os.path.join(_IMAGES_DIR, 'loss_curve.png'), dpi=150, bbox_inches='tight')
    plt.show()
    print(f"损失曲线已保存为 {os.path.join(_IMAGES_DIR, 'loss_curve.png')}")


# ============================================================================
# 第五部分：主程序
# ============================================================================

def main():
    """
    主函数：演示逻辑回归的二分类和多分类完整流程。
    """
    print("=" * 60)
    print("逻辑回归从零实现 — s03_logistic_regression")
    print("=" * 60)

    # ---- 1. Sigmoid 函数可视化 ----
    print("\n[步骤 1] 可视化 Sigmoid 函数...")
    plot_sigmoid_curve()

    # ---- 2. 加载 Iris 数据集 ----
    print("\n[步骤 2] 加载 Iris 数据集...")
    iris = load_iris()
    X_full = iris.data  # 全部 4 个特征
    y_full = iris.target  # 全部 3 个类别
    print(f"数据集: {X_full.shape[0]} 样本, {X_full.shape[1]} 特征, "
          f"{len(np.unique(y_full))} 个类别")
    print(f"类别名称: {iris.target_names}")
    print(f"特征名称: {iris.feature_names}")

    # ---- 3. 二分类 ----
    print("\n" + "=" * 40)
    print("[步骤 3] 二分类逻辑回归（类别 0 vs 类别 1）")
    print("=" * 40)

    # 取前两个类别和前两个特征（便于可视化）
    mask_binary = (y_full == 0) | (y_full == 1)  # 只取类别 0 和 1
    X_binary = X_full[mask_binary][:, :2]  # 取前两个特征用于二维可视化
    y_binary = y_full[mask_binary]
    print(f"二分类数据: {X_binary.shape[0]} 样本, 类别 0: {np.sum(y_binary == 0)}, "
          f"类别 1: {np.sum(y_binary == 1)}")

    # 划分训练/测试集（80/20）
    X_tr, X_te, y_tr, y_te = train_test_split(
        X_binary, y_binary, test_size=0.2, random_state=42
    )

    # 训练逻辑回归模型
    model_binary = LogisticRegression(learning_rate=0.5, max_epochs=2000)
    model_binary.fit(X_tr, y_tr)

    # 评估
    y_pred = model_binary.predict(X_te)
    accuracy = np.mean(y_pred == y_te)
    print(f"\n测试集准确率: {accuracy:.2%}")

    # 混淆矩阵
    tp = np.sum((y_pred == 1) & (y_te == 1))  # 真正例
    tn = np.sum((y_pred == 0) & (y_te == 0))  # 真负例
    fp = np.sum((y_pred == 1) & (y_te == 0))  # 假正例
    fn = np.sum((y_pred == 0) & (y_te == 1))  # 假负例
    print(f"混淆矩阵: TP={tp}, TN={tn}, FP={fp}, FN={fn}")

    # 可视化
    plot_decision_boundary(model_binary, X_tr, y_tr,
                           title='Binary Logistic Regression - Decision Boundary & Probability Heatmap')
    plot_loss_curve(model_binary.loss_history, title='Binary Logistic Regression Training Loss')

    # ---- 4. 多分类（Softmax 回归） ----
    print("\n" + "=" * 40)
    print("[步骤 4] 多分类 Softmax 回归（全部 3 个类别）")
    print("=" * 40)

    # 取前两个特征用于训练和可视化
    X_multi = X_full[:, :2]
    y_multi = y_full

    X_tr_m, X_te_m, y_tr_m, y_te_m = train_test_split(
        X_multi, y_multi, test_size=0.2, random_state=42
    )

    # 训练 Softmax 回归
    model_softmax = SoftmaxRegression(learning_rate=0.5, max_epochs=5000)
    model_softmax.fit(X_tr_m, y_tr_m)

    # 评估
    y_pred_m = model_softmax.predict(X_te_m)
    accuracy_m = np.mean(y_pred_m == y_te_m)
    print(f"\n测试集准确率: {accuracy_m:.2%}")

    # 可视化多分类决策边界
    x_min, x_max = X_tr_m[:, 0].min() - 0.5, X_tr_m[:, 0].max() + 0.5
    y_min, y_max = X_tr_m[:, 1].min() - 0.5, X_tr_m[:, 1].max() + 0.5
    xx, yy = np.meshgrid(np.linspace(x_min, x_max, 300),
                         np.linspace(y_min, y_max, 300))
    grid_points = np.c_[xx.ravel(), yy.ravel()]
    Z_m = model_softmax.predict(grid_points).reshape(xx.shape)

    fig, ax = plt.subplots(figsize=(9, 7))
    ax.contourf(xx, yy, Z_m, levels=np.arange(-0.5, 3.5, 1),
                colors=['#E8F5E9', '#FFF3E0', '#E3F2FD'], alpha=0.6)
    scatter = ax.scatter(X_tr_m[:, 0], X_tr_m[:, 1], c=y_tr_m,
                         cmap='viridis', edgecolors='k', s=60, alpha=0.8)
    legend = ax.legend(*scatter.legend_elements(), title='Class',
                        fontsize=10, title_fontsize=11)
    ax.set_xlabel('Feature x1 (Sepal Length)', fontsize=13)
    ax.set_ylabel('Feature x2 (Sepal Width)', fontsize=13)
    ax.set_title('Softmax Multi-Class - Decision Regions', fontsize=14)
    ax.grid(True, alpha=0.2)
    plt.tight_layout()
    plt.savefig(os.path.join(_IMAGES_DIR, 'softmax_multiclass_boundary.png'), dpi=150, bbox_inches='tight')
    plt.show()
    print(f"多分类决策区域图已保存为 {os.path.join(_IMAGES_DIR, 'softmax_multiclass_boundary.png')}")

    print("\n" + "=" * 60)
    print("演示完成！")
    print("=" * 60)


if __name__ == '__main__':
    main()

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628

s03 逻辑回归 — demo.py 代码详解 ​

运行方式 ​

代码逐段详解 ​

第1步：导入库 — 每个库是做什么的 ​

第2步：Sigmoid 函数 — 从实数到概率 ​

第3步：Softmax 函数 — 从得分到概率分布 ​

第4步：二分类逻辑回归模型 ​

4.1 模型定义 ​

4.2 损失函数：二元交叉熵 ​

4.3 梯度：数学之美的巅峰 ​

4.4 训练循环 ​

第5步：多分类 Softmax 回归 ​

第6步：可视化 ​

6.1 Sigmoid 函数曲线 ​

6.2 决策边界与概率热力图 ​

第7步：模型评估 — 混淆矩阵 ​

关键概念速查表 ​

完整代码 ​

s03 逻辑回归 — demo.py 代码详解

运行方式

代码逐段详解

第1步：导入库 — 每个库是做什么的

第2步：Sigmoid 函数 — 从实数到概率

第3步：Softmax 函数 — 从得分到概率分布

第4步：二分类逻辑回归模型

4.1 模型定义

4.2 损失函数：二元交叉熵

4.3 梯度：数学之美的巅峰

4.4 训练循环

第5步：多分类 Softmax 回归

第6步：可视化

6.1 Sigmoid 函数曲线

6.2 决策边界与概率热力图

第7步：模型评估 — 混淆矩阵

关键概念速查表

完整代码