从公式到代码：手把手构建你的第一个线性回归模型

很多朋友在学吴恩达老师的机器学习课程时，都有过类似的困惑：那些优雅的数学公式，比如代价函数、梯度下降，听起来逻辑清晰，可一旦打开编辑器，面对空白的Python文件，却不知从何下手。理论上的“理解”和实际能“跑起来”的代码之间，似乎隔着一道鸿沟。这篇文章，就是为你准备的桥梁。我们不谈空泛的概念，而是聚焦于一个最经典也最基础的算法——线性回归，用Python一步步把它从数学符号变成可运行、可调试的代码。无论你是想巩固理论基础，还是为面试积累实战项目，这个过程都将让你对机器学习如何“工作”有更扎实的把握。

1. 环境搭建与数据初探

在动手写任何算法之前，一个干净、可复现的环境是高效工作的基石。我强烈建议使用Anaconda来管理你的Python环境，它能很好地处理包依赖问题。对于这个项目，我们需要的核心库并不多：

# 创建并激活一个专用于本项目的虚拟环境（可选但推荐）
conda create -n linear_regression_demo python=3.9
conda activate linear_regression_demo

# 安装必要库
pip install numpy pandas matplotlib scikit-learn jupyter

• numpy 是进行矩阵和数值计算的基石，没有它，机器学习中的向量化操作将寸步难行。
• pandas 用于数据加载和初步的清洗、查看，它让处理表格数据变得异常轻松。
• matplotlib 是我们的可视化工具，无论是绘制原始数据散点图，还是观察梯度下降过程，都离不开它。
• scikit-learn 在这里并非用于直接调用它的线性回归模型，而是借用其丰富的数据集和工具函数，例如数据划分。
• jupyter 提供了一个交互式的探索环境，非常适合一步步调试和观察中间结果。

接下来，我们加载数据。为了完全复现课程作业的感觉，我们可以使用一个经典的数据集：波士顿房价数据集（虽然该数据集因伦理问题已从sklearn最新版本中移除，但其历史地位无可替代，我们可用其他类似数据集或自行构造）。这里，为了纯粹聚焦算法，我们自己构造一份模拟数据。这有一个额外好处：你清楚地知道数据是怎么来的，便于理解算法如何拟合。

import numpy as np
import matplotlib.pyplot as plt

# 设置随机种子，确保每次运行结果一致
np.random.seed(42)

# 生成特征X和目标y
# 假设真实关系为 y = 2.5 * X + 1.0 + 噪声
m = 100  # 样本数量
X = 2 * np.random.rand(m, 1)  # 生成0到2之间的均匀分布，形状(100, 1)
y = 2.5 * X + 1.0 + np.random.randn(m, 1)  # 加入高斯噪声

# 快速可视化，看看数据长什么样
plt.figure(figsize=(10, 6))
plt.scatter(X, y, alpha=0.7, edgecolors='k', label='Training Data')
plt.xlabel('Feature X (e.g., House Size)')
plt.ylabel('Target y (e.g., Price)')
plt.title('Simulated Linear Regression Dataset')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.show()

运行这段代码，你应该能看到大致呈线性趋势的散点图。我们的目标，就是找到一条直线，最好地穿过这些点。记住这个图，后面我们会把算法找到的直线画上去做对比。

注意：在实际项目中，拿到数据后的第一步永远是探索性数据分析(EDA)。你需要查看数据的基本统计信息（均值、标准差）、检查缺失值、观察特征分布以及特征与目标之间的关系。对于线性回归，尤其要留意异常值，因为它们会对最小二乘估计产生巨大影响。

2. 线性回归的核心：代价函数与梯度下降

现在进入正题。线性回归的目标是找到参数 $\theta_0$ (截距) 和 $\theta_1$ (斜率)，使得我们的预测值 $\hat{y} = \theta_0 + \theta_1 X$ 与真实值 $y$ 之间的差距最小。这个差距用代价函数（Cost Function） 来衡量，最常用的是均方误差（MSE）。

2.1 实现代价函数

代价函数 $J(\theta_0, \theta_1)$ 的数学定义是：

$$J(\theta_0, \theta_1) = \frac{1}{2m} \sum_{i=1}^{m} (\hat{y}^{(i)} - y^{(i)})^2 = \frac{1}{2m} \sum_{i=1}^{m} (\theta_0 + \theta_1 X^{(i)} - y^{(i)})^2$$

这里除以 $2m$ 而不是 $m$ 是为了后续求导时形式更整洁（系数2会被消掉）。让我们用Python实现它，并采用向量化的写法，这比用for循环高效得多。

def compute_cost(X, y, theta):
    """
    计算线性回归的代价函数（MSE）。
    参数：
    X : 特征矩阵，形状 (m, n)，其中n是特征数（本例中n=1，但函数支持多元）
    y : 目标向量，形状 (m, 1)
    theta : 参数向量，形状 (n+1, 1)
    返回：
    cost : 标量，代价函数值
    """
    m = len(y)  # 样本数量
    # 计算预测值。注意：X需要添加一列1以对应theta_0 (偏置项)
    # 这里假设传入的X是原始特征，未添加偏置列。我们在函数内部处理。
    # 更通用的做法是，在外部将X准备好为 [1, X] 的形式。
    # 为了清晰，我们假设X已经是添加了偏置列的形式，即第一列全为1。
    predictions = X.dot(theta)  # (m, n+1) dot (n+1, 1) -> (m, 1)
    errors = predictions - y    # (m, 1)
    # 向量化计算平方误差和
    cost = (1/(2*m)) * np.sum(errors**2)
    return cost

# 让我们测试一下这个函数。首先，为我们的数据X添加偏置列。
X_b = np.c_[np.ones((m, 1)), X]  # 将形状(100,1)的X变为(100,2)，第一列为1
print("X_b的形状:", X_b.shape)  # 应该输出 (100, 2)

# 初始化参数theta，比如全设为0
theta_initial = np.zeros((2, 1))
print("初始theta:", theta_initial.T)  # 转置一下方便看

# 计算初始代价
initial_cost = compute_cost(X_b, y, theta_initial)
print(f"当theta=[0, 0]时，代价函数值: {initial_cost:.4f}")

如果一切正常，你会得到一个较大的代价函数值，因为用一条水平线（斜率为0）去拟合数据，误差当然很大。

2.2 实现梯度下降算法

代价函数告诉我们当前参数有多“差”，而梯度下降（Gradient Descent） 则告诉我们如何改进参数。其核心思想是沿着代价函数下降最快的方向（负梯度方向）更新参数。

更新规则（同时更新所有参数）： $$\theta_j := \theta_j - \alpha \frac{1}{m} \sum_{i=1}^{m} (\hat{y}^{(i)} - y^{(i)}) X_j^{(i)} \quad \text{for } j = 0, 1$$

其中 $\alpha$ 是学习率（Learning Rate），控制着每一步更新的幅度。

def gradient_descent(X, y, theta, alpha, num_iters):
    """
    执行批量梯度下降以学习参数theta。
    参数：
    X : 特征矩阵（已添加偏置列），形状 (m, n+1)
    y : 目标向量，形状 (m, 1)
    theta : 初始参数向量，形状 (n+1, 1)
    alpha : 学习率
    num_iters : 迭代次数
    返回：
    theta : 学习到的参数，形状 (n+1, 1)
    cost_history : 每次迭代的代价函数值列表，用于可视化
    """
    m = len(y)
    cost_history = []
    
    for i in range(num_iters):
        # 计算预测值
        predictions = X.dot(theta)  # (m, 1)
        # 计算误差
        errors = predictions - y    # (m, 1)
        # 计算梯度。注意：这是向量化形式，一次性计算所有theta的梯度。
        # X.T.dot(errors) 得到形状 (n+1, 1) 的梯度向量
        gradients = (1/m) * X.T.dot(errors)
        # 更新参数
        theta = theta - alpha * gradients
        # 记录本次迭代的代价
        cost = compute_cost(X, y, theta)
        cost_history.append(cost)
        
        # 每100次迭代打印一次进度（可选）
        if i % 100 == 0:
            print(f"Iteration {i}: Cost = {cost:.6f}")
    
    return theta, cost_history

# 设置超参数并运行梯度下降
alpha = 0.1  # 学习率
num_iters = 1000  # 迭代次数

# 再次使用初始化为0的参数
theta_initial = np.zeros((2, 1))
theta_learned, cost_hist = gradient_descent(X_b, y, theta_initial, alpha, num_iters)

print(f"\n学习到的参数 theta: {theta_learned.flatten()}")
print(f"真实的生成参数是: 截距=1.0, 斜率=2.5")

运行后，你应该看到代价函数随着迭代不断下降，最终学习到的参数应该接近我们生成数据时使用的真实参数（1.0和2.5）。由于我们加入了噪声，结果不会完全一致，但会非常接近。

2.3 可视化学习过程

理解算法如何工作，可视化是最有力的工具。我们来画两张图：一是代价函数的下降曲线，二是梯度下降过程中拟合直线的演变。

# 图1：代价函数随迭代次数的变化
plt.figure(figsize=(12, 5))

plt.subplot(1, 2, 1)
plt.plot(range(num_iters), cost_hist, 'b-', linewidth=2)
plt.xlabel('Iteration Number')
plt.ylabel('Cost J')
plt.title('Gradient Descent: Cost vs. Iterations')
plt.grid(True)

# 图2：最终拟合结果
plt.subplot(1, 2, 2)
plt.scatter(X, y, alpha=0.7, label='Training Data')
# 生成预测线
X_range = np.array([[0], [2.5]])  # 取X的最小和最大值范围
X_range_b = np.c_[np.ones((2, 1)), X_range]  # 添加偏置列
y_predict = X_range_b.dot(theta_learned)
plt.plot(X_range, y_predict, 'r-', linewidth=3, label=f'Fit: y={theta_learned[1][0]:.2f}x + {theta_learned[0][0]:.2f}')
plt.xlabel('Feature X')
plt.ylabel('Target y')
plt.title('Linear Regression Fit')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

第一张图展示了优化过程是否健康。一个平滑、持续下降的曲线表明学习率设置合适。如果曲线震荡剧烈或上升，说明学习率可能太大了。如果下降极其缓慢，则学习率可能太小。

3. 关键技巧与深度解析

仅仅实现基础版本还不够。在实际应用中，有几个技巧能显著提升模型的性能和训练效率。这部分往往是理论课程中一笔带过，但实践中至关重要的。

3.1 特征缩放：为什么以及怎么做

当特征之间的尺度差异巨大时（例如，一个特征是房屋面积（100-300平方米），另一个特征是房间数（2-6）），梯度下降会变得非常低效。代价函数的等高线图会变得又高又窄，梯度下降需要很多次曲折的“之字形”移动才能到达最低点。

特征缩放就是将不同特征的值缩放到相近的范围内。两种最常用的方法是：

1. 均值归一化（Mean Normalization）: $X' = \frac{X - \mu}{\text{range}}$ 或 $\frac{X - \mu}{\sigma}$
2. 标准化（Standardization / Z-score Normalization）: $X' = \frac{X - \mu}{\sigma}$

其中 $\mu$ 是均值，$\sigma$ 是标准差，range是极差（最大值减最小值）。标准化通常更常用，因为它处理了数据的分布，使得特征值大致服从标准正态分布。

def feature_standardization(X):
    """
    对特征矩阵X的每一列（每个特征）进行标准化。
    参数：
    X : 特征矩阵，形状 (m, n)，每列是一个特征
    返回：
    X_norm : 标准化后的特征矩阵
    mu : 每个特征的均值，形状 (1, n)
    sigma : 每个特征的标准差，形状 (1, n)
    """
    mu = np.mean(X, axis=0)  # 沿样本轴（第0轴）求均值
    sigma = np.std(X, axis=0)
    # 防止除零，如果某个特征标准差为0（所有值相同），则设为1
    sigma = np.where(sigma == 0, 1, sigma)
    X_norm = (X - mu) / sigma
    return X_norm, mu, sigma

# 注意：对于我们的单特征例子，特征缩放效果不明显，但为了演示，我们实现它。
# 假设我们有一个新的多元特征数据集
np.random.seed(123)
m2 = 50
X_multi = np.random.randn(m2, 2) * np.array([10, 0.1]) + np.array([5, 2])  # 两个尺度差异巨大的特征
print("原始特征示例（前5行）:")
print(X_multi[:5])
print(f"\n特征1的均值: {np.mean(X_multi[:,0]):.2f}, 标准差: {np.std(X_multi[:,0]):.2f}")
print(f"特征2的均值: {np.mean(X_multi[:,1]):.2f}, 标准差: {np.std(X_multi[:,1]):.2f}")

X_norm, mu, sigma = feature_standardization(X_multi)
print("\n标准化后特征示例（前5行）:")
print(X_norm[:5])
print(f"\n标准化后特征1的均值: {np.mean(X_norm[:,0]):.4f}, 标准差: {np.std(X_norm[:,0]):.4f}")
print(f"标准化后特征2的均值: {np.mean(X_norm[:,1]):.4f}, 标准差: {np.std(X_norm[:,1]):.4f}")

重要提示：计算均值mu和标准差sigma时，必须只使用训练集数据。然后用训练集计算出的mu和sigma去标准化验证集和测试集。绝对不能用整个数据集（包含测试集）来计算这些统计量，否则就造成了数据泄露，会严重高估模型性能。

3.2 学习率的选择与调试

学习率 $\alpha$ 是梯度下降中最重要的超参数。如何选择？

• $\alpha$ 太小：收敛速度极慢，需要大量迭代。
• $\alpha$ 太大：代价函数可能不收敛，甚至发散（代价函数值越来越大）。

一个实用的调试方法是：尝试一系列呈指数增长或衰减的 $\alpha$ 值，例如 0.001, 0.003, 0.01, 0.03, 0.1, 0.3。观察代价函数下降曲线。

def plot_learning_rates(X, y, thetas_init, alphas, num_iters=50):
    """
    绘制不同学习率下代价函数的下降曲线。
    """
    plt.figure(figsize=(10, 6))
    for alpha in alphas:
        theta_temp = thetas_init.copy()
        cost_hist = []
        for i in range(num_iters):
            predictions = X.dot(theta_temp)
            errors = predictions - y
            gradients = (1/len(y)) * X.T.dot(errors)
            theta_temp = theta_temp - alpha * gradients
            cost_hist.append(compute_cost(X, y, theta_temp))
        plt.plot(range(num_iters), cost_hist, label=f'alpha={alpha}')
    
    plt.xlabel('Iterations')
    plt.ylabel('Cost')
    plt.title('Effect of Learning Rate on Convergence')
    plt.legend()
    plt.grid(True)
    plt.yscale('log')  # 使用对数坐标可以更清晰地看到差异
    plt.show()

# 为我们的简单数据（已添加偏置列）测试不同学习率
alphas_to_try = [0.001, 0.01, 0.1, 0.5, 1.0]
plot_learning_rates(X_b, y, np.zeros((2, 1)), alphas_to_try)

从图中，你可以直观地看到：过小的学习率（如0.001）下降缓慢；合适的学习率（如0.01, 0.1）平滑快速下降；过大的学习率（如0.5, 1.0）可能导致震荡甚至发散（代价上升）。

3.3 正规方程：另一种求解方式

对于线性回归，我们其实可以不使用迭代的梯度下降，而通过解析方法直接求出最优参数。这个方法叫做正规方程（Normal Equation）。

公式为：$\theta = (X^T X)^{-1} X^T y$

其推导来自于令代价函数对 $\theta$ 的导数为零。它的优点是无需选择学习率，无需迭代，一次计算得到最优解。缺点是当特征数量 $n$ 很大时（例如 > 10000），计算 $(X^T X)^{-1}$ 的逆矩阵会非常慢（时间复杂度约为 $O(n^3)$），而梯度下降在特征数量很大时依然表现良好。

def normal_equation(X, y):
    """
    使用正规方程求解线性回归参数。
    参数：
    X : 特征矩阵（已添加偏置列），形状 (m, n+1)
    y : 目标向量，形状 (m, 1)
    返回：
    theta : 最优参数，形状 (n+1, 1)
    """
    # 计算 theta = np.linalg.inv(X.T @ X) @ X.T @ y
    # 使用np.linalg.pinv求伪逆，数值上更稳定，即使X.T @ X不可逆也能处理
    theta = np.linalg.pinv(X.T.dot(X)).dot(X.T).dot(y)
    return theta

# 用正规方程求解我们模拟数据的最优参数
theta_normal_eq = normal_equation(X_b, y)
print(f"正规方程求得的 theta: {theta_normal_eq.flatten()}")
print(f"梯度下降求得的 theta: {theta_learned.flatten()}")
print(f"两者差异: {np.abs(theta_normal_eq - theta_learned).flatten()}")

你会发现，正规方程和梯度下降（在充分迭代后）得到的结果几乎相同。对于小数据集，正规方程非常方便快捷。

方法	优点	缺点	适用场景
梯度下降	当特征数量n很大时也能高效工作；适用于各种模型（逻辑回归、神经网络）	需要选择学习率α；需要多次迭代	大规模数据集（m很大），特征数量多（n大）
正规方程	无需选择学习率；无需迭代；直接得到解析解	需要计算$(X^T X)^{-1}$，复杂度$O(n^3)$；当n很大时非常慢；若$X^T X$不可逆需特殊处理	小规模数据集（m适中），特征数量少（n < 1000）

4. 从一元到多元：扩展你的模型

现实问题中，影响结果的因素很少只有一个。房价不仅取决于面积，还有房间数、房龄、地段等。这就需要多元线性回归。好消息是，我们之前写的代码几乎不需要改动就能支持多元！因为我们已经使用了向量化的形式。

假设现在我们有一个包含两个特征的数据集：$X_1$（面积）和 $X_2$（房间数）。我们的假设函数变为： $\hat{y} = \theta_0 + \theta_1 X_1 + \theta_2 X_2$

向量化表示中，$X$ 是一个 $m \times 3$ 的矩阵（第一列是1），$\theta$ 是一个 $3 \times 1$ 的向量。代价函数和梯度下降的公式在形式上完全一致。

让我们用代码演示一下：

# 生成一个二元特征的模拟数据集
np.random.seed(2023)
m_multi = 200
# 特征1：面积（单位：百平方英尺），范围在10到30之间
X1 = 10 + 20 * np.random.rand(m_multi, 1)
# 特征2：房间数，范围在2到6之间
X2 = 2 + 4 * np.random.rand(m_multi, 1)
# 真实关系：价格 = 30 + 5*面积 + 10*房间数 + 噪声
y_multi = 30 + 5*X1 + 10*X2 + np.random.randn(m_multi, 1) * 5

# 组合特征
X_raw = np.hstack([X1, X2])  # 形状 (200, 2)
print("原始特征矩阵形状:", X_raw.shape)

# 1. 特征缩放（非常重要！）
X_norm, mu_multi, sigma_multi = feature_standardization(X_raw)

# 2. 添加偏置列
X_multi_b = np.c_[np.ones((m_multi, 1)), X_norm]  # 形状 (200, 3)

# 3. 初始化参数 (现在是3个：theta0, theta1, theta2)
theta_init_multi = np.zeros((3, 1))

# 4. 运行梯度下降
alpha_multi = 0.1
iters_multi = 500
theta_multi_learned, cost_hist_multi = gradient_descent(X_multi_b, y_multi, theta_init_multi, alpha_multi, iters_multi)

print(f"\n学习到的参数 (标准化后特征空间): {theta_multi_learned.flatten()}")

现在得到的 $\theta$ 是针对标准化后的特征 $X'$ 的。如果我们想得到针对原始特征 $X$ 的模型参数，需要进行转换。回忆一下，标准化是 $X' = (X - \mu) / \sigma$。我们的模型是：

$\hat{y} = \theta_0' + \theta_1' X_1' + \theta_2' X_2' = \theta_0' + \theta_1' \frac{X_1 - \mu_1}{\sigma_1} + \theta_2' \frac{X_2 - \mu_2}{\sigma_2}$

将其整理成 $\hat{y} = \theta_0 + \theta_1 X_1 + \theta_2 X_2$ 的形式：

$\theta_1 = \theta_1' / \sigma_1$ $\theta_2 = \theta_2' / \sigma_2$ $\theta_0 = \theta_0' - \theta_1' \frac{\mu_1}{\sigma_1} - \theta_2' \frac{\mu_2}{\sigma_2}$

# 将标准化特征空间的参数转换回原始特征空间
theta1_original = theta_multi_learned[1] / sigma_multi[0]
theta2_original = theta_multi_learned[2] / sigma_multi[1]
theta0_original = theta_multi_learned[0] - theta_multi_learned[1] * mu_multi[0]/sigma_multi[0] - theta_multi_learned[2] * mu_multi[1]/sigma_multi[1]

theta_original = np.array([[theta0_original[0]], [theta1_original[0]], [theta2_original[0]]])
print(f"转换回原始特征空间的参数: {theta_original.flatten()}")
print(f"真实生成参数: [30, 5, 10]")

对比一下，转换后的参数应该非常接近我们生成数据时使用的真实参数 [30, 5, 10]。这个转换步骤在实际部署模型时必不可少，因为你预测新数据时，输入的是原始特征值。

实现多元线性回归后，一个很自然的问题是：如何评估模型的好坏？除了看代价函数，在训练集上我们还可以计算决定系数 $R^2$，它表示模型对目标变量方差的解释比例，越接近1越好。

def r2_score(y_true, y_pred):
    """
    计算决定系数 R^2。
    """
    ss_res = np.sum((y_true - y_pred) ** 2)
    ss_tot = np.sum((y_true - np.mean(y_true)) ** 2)
    r2 = 1 - (ss_res / ss_tot)
    return r2

# 计算训练集上的预测和R^2
y_pred_train = X_multi_b.dot(theta_multi_learned)
train_r2 = r2_score(y_multi, y_pred_train)
print(f"模型在训练集上的 R^2 分数: {train_r2:.4f}")

一个接近1的 $R^2$ 分数（例如 > 0.9）表明模型拟合得很好。但要注意，这只是在训练集上的表现。为了评估模型的泛化能力，我们必须使用未参与训练的验证集或测试集。这引出了机器学习中另一个核心概念——训练/验证/测试集划分，以及防止过拟合的正则化技术，这些将是构建健壮模型的下一个关键步骤。不过，基于我们当前的目标，你已经成功地将线性回归的理论完整地实现为了可运行的代码，并理解了从数据准备、算法实现到结果评估的完整链路。下次当你看到线性回归的公式时，脑海中浮现的将不再只是符号，而是一行行可以构建出预测世界的代码。