s02 线性回归 — demo.py 代码详解

Download demo.py

运行方式

cd s02_linear_regression/code
python demo.py

代码逐段详解

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

import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from mpl_toolkits.mplot3d import Axes3D
from sklearn.linear_model import LinearRegression as SklearnLR

• os：操作文件路径，用于创建 images/ 目录
• numpy：数值计算核心。关键用法：
  • np.random.uniform() 生成均匀分布的随机数
  • np.mean() 计算平均值（MSE 损失的核心）
  • np.sum() 求和（梯度计算中累加误差）
  • np.linalg.inv() 矩阵求逆（正规方程求解）
  • np.column_stack() 按列拼接（为正规方程构建增广矩阵）
  • np.meshgrid() 创建网格（绘制损失函数等高线）
• matplotlib：绘图，包括 2D 散点/线图、3D 曲面、等高线图
• Axes3D（来自 mpl_toolkits.mplot3d）：支持 3D 绘图，用于展示损失函数的 3D 曲面
• sklearn.linear_model.LinearRegression：scikit-learn 的标准线性回归实现，用作基准对比，验证我们从头实现的正确性

第2步：数据生成

def generate_regression_data(n_samples=100, noise_std=3.0, true_w=2.0, true_b=5.0, random_seed=42):

生成回归任务的合成数据。真实函数为：

$$y=2x+5+ϵ,ϵ\sim \mathcal{N}(0,3^2)$$

具体步骤：

1. np.random.uniform(0, 10, n_samples) 在 $[0,10]$ 区间均匀采样 100 个 $x$ 值
2. np.random.randn(n_samples) * noise_std 生成标准差为 3.0 的高斯噪声
3. 按 $y=2x+5+\text{noise}$ 计算目标值

噪声标准差 3.0 相对较大（真实值范围约 $[5,25]$），这使得数据点较为分散，更接近真实场景中的数据。

第3步：线性回归模型 — 梯度下降法

3.1 模型假设

线性回归假设输出是输入的线性函数：

$$\hat{y}=wx+b$$

其中 $w$ 是权重（斜率），$b$ 是偏置（截距）。这两个参数就是模型需要从数据中"学习"的内容。

3.2 损失函数：均方误差 MSE

$$J(w,b)=\frac{1}{n}\sum _{i=1}^n({\hat{y}}_i-y_i{)}^2$$

为什么用平方而不是绝对值？ 四个关键原因：

1. 处处可导：平方函数光滑，绝对值在 $x=0$ 处不可导——梯度下降需要处处有梯度
2. 对大误差更敏感：误差为 10 的样本，MSE 中惩罚 100，MAE 中仅惩罚 10——这鼓励模型优先修正大偏差
3. 概率解释：若假设误差服从正态分布 $ϵ\sim \mathcal{N}(0,{\sigma }^2)$，最小化 MSE 等价于最大似然估计（MLE）
4. 凸函数：MSE 关于 $(w,b)$ 是凸函数，只有一个全局最小值，不会被卡在局部最优

代码实现：

def _compute_loss(self, X, y):
    y_pred = self.predict(X)            # 计算所有预测值
    return np.mean((y_pred - y) ** 2)   # 平均平方误差

np.mean((y_pred - y) ** 2) 是向量化操作：先逐元素计算平方差，再取平均。一条语句完成全部 $n$ 个样本的 MSE 计算。

3.3 梯度推导与计算

MSE 损失对 $w$ 和 $b$ 的偏导数（通过链式法则推导）：

$$\frac{\partial J}{\partial w}=\frac{2}{n}\sum _{i=1}^n({\hat{y}}_i-y_i)\cdot x_i$$$$\frac{\partial J}{\partial b}=\frac{2}{n}\sum _{i=1}^n({\hat{y}}_i-y_i)$$

代码实现：

def _compute_gradients(self, X, y):
    y_pred = self.predict(X)          # 计算预测值 (n,)
    n = len(y)
    errors = y_pred - y               # 误差向量 (n,)
    dw = (2.0 / n) * np.sum(errors * X)  # ∂J/∂w: 逐元素乘后求和
    db = (2.0 / n) * np.sum(errors)      # ∂J/∂b: 误差求和
    return dw, db

errors * X 是逐元素乘法：对每个样本 $i$，计算 $({\hat{y}}_i-y_i)\cdot x_i$，然后用 np.sum() 累加。系数 $2/n$ 中，因子 2 来自平方函数的导数，$1/n$ 来自平均。

3.4 梯度下降更新规则

$$w\leftarrow w-\eta \cdot \frac{\partial J}{\partial w}$$$$b\leftarrow b-\eta \cdot \frac{\partial J}{\partial b}$$

代码：

self.w -= self.learning_rate * dw
self.b -= self.learning_rate * db

注意这里是 -=（减等于），因为我们沿梯度的反方向走——梯度指向函数值上升最快的方向，我们要下降所以要反向。

收敛条件：当相邻两轮的损失变化小于 tolerance（默认 ${10}^{-6}$）时，认为已收敛，提前终止训练。

参数初始化：

self.w = np.random.randn() * 0.1
self.b = np.random.randn() * 0.1

用标准正态分布的小随机数初始化。为什么不是全零？对于单变量线性回归，全零初始化也可以工作（因为只有一个 $w$，不存在对称性问题），但用随机初始化是更好的通用实践。乘 0.1 是为了让初始值接近 0 但不完全为 0。

第4步：正规方程 — 线性回归的解析解

对于线性回归，我们不仅能梯度下降，还能直接求出闭式解。这是少数能写出解析解的机器学习模型之一。

推导：将 MSE 写成矩阵形式 $J(\theta )=\frac{1}{n}\parallel X\theta -y{\parallel }^2$，令 ${\mathrm{\nabla }}_{\theta }J=0$，解得：

$${\theta }^{\ast }=(X^{T}X{)}^{-1}{X}^{T}y$$

这就是正规方程（Normal Equation）。

代码实现：

def normal_equation_solution(X, y):
    n = len(X)
    # 构建增广矩阵: X_aug = [x, 1]，第一列是x，第二列是全1（对应偏置）
    X_aug = np.column_stack([X, np.ones(n)])
    # θ = (X^T X)^{-1} X^T y
    theta = np.linalg.inv(X_aug.T @ X_aug) @ X_aug.T @ y
    w = theta[0]  # 权重
    b = theta[1]  # 偏置
    return w, b

关键细节：

• np.column_stack([X, np.ones(n)]) 构建增广矩阵 $[x,\mathbf{1}]$，形状 $(n,2)$。第二列全 1 对应偏置项——这样 $\theta =[w,b{]}^T$ 就能用一个矩阵方程求解。
• @ 是 Python 3.5+ 的矩阵乘法运算符，等价于 np.matmul()。X_aug.T @ X_aug 计算 $X^{T}X$，形状 $(2,2)$。
• np.linalg.inv() 计算矩阵的逆。对于 $2\times 2$ 矩阵，求逆非常快。

正规方程 vs 梯度下降：

• 正规方程：一步到位得到精确解，无需选择学习率，无需迭代。但 $X^{T}X$ 求逆的复杂度是 $O(d^3)$，当特征维度 $d$ 很大时不可行。
• 梯度下降：需要选择学习率，需要多轮迭代，但复杂度是 $O(nd)$ 每轮，适合大规模数据和深度学习。

第5步：可视化 — 四合一分析图

代码生成了一个 $2\times 2$ 的子图布局：

子图 1：数据散点和拟合直线

三条直线分别来自梯度下降法（红色实线）、正规方程（绿色虚线）和 sklearn（蓝色点划线）。理想情况下三条线几乎重合——三种方法给出的 $w,b$ 非常接近，互相验证了实现的正确性。

子图 2：训练损失曲线

横轴是 epoch，纵轴是 MSE 损失（对数刻度）。损失通常在前几十轮快速下降，之后趋于平稳。对数刻度让早期快速下降和后期精细调整都能看清楚。

子图 3：损失函数等高线 + 梯度下降轨迹

这是在 $(w,b)$ 参数空间中绘制 MSE 的等高线图。同心椭圆（或近似椭圆）是等损失线——同一椭圆上的 $(w,b)$ 组合给出相同的 MSE。

红色轨迹是梯度下降的优化路径：从蓝色起点出发，沿着局部梯度方向一步步走向红色星形的最优点。轨迹垂直于等高线（因为梯度方向垂直于等高线），且步长越来越小（接近最优点时梯度趋近于零）。

生成等高线用的 np.meshgrid() 在 $(w,b)$ 平面上创建了 $100\times 100$ 的网格，对每个网格点计算 MSE。损失函数的等高线呈椭圆形状，说明 $w$ 和 $b$ 的最优值之间存在一定的耦合关系。

子图 4：三种方法的参数对比

柱状图直观比较梯度下降、正规方程和 sklearn 三种方法得到 $w$ 和 $b$ 值。灰色虚线标注了真实值 $w=2.0$ 和 $b=5.0$。

第6步：学习率对比实验

compare_learning_rates() 函数用三种不同的学习率（0.001、0.01、0.05）训练模型，在一张图上对比它们的损失曲线：

• $\eta =0.001$：收敛最慢，200 轮后损失仍较高
• $\eta =0.01$：适中，在大约 300 轮时收敛到较小的损失值
• $\eta =0.05$：最快收敛，但如果太大可能震荡或发散

这个实验直观展示了学习率作为"步长"的含义：步长太小，下山太慢；步长适中，高效到达谷底；步长太大，可能在谷底来回跳跃甚至发散。

第7步：模型评估 — MSE 和 $R^2$

在 main() 中，通过计算 MSE 和 $R^2$ 来评估模型：

• MSE（均方误差）：$\text{MSE}=\frac{1}{n}\sum ({\hat{y}}_i-y_i{)}^2$，单位是目标值的平方。越小越好，但绝对值依赖于数据的尺度。
• $R^2$（决定系数）：$R^2=1-\frac{S{S}_{\text{res}}}{S{S}_{\text{tot}}}$，衡量模型解释了数据中多少比例的方差。越接近 1 表示模型解释力越强。$S{S}_{\text{res}}=\sum (y_i-{\hat{y}}_i{)}^2$ 是残差平方和，$S{S}_{\text{tot}}=\sum (y_i-\overline{y}{)}^2$ 是总平方和（$\overline{y}$ 是 $y$ 的均值）。

$R^2$ 可以被理解为"相比于只用均值来预测，模型减少了多少比例的误差"。$R^2=0$ 表示模型不比均值好，$R^2=1$ 表示完美预测。

关键概念速查表

概念	数学形式	代码位置	关键说明
线性模型	$\hat{y}=wx+b$	predict()	最简单的参数化模型
MSE 损失	$J=\frac{1}{n}\sum ({\hat{y}}_i-y_i{)}^2$	_compute_loss()	处处可导，对大误差惩罚重
梯度 $\partial J/\partial w$	$\frac{2}{n}\sum ({\hat{y}}_i-y_i)x_i$	_compute_gradients()	链式法则推导
梯度 $\partial J/\partial b$	$\frac{2}{n}\sum ({\hat{y}}_i-y_i)$	_compute_gradients()	比 $w$ 的梯度少乘 $x_i$
梯度下降更新	$\theta \leftarrow \theta -\eta {\mathrm{\nabla }}_{\theta }J$	fit()	self.w -= lr * dw
正规方程	${\theta }^{\ast }=(X^{T}X{)}^{-1}{X}^{T}y$	normal_equation_solution()	闭式解，$O(d^3)$
学习率 $\eta$	步长	__init__()	太小慢，太大震荡
$R^2$	$1-\frac{S{S}_{\text{res}}}{S{S}_{\text{tot}}}$	main()	模型解释方差比例

完整代码

# -*- coding: utf-8 -*-
"""
===============================================================================
s02_linear_regression/code/demo.py — 线性回归从零实现
===============================================================================
本演示从零实现线性回归，涵盖梯度下降法和正规方程两种求解方式，
并可视化数据、拟合直线、损失曲线等内容。

通过本演示，你将理解：
  1. 线性模型 ŷ = wx + b 的数学形式和参数含义
  2. MSE 损失函数及其梯度的推导和计算
  3. 梯度下降法：如何沿着梯度方向一步步逼近最优解
  4. 正规方程：线性回归的封闭形式解析解
  5. 学习率对收敛速度的影响
  6. 与 sklearn 标准实现的对比验证

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

import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from mpl_toolkits.mplot3d import Axes3D  # 用于绘制 3D 损失曲面
from sklearn.linear_model import LinearRegression 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 generate_regression_data(n_samples: int = 100, noise_std: float = 3.0,
                             true_w: float = 2.0, true_b: float = 5.0,
                             random_seed: int = 42):
    """
    生成线性回归的合成数据。

    按照 y = true_w * x + true_b + noise 生成数据，
    其中 noise 是服从 N(0, noise_std²) 的高斯噪声。

    参数:
        n_samples: int, 样本数量
        noise_std: float, 噪声的标准差（越大数据越散）
        true_w: float, 真实的斜率
        true_b: float, 真实的截距
        random_seed: int, 随机种子

    返回:
        X: np.ndarray, 形状 (n_samples,) 的特征向量
        y: np.ndarray, 形状 (n_samples,) 的目标值向量
    """
    np.random.seed(random_seed)  # 固定随机种子，保证可复现

    # 在 [0, 10] 范围内均匀采样特征 x
    X = np.random.uniform(low=0.0, high=10.0, size=n_samples)

    # 生成噪声：从 N(0, noise_std²) 中采样
    noise = np.random.randn(n_samples) * noise_std

    # 生成目标值：y = 2x + 5 + noise
    y = true_w * X + true_b + noise

    return X, y


# ============================================================================
# 第二部分：线性回归模型（梯度下降法）
# ============================================================================

class LinearRegressionGD:
    """
    使用梯度下降法求解的线性回归模型。

    模型假设: ŷ = w * x + b
    损失函数: J(w,b) = (1/n) * Σ (ŷ_i - y_i)²  (MSE)
    优化方法: 批量梯度下降

    梯度推导:
        ∂J/∂w = (2/n) * Σ (ŷ_i - y_i) * x_i
        ∂J/∂b = (2/n) * Σ (ŷ_i - y_i)

    属性:
        w: float, 权重（斜率）
        b: float, 偏置（截距）
        loss_history: list, 每个 epoch 的损失值（用于绘制损失曲线）
        params_history: list, 每个 epoch 的 (w, b) 值（用于绘制优化轨迹）
    """

    def __init__(self, learning_rate: float = 0.01, 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  # 权重（斜率），将在 fit 中初始化
        self.b = None  # 偏置（截距），将在 fit 中初始化
        self.loss_history = []  # 记录每轮的损失值
        self.params_history = []  # 记录每轮的参数值

    def predict(self, X: np.ndarray) -> np.ndarray:
        """
        对输入 X 进行预测。

        预测公式: ŷ = w * X + b

        参数:
            X: np.ndarray, 形状 (n_samples,) 的特征

        返回:
            np.ndarray, 形状 (n_samples,) 的预测值
        """
        return self.w * X + self.b  # 向量化运算，直接对所有样本预测

    def _compute_loss(self, X: np.ndarray, y: np.ndarray) -> float:
        """
        计算 MSE 损失。

        J(w,b) = (1/n) * Σ (ŷ_i - y_i)²

        参数:
            X: np.ndarray, 特征
            y: np.ndarray, 真实目标值

        返回:
            float, 当前参数下的 MSE 损失值
        """
        y_pred = self.predict(X)  # 计算预测值
        n = len(y)  # 样本数
        return np.mean((y_pred - y) ** 2)  # MSE = 平均平方误差

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

        ∂J/∂w = (2/n) * Σ (ŷ_i - y_i) * x_i
        ∂J/∂b = (2/n) * Σ (ŷ_i - y_i)

        参数:
            X: np.ndarray, 特征
            y: np.ndarray, 真实目标值

        返回:
            dw: float, 损失对 w 的偏导数
            db: float, 损失对 b 的偏导数
        """
        y_pred = self.predict(X)  # 预测值
        n = len(y)  # 样本数
        errors = y_pred - y  # 预测误差向量

        # ∂J/∂w = (2/n) * Σ (ŷ_i - y_i) * x_i
        dw = (2.0 / n) * np.sum(errors * X)

        # ∂J/∂b = (2/n) * Σ (ŷ_i - y_i)
        db = (2.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. 计算当前预测 ŷ = w*x + b
                b. 计算损失 J = MSE(ŷ, y)
                c. 计算梯度 ∂J/∂w 和 ∂J/∂b
                d. 更新参数: w = w - η * ∂J/∂w, b = b - η * ∂J/∂b

        参数:
            X: np.ndarray, 特征
            y: np.ndarray, 真实目标值
            verbose: bool, 是否打印训练日志
        """
        # 用正态分布小随机数初始化 w 和 b，使初始参数接近 0 但不全为 0
        self.w = np.random.randn() * 0.1  # 权重初始化为小随机数
        self.b = np.random.randn() * 0.1  # 偏置初始化为小随机数
        self.loss_history = []  # 清空损失历史
        self.params_history = []  # 清空参数历史

        for epoch in range(self.max_epochs):
            # 步骤 1: 计算当前损失
            loss = self._compute_loss(X, y)
            self.loss_history.append(loss)  # 记录损失
            self.params_history.append((self.w, self.b))  # 记录当前参数

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

            # 步骤 3: 使用梯度下降更新参数
            self.w -= self.learning_rate * dw  # w ← w - η * ∂J/∂w
            self.b -= self.learning_rate * db  # b ← b - η * ∂J/∂b

            # 步骤 4: 检查收敛条件
            if len(self.loss_history) > 1:
                loss_change = abs(self.loss_history[-2] - loss)  # 当前损失与上一轮的差值
                if loss_change < self.tolerance:  # 损失变化很小，认为已收敛
                    if verbose:
                        print(f"第 {epoch + 1} 轮收敛！损失变化 {loss_change:.8f} < {self.tolerance}")
                    break

            # 每 50 轮打印一次训练进度
            if verbose and (epoch + 1) % 50 == 0:
                print(f"Epoch {epoch + 1:4d}: loss={loss:.6f}, w={self.w:.4f}, b={self.b:.4f}")

        if verbose:
            print(f"\n训练完成，共 {epoch + 1} 轮")
            print(f"学到的参数: w = {self.w:.4f}, b = {self.b:.4f}")
            print(f"最终损失: {self.loss_history[-1]:.6f}")


# ============================================================================
# 第三部分：正规方程求解
# ============================================================================

def normal_equation_solution(X: np.ndarray, y: np.ndarray):
    """
    使用正规方程直接求解线性回归的解析解。

    正规方程: θ* = (X^T X)^(-1) X^T y

    其中 X 是包含偏置列的特征矩阵（X 被添加了一列全 1），
    θ 是包含 [w, b] 的参数向量。

    解析解的推导过程:
        J(θ) = (1/n) * ||Xθ - y||²
        令 ∂J/∂θ = 0
        => X^T X θ = X^T y
        => θ = (X^T X)^(-1) X^T y

    参数:
        X: np.ndarray, 形状 (n_samples,) 的特征向量
        y: np.ndarray, 形状 (n_samples,) 的目标值向量

    返回:
        w: float, 权重（斜率）
        b: float, 偏置（截距）
    """
    n = len(X)
    # 构建增广特征矩阵 X_aug = [x, 1]，形状 (n, 2)
    # 第一列是特征 x，第二列是全 1（用于计算偏置）
    X_aug = np.column_stack([X, np.ones(n)])

    # 使用正规方程求解: θ = (X^T X)^(-1) X^T y
    # @ 是 Python 3.5+ 的矩阵乘法运算符，等价于 np.matmul()
    theta = np.linalg.inv(X_aug.T @ X_aug) @ X_aug.T @ y

    w = theta[0]  # 权重（斜率）
    b = theta[1]  # 偏置（截距）
    return w, b


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

def plot_results(X, y, model_gd, w_ne, b_ne, w_sk, b_sk):
    """
    全面可视化线性回归的结果。

    生成一个包含 2 行 2 列的复合图，展示:
      (1) 数据散点与拟合直线
      (2) 训练损失曲线
      (3) 3D 损失曲面上的梯度下降轨迹
      (4) 不同方法的参数对比

    参数:
        X: np.ndarray, 特征
        y: np.ndarray, 目标值
        model_gd: LinearRegressionGD, 梯度下降训练好的模型
        w_ne, b_ne: float, 正规方程解出的参数
        w_sk, b_sk: float, sklearn 解出的参数
    """
    fig = plt.figure(figsize=(16, 12))

    # ---- 子图 1: 数据散点和拟合直线 ----
    ax1 = fig.add_subplot(2, 2, 1)
    ax1.scatter(X, y, c='steelblue', alpha=0.7, s=40, label='Training Data',
                edgecolors='white', linewidth=0.5)

    # 生成一条平滑的 x 序列用于画拟合直线
    X_line = np.linspace(X.min(), X.max(), 200)
    # 梯度下降法拟合的直线（红色）
    ax1.plot(X_line, model_gd.predict(X_line), 'r-', linewidth=2,
             label=f'Gradient Descent: y = {model_gd.w:.2f}x + {model_gd.b:.2f}')
    # 正规方程拟合的直线（绿色虚线）
    ax1.plot(X_line, w_ne * X_line + b_ne, 'g--', linewidth=2,
             label=f'Normal Equation: y = {w_ne:.2f}x + {b_ne:.2f}')
    # sklearn 拟合的直线（蓝色点划线）
    ax1.plot(X_line, w_sk * X_line + b_sk, 'b-.', linewidth=1.5, alpha=0.6,
             label=f'sklearn: y = {w_sk:.2f}x + {b_sk:.2f}')

    ax1.set_xlabel('Feature x', fontsize=12)
    ax1.set_ylabel('Target y', fontsize=12)
    ax1.set_title('Linear Regression: Data and Fitted Lines', fontsize=14)
    ax1.legend(fontsize=9, loc='upper left')
    ax1.grid(True, alpha=0.3)

    # ---- 子图 2: 训练损失曲线 ----
    ax2 = fig.add_subplot(2, 2, 2)
    epochs = range(1, len(model_gd.loss_history) + 1)
    ax2.plot(epochs, model_gd.loss_history, 'b-', linewidth=1.5)
    ax2.set_xlabel('Epoch', fontsize=12)
    ax2.set_ylabel('MSE Loss', fontsize=12)
    ax2.set_title('Loss During Gradient Descent Training', fontsize=14)
    ax2.grid(True, alpha=0.3)
    # 用对数刻度显示 y 轴（因为损失在早期下降很快，后期趋于平稳）
    ax2.set_yscale('log')
    ax2.annotate(f'Initial Loss: {model_gd.loss_history[0]:.2f}',
                 xy=(1, model_gd.loss_history[0]),
                 fontsize=9, color='red')
    ax2.annotate(f'Final Loss: {model_gd.loss_history[-1]:.3f}',
                 xy=(len(epochs), model_gd.loss_history[-1]),
                 fontsize=9, color='green')

    # ---- 子图 3: 损失函数等高线图与优化轨迹 ----
    ax3 = fig.add_subplot(2, 2, 3)

    # 在 (w, b) 平面上计算损失函数的网格值
    w_range = np.linspace(model_gd.w - 1.5, model_gd.w + 1.5, 100)
    b_range = np.linspace(model_gd.b - 3.0, model_gd.b + 3.0, 100)
    W_grid, B_grid = np.meshgrid(w_range, b_range)  # 创建网格

    # 对每个 (w, b) 网格点计算损失
    Z_grid = np.zeros_like(W_grid)
    for i in range(len(b_range)):
        for j in range(len(w_range)):
            w_val = W_grid[i, j]
            b_val = B_grid[i, j]
            y_pred = w_val * X + b_val  # 用当前参数预测
            Z_grid[i, j] = np.mean((y_pred - y) ** 2)  # 计算 MSE

    # 绘制损失函数等高线
    contour = ax3.contour(W_grid, B_grid, Z_grid, levels=20, cmap='viridis', alpha=0.7)
    ax3.clabel(contour, inline=True, fontsize=7)  # 在等高线上标注数值

    # 绘制梯度下降的优化轨迹
    params_arr = np.array(model_gd.params_history)
    ax3.plot(params_arr[:, 0], params_arr[:, 1], 'r.-', markersize=2, linewidth=1,
             label='GD Trajectory')

    # 标记起点和终点
    ax3.scatter(params_arr[0, 0], params_arr[0, 1], c='blue', s=100, marker='o',
                zorder=5, label=f'Start (w={params_arr[0,0]:.2f}, b={params_arr[0,1]:.2f})')
    ax3.scatter(params_arr[-1, 0], params_arr[-1, 1], c='red', s=100, marker='*',
                zorder=5, label=f'End (w={params_arr[-1,0]:.2f}, b={params_arr[-1,1]:.2f})')

    ax3.set_xlabel('Weight w', fontsize=12)
    ax3.set_ylabel('Bias b', fontsize=12)
    ax3.set_title('Loss Contour and Gradient Descent Trajectory', fontsize=14)
    ax3.legend(fontsize=8, loc='upper right')

    # ---- 子图 4: 方法对比条形图 ----
    ax4 = fig.add_subplot(2, 2, 4)
    methods = ['Gradient Descent', 'Normal Equation', 'sklearn']
    w_values = [model_gd.w, w_ne, w_sk]
    b_values = [model_gd.b, b_ne, b_sk]
    x_pos = np.arange(len(methods))
    width = 0.35

    bars1 = ax4.bar(x_pos - width/2, w_values, width, label='Weight w', color='steelblue', alpha=0.8)
    bars2 = ax4.bar(x_pos + width/2, b_values, width, label='Bias b', color='coral', alpha=0.8)

    # 在每个柱状图上方标注数值
    for bar in bars1:
        ax4.text(bar.get_x() + bar.get_width()/2., bar.get_height() + 0.01,
                 f'{bar.get_height():.3f}', ha='center', va='bottom', fontsize=9)
    for bar in bars2:
        ax4.text(bar.get_x() + bar.get_width()/2., bar.get_height() + 0.01,
                 f'{bar.get_height():.3f}', ha='center', va='bottom', fontsize=9)

    ax4.set_xticks(x_pos)
    ax4.set_xticklabels(methods, fontsize=11)
    ax4.set_ylabel('Parameter Value', fontsize=12)
    ax4.set_title('Parameter Comparison Across Three Methods', fontsize=14)
    ax4.legend(fontsize=10)
    ax4.axhline(y=2.0, color='gray', linestyle='--', alpha=0.5, label='True w=2.0')
    ax4.axhline(y=5.0, color='gray', linestyle=':', alpha=0.5, label='True b=5.0')

    plt.tight_layout()
    plt.savefig(os.path.join(_IMAGES_DIR, 'linear_regression_results.png'), dpi=150, bbox_inches='tight')
    plt.show()
    print(f"\n图片已保存为 {os.path.join(_IMAGES_DIR, 'linear_regression_results.png')}")


def compare_learning_rates(X, y):
    """
    比较不同学习率对梯度下降收敛的影响。

    分别用 3 种不同的学习率训练模型，并在同一张图上对比它们的损失曲线。

    参数:
        X: np.ndarray, 特征
        y: np.ndarray, 目标值
    """
    rates = [0.001, 0.01, 0.05]  # 三种学习率
    colors = ['blue', 'green', 'red']  # 对应的颜色

    fig, ax = plt.subplots(figsize=(10, 5))

    for lr, color in zip(rates, colors):
        model = LinearRegressionGD(learning_rate=lr, max_epochs=200)
        model.fit(X, y, verbose=False)
        epochs = range(1, len(model.loss_history) + 1)
        ax.plot(epochs, model.loss_history, color=color, linewidth=1.5,
                label=f'lr={lr} (loss={model.loss_history[-1]:.2e})')

    ax.set_xlabel('Epoch', fontsize=12)
    ax.set_ylabel('MSE Loss', fontsize=12)
    ax.set_title('Effect of Learning Rate on Convergence Speed', fontsize=14)
    ax.legend(fontsize=8)
    ax.grid(True, alpha=0.3)
    ax.set_yscale('log')
    from matplotlib.ticker import ScalarFormatter
    ax.yaxis.set_major_formatter(ScalarFormatter())
    ax.ticklabel_format(axis='y', style='sci', scilimits=(-2, 3))

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


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

def main():
    """
    主函数：串联整个线性回归的教学演示流程。
    """
    print("=" * 60)
    print("线性回归从零实现 — s02_linear_regression")
    print("=" * 60)

    # 1. 生成合成数据（真实参数: w=2, b=5）
    print("\n[步骤 1] 生成合成数据 (y = 2x + 5 + noise)...")
    X, y = generate_regression_data(n_samples=100, true_w=2.0, true_b=5.0,
                                     noise_std=3.0, random_seed=42)
    print(f"数据形状: X={X.shape}, y={y.shape}")
    print(f"X 范围: [{X.min():.2f}, {X.max():.2f}]")
    print(f"y 范围: [{y.min():.2f}, {y.max():.2f}]")

    # 2. 使用梯度下降法训练
    print("\n[步骤 2] 使用梯度下降法训练线性回归模型...")
    model_gd = LinearRegressionGD(learning_rate=0.01, max_epochs=500)
    model_gd.fit(X, y, verbose=True)

    # 3. 使用正规方程求解
    print("\n[步骤 3] 使用正规方程求解...")
    w_ne, b_ne = normal_equation_solution(X, y)
    print(f"正规方程解: w = {w_ne:.4f}, b = {b_ne:.4f}")

    # 4. 使用 sklearn 求解（对比验证）
    print("\n[步骤 4] 使用 sklearn 求解（作为基准）...")
    X_sk = X.reshape(-1, 1)  # sklearn 要求特征为二维: (n_samples, n_features)
    model_sk = SklearnLR()
    model_sk.fit(X_sk, y)
    w_sk = model_sk.coef_[0]  # sklearn 的权重在 coef_ 属性中
    b_sk = model_sk.intercept_  # sklearn 的偏置在 intercept_ 属性中
    print(f"sklearn 解: w = {w_sk:.4f}, b = {b_sk:.4f}")

    # 5. 评估模型
    print("\n[步骤 5] 评估模型...")
    y_pred = model_gd.predict(X)  # 用梯度下降模型预测
    mse = np.mean((y_pred - y) ** 2)  # 计算 MSE
    # R² 分数: 1 - SS_res / SS_tot，越接近 1 表示模型解释力越强
    ss_res = np.sum((y - y_pred) ** 2)  # 残差平方和
    ss_tot = np.sum((y - np.mean(y)) ** 2)  # 总平方和
    r2 = 1 - ss_res / ss_tot
    print(f"MSE = {mse:.4f}")
    print(f"R²  = {r2:.4f}")

    # 6. 综合可视化
    print("\n[步骤 6] 综合可视化...")
    plot_results(X, y, model_gd, w_ne, b_ne, w_sk, b_sk)

    # 7. 不同学习率对比
    print("\n[步骤 7] 比较不同学习率的效果...")
    compare_learning_rates(X, y)

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


if __name__ == '__main__':
    main()