第3章：线性模型（Linear Models for Regression）
一、导读
3.1 本章学习目标
第3章是全书的核心技术章节之一，系统介绍了线性模型在回归问题中的应用。线性模型看似简单，但它蕴含了机器学习中最重要的思想：基函数展开、正则化、贝叶斯推断、证据近似。掌握本章内容，就掌握了理解后续更复杂模型（神经网络、核方法、高斯过程）的钥匙。
本章的核心学习目标包括：
- 掌握线性模型的基函数展开：如何将非线性问题转化为线性问题？常见的基函数有哪些？
- 理解最小二乘法的几何直觉：为什么最小二乘法等价于投影？正规方程的几何含义是什么？
- 掌握正则化的理论与实现：岭回归、Lasso回归的数学推导和区别。
- 理解贝叶斯线性回归：如何将贝叶斯方法应用于线性模型？后验分布如何计算？
- 掌握证据近似（Type-II Maximum Likelihood）：如何自动确定正则化系数？
- 理解等价核（Equivalent Kernel）：线性回归的预测可以理解为训练目标的加权平均。
3.2 为什么本章如此重要
线性模型是机器学习的"原子模型"。几乎所有复杂的模型都可以分解为线性模型的组合：神经网络是多层线性模型的非线性组合，核方法是无穷维线性模型，高斯过程是贝叶斯线性回归的对偶形式。
更重要的是，本章引入的贝叶斯线性回归是理解现代机器学习中不确定性量化的基础。在安全关键的应用（如自动驾驶、医疗诊断）中，模型不仅需要给出预测，还需要给出预测的置信度。贝叶斯线性回归提供了这种能力。
3.3 核心问题引导
在阅读本章时，请带着以下问题思考：
为什么线性模型可以处理非线性问题？基函数展开的本质是什么？
岭回归和Lasso回归的正则化效果有什么不同？为什么Lasso能产生稀疏解？
贝叶斯线性回归中，后验预测方差由哪两部分组成？各自的含义是什么？
证据近似（Type-II ML）是如何自动确定正则化系数的？
等价核揭示了线性回归的什么本质？它与核密度估计有什么联系？
当特征数大于样本数时，最小二乘法为什么会失效？如何解决？
二、核心概念详解
3.1 基函数线性模型
3.1.1 模型定义
线性模型的一般形式为：
$$y(\mathbf{x}, \mathbf{w}) = \sum_{j=0}^{M-1} w_j \phi_j(\mathbf{x}) = \mathbf{w}^T \boldsymbol{\phi}(\mathbf{x})$$
其中 $\phi_j(\mathbf{x})$ 是基函数（Basis Functions），$\mathbf{w}$ 是权重向量。
"线性"指的是模型对参数 $\mathbf{w}$ 是线性的，而不是对输入 $\mathbf{x}$ 是线性的。通过使用非线性基函数，线性模型可以拟合任意复杂的非线性函数。
3.1.2 常见基函数
多项式基函数：$\phi_j(x) = x^j$
高斯基函数：$\phi_j(x) = \exp\left\{-\frac{(x-\mu_j)^2}{2s^2}\right\}$
Sigmoid基函数：$\phi_j(x) = \sigma(a_j(x - b_j)) = \frac{1}{1 + \exp(-a_j(x-b_j))}$
import numpy as np
import matplotlib.pyplot as plt
def polynomial_basis(x, M):
"""多项式基函数"""
return np.column_stack([x ** j for j in range(M)])
def gaussian_basis(x, centers, width):
"""高斯基函数"""
return np.exp(-0.5 * ((x[:, None] - centers[None, :]) / width) ** 2)
def sigmoid_basis(x, centers, slopes):
"""Sigmoid基函数"""
return 1 / (1 + np.exp(-slopes[None, :] * (x[:, None] - centers[None, :])))
# 可视化不同基函数
x = np.linspace(-3, 3, 200)
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
# 多项式基函数
for M in range(5):
axes[0].plot(x, x ** M, label=f'x^{M}')
axes[0].set_title('Polynomial Basis')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
axes[0].set_ylim(-5, 5)
# 高斯基函数
centers = np.linspace(-2, 2, 5)
width = 0.5
for i, c in enumerate(centers):
axes[1].plot(x, np.exp(-0.5 * ((x - c) / width) ** 2), label=f'μ={c:.1f}')
axes[1].set_title('Gaussian Basis')
axes[1].legend(fontsize=8)
axes[1].grid(True, alpha=0.3)
# Sigmoid基函数
for i, c in enumerate(centers):
axes[2].plot(x, 1 / (1 + np.exp(-5 * (x - c))), label=f'b={c:.1f}')
axes[2].set_title('Sigmoid Basis')
axes[2].legend(fontsize=8)
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('basis_functions.png', dpi=150)
plt.close()3.1.3 设计矩阵与正规方程
给定训练数据 $\{(\mathbf{x}_n, t_n)\}_{n=1}^N$，定义设计矩阵 $\mathbf{\Phi}$：
$$\mathbf{\Phi} = \begin{pmatrix} \phi_0(\mathbf{x}_1) & \phi_1(\mathbf{x}_1) & \cdots & \phi_{M-1}(\mathbf{x}_1) \\ \phi_0(\mathbf{x}_2) & \phi_1(\mathbf{x}_2) & \cdots & \phi_{M-1}(\mathbf{x}_2) \\ \vdots & \vdots & \ddots & \vdots \\ \phi_0(\mathbf{x}_N) & \phi_1(\mathbf{x}_N) & \cdots & \phi_{M-1}(\mathbf{x}_N) \end{pmatrix}$$
最小二乘法的解（正规方程）为：
$$\mathbf{w}_{ML} = (\mathbf{\Phi}^T \mathbf{\Phi})^{-1} \mathbf{\Phi}^T \mathbf{t}$$
几何直觉：$\mathbf{\Phi} \mathbf{w}_{ML}$ 是目标向量 $\mathbf{t}$ 在设计矩阵列空间上的正交投影。残差向量 $\mathbf{t} - \mathbf{\Phi} \mathbf{w}_{ML}$ 与列空间正交。
def linear_regression_basis(x_train, t_train, x_test, basis_func, **basis_kwargs):
"""使用基函数的线性回归"""
Phi_train = basis_func(x_train, **basis_kwargs)
Phi_test = basis_func(x_test, **basis_kwargs)
# 正规方程
w_ml = np.linalg.inv(Phi_train.T @ Phi_train) @ Phi_train.T @ t_train
# 预测
y_test = Phi_test @ w_ml
return y_test, w_ml
# 示例：使用高斯基函数拟合正弦函数
np.random.seed(42)
N = 25
x_train = np.sort(np.random.uniform(-3, 3, N))
t_train = np.sin(x_train) + np.random.normal(0, 0.2, N)
x_test = np.linspace(-3, 3, 200)
centers = np.linspace(-3, 3, 10)
y_test, w = linear_regression_basis(x_train, t_train, x_test, gaussian_basis,
centers=centers, width=0.5)
print(f"训练误差: {np.sqrt(np.mean((gaussian_basis(x_train, centers, 0.5) @ w - t_train)**2)):.4f}")
print(f"权重范围: [{w.min():.3f}, {w.max():.3f}]")3.2 正则化
3.2.1 岭回归（L2正则化）
岭回归在误差函数中添加 $L_2$ 惩罚项：
$$\tilde{E}(\mathbf{w}) = \frac{1}{2} \sum_{n=1}^{N} \{y(\mathbf{x}_n, \mathbf{w}) - t_n\}^2 + \frac{\lambda}{2} \|\mathbf{w}\|^2$$
解为：
$$\mathbf{w}_{Ridge} = (\mathbf{\Phi}^T \mathbf{\Phi} + \lambda \mathbf{I})^{-1} \mathbf{\Phi}^T \mathbf{t}$$
岭回归的效果：
- 使设计矩阵 $\mathbf{\Phi}^T \mathbf{\Phi}$ 的条件数改善，数值更稳定
- 收缩权重值，防止过拟合
- 所有权重都趋向于非零的小值
3.2.2 Lasso回归（L1正则化）
Lasso回归使用 $L_1$ 惩罚项：
$$\tilde{E}(\mathbf{w}) = \frac{1}{2} \sum_{n=1}^{N} \{y(\mathbf{x}_n, \mathbf{w}) - t_n\}^2 + \lambda \|\mathbf{w}\|_1$$
Lasso的关键特性是产生稀疏解——部分权重精确为零。这使得Lasso同时进行特征选择和参数估计。
Lasso产生稀疏解的几何直觉：$L_1$ 约束区域是超立方体，其尖角使得误差函数的等高线更容易在坐标轴上与之相切。
def ridge_regression(Phi, t, lam):
"""岭回归"""
M = Phi.shape[1]
I = np.eye(M)
I[0, 0] = 0 # 不正则化偏置
w = np.linalg.inv(Phi.T @ Phi + lam * I) @ Phi.T @ t
return w
def lasso_regression_ista(Phi, t, lam, n_iter=1000, lr=0.001):
"""Lasso回归（ISTA算法）"""
M = Phi.shape[1]
w = np.zeros(M)
for _ in range(n_iter):
# 梯度步
grad = Phi.T @ (Phi @ w - t)
w = w - lr * grad
# 近端算子（软阈值）
w = np.sign(w) * np.maximum(np.abs(w) - lr * lam, 0)
w[0] = np.mean(t - Phi[:, 1:] @ w[1:]) # 偏置不正则化
return w
# 比较岭回归和Lasso
np.random.seed(42)
N, M = 50, 20
x = np.random.randn(N, M)
true_w = np.zeros(M)
true_w[:5] = [2, -1, 3, 0.5, -2] # 只有前5个特征相关
t = x @ true_w + 0.5 * np.random.randn(N)
Phi = np.column_stack([np.ones(N), x])
# 岭回归
lam = 1.0
w_ridge = ridge_regression(Phi, t, lam)
# Lasso
w_lasso = lasso_regression_ista(Phi, t, lam=0.1, n_iter=5000, lr=0.001)
print("真实权重:", np.round(np.concatenate([[0], true_w]), 2))
print("岭回归权重:", np.round(w_ridge, 2))
print("Lasso权重:", np.round(w_lasso, 2))
print(f"\nLasso非零权重数: {np.sum(np.abs(w_lasso) > 0.01)} / {M+1}")3.3 贝叶斯线性回归
3.3.1 模型设定
贝叶斯线性回归的完整概率模型：
- 似然：$p(\mathbf{t}|\mathbf{w}, \beta) = \prod_{n=1}^{N} \mathcal{N}(t_n|\mathbf{w}^T \boldsymbol{\phi}(\mathbf{x}_n), \beta^{-1})$
- 先验：$p(\mathbf{w}|\alpha) = \mathcal{N}(\mathbf{w}|\mathbf{0}, \alpha^{-1} \mathbf{I})$
- 后验：$p(\mathbf{w}|\mathbf{t}, \alpha, \beta) = \mathcal{N}(\mathbf{w}|\mathbf{m}_N, \mathbf{S}_N)$
后验参数：
- $\mathbf{S}_N^{-1} = \alpha \mathbf{I} + \beta \mathbf{\Phi}^T \mathbf{\Phi}$
- $\mathbf{m}_N = \beta \mathbf{S}_N \mathbf{\Phi}^T \mathbf{t}$
注意：$\alpha$ 和 $\lambda$ 的关系——$\alpha/\beta = \lambda$。先验精度 $\alpha$ 越大，正则化越强。
3.3.2 预测分布
贝叶斯预测考虑了参数的不确定性：
$$p(t|\mathbf{x}, \mathbf{t}, \alpha, \beta) = \int p(t|\mathbf{x}, \mathbf{w}, \beta) p(\mathbf{w}|\mathbf{t}, \alpha, \beta) d\mathbf{w} = \mathcal{N}(t|\mathbf{m}_N^T \boldsymbol{\phi}(\mathbf{x}), \sigma_N^2(\mathbf{x}))$$
预测方差：
$$\sigma_N^2(\mathbf{x}) = \frac{1}{\beta} + \boldsymbol{\phi}(\mathbf{x})^T \mathbf{S}_N \boldsymbol{\phi}(\mathbf{x})$$
- $\frac{1}{\beta}$：数据噪声（不可约）
- $\boldsymbol{\phi}(\mathbf{x})^T \mathbf{S}_N \boldsymbol{\phi}(\mathbf{x})$：模型不确定性（随数据增多而减小）
def bayesian_linear_regression(x_train, t_train, x_test, basis_func, alpha, beta, **basis_kwargs):
"""贝叶斯线性回归"""
Phi_train = basis_func(x_train, **basis_kwargs)
Phi_test = basis_func(x_test, **basis_kwargs)
M = Phi_train.shape[1]
# 后验参数
S_N_inv = alpha * np.eye(M) + beta * Phi_train.T @ Phi_train
S_N = np.linalg.inv(S_N_inv)
m_N = beta * S_N @ Phi_train.T @ t_train
# 预测分布
pred_mean = Phi_test @ m_N
pred_var = np.array([1/beta + Phi_test[i] @ S_N @ Phi_test[i] for i in range(len(x_test))])
return pred_mean, pred_var, m_N, S_N
# 完整示例
np.random.seed(42)
N = 15
x_train = np.sort(np.random.uniform(-3, 3, N))
t_train = np.sin(x_train) + np.random.normal(0, 0.3, N)
x_test = np.linspace(-4, 4, 200)
centers = np.linspace(-3, 3, 8)
width = 0.8
alpha = 2.0
beta = 10.0
pred_mean, pred_var, m_N, S_N = bayesian_linear_regression(
x_train, t_train, x_test, gaussian_basis, alpha, beta,
centers=centers, width=width
)
pred_std = np.sqrt(pred_var)
# 可视化
fig, ax = plt.subplots(figsize=(10, 6))
ax.scatter(x_train, t_train, c='red', s=50, zorder=5, label='Training data')
ax.plot(x_test, pred_mean, 'b-', linewidth=2, label='Predictive mean')
ax.fill_between(x_test, pred_mean - 2*pred_std, pred_mean + 2*pred_std,
alpha=0.2, color='blue', label='95% CI')
ax.plot(x_test, np.sin(x_test), 'g--', linewidth=2, label='True function')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('Bayesian Linear Regression with Gaussian Basis')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('bayesian_linear_regression.png', dpi=150)
plt.close()
print(f"后验均值范围: [{m_N.min():.3f}, {m_N.max():.3f}]")
print(f"预测标准差范围: [{pred_std.min():.3f}, {pred_std.max():.3f}]")3.4 证据近似（Type-II Maximum Likelihood）
3.4.1 动机
贝叶斯线性回归需要指定超参数 $\alpha$（先验精度）和 $\beta$（噪声精度）。手动选择这些超参数既繁琐又不够客观。证据近似提供了一种数据驱动的方法来确定超参数。
3.4.2 边缘似然
边缘似然（Evidence）是通过积分掉参数 $\mathbf{w}$ 得到的：
$$p(\mathbf{t}|\alpha, \beta) = \int p(\mathbf{t}|\mathbf{w}, \beta) p(\mathbf{w}|\alpha) d\mathbf{w}$$
对于线性高斯模型，边缘似然有解析解：
$$\ln p(\mathbf{t}|\alpha, \beta) = \frac{M}{2} \ln \alpha + \frac{N}{2} \ln \beta - E(\mathbf{m}_N) - \frac{1}{2} \ln |\mathbf{S}_N^{-1}| - \frac{N}{2} \ln(2\pi)$$
其中 $E(\mathbf{m}_N) = \frac{\beta}{2} \|\mathbf{t} - \mathbf{\Phi} \mathbf{m}_N\|^2 + \frac{\alpha}{2} \mathbf{m}_N^T \mathbf{m}_N$。
3.4.3 超参数的优化
通过对边缘似然关于 $\alpha$ 和 $\beta$ 求导并令其为零，得到迭代更新公式：
- $\alpha = \frac{\gamma}{\mathbf{m}_N^T \mathbf{m}_N}$
- $\beta = \frac{N - \gamma}{\|\mathbf{t} - \mathbf{\Phi} \mathbf{m}_N\|^2}$
- $\gamma = \sum_{i=1}^{M} \frac{\lambda_i}{\alpha + \lambda_i}$（$\lambda_i$ 是 $\beta \mathbf{\Phi}^T \mathbf{\Phi}$ 的特征值）
def evidence_approximation(x_train, t_train, basis_func, n_iter=100, **basis_kwargs):
"""证据近似确定超参数"""
Phi = basis_func(x_train, **basis_kwargs)
N, M = Phi.shape
# 初始化
alpha = 0.01
beta = 1.0
for iteration in range(n_iter):
# 后验参数
S_N_inv = alpha * np.eye(M) + beta * Phi.T @ Phi
S_N = np.linalg.inv(S_N_inv)
m_N = beta * S_N @ Phi.T @ t_train
# 特征值
eigenvalues = np.linalg.eigvalsh(beta * Phi.T @ Phi)
# gamma
gamma = np.sum(eigenvalues / (alpha + eigenvalues))
# 更新超参数
alpha_new = gamma / (m_N @ m_N + 1e-10)
beta_new = (N - gamma) / (np.sum((t_train - Phi @ m_N) ** 2) + 1e-10)
# 检查收敛
if abs(alpha_new - alpha) < 1e-6 and abs(beta_new - beta) < 1e-6:
print(f"收敛于第 {iteration+1} 次迭代")
break
alpha = alpha_new
beta = beta_new
# 边缘似然
E_mN = beta/2 * np.sum((t_train - Phi @ m_N)**2) + alpha/2 * m_N @ m_N
log_evidence = M/2 * np.log(alpha) + N/2 * np.log(beta) - E_mN \
- 0.5 * np.log(np.linalg.det(S_N_inv)) - N/2 * np.log(2*np.pi)
print(f"最优 α = {alpha:.4f}")
print(f"最优 β = {beta:.4f}")
print(f"等效 λ = α/β = {alpha/beta:.4f}")
print(f"对数边缘似然 = {log_evidence:.4f}")
return alpha, beta, m_N, S_N
# 示例
np.random.seed(42)
N = 25
x_train = np.sort(np.random.uniform(-3, 3, N))
t_train = np.sin(x_train) + np.random.normal(0, 0.3, N)
centers = np.linspace(-3, 3, 10)
alpha, beta, m_N, S_N = evidence_approximation(x_train, t_train, gaussian_basis,
n_iter=100, centers=centers, width=0.5)3.5 等价核
3.5.1 概念
贝叶斯线性回归的预测可以写成训练目标的加权平均：
$$y(\mathbf{x}) = \mathbf{m}_N^T \boldsymbol{\phi}(\mathbf{x}) = \beta \boldsymbol{\phi}(\mathbf{x})^T \mathbf{S}_N \mathbf{\Phi}^T \mathbf{t} = \sum_{n=1}^{N} k(\mathbf{x}, \mathbf{x}_n) t_n$$
其中等价核（Equivalent Kernel）为：
$$k(\mathbf{x}, \mathbf{x}') = \beta \boldsymbol{\phi}(\mathbf{x})^T \mathbf{S}_N \boldsymbol{\phi}(\mathbf{x}')$$
等价核揭示了线性回归的本质：预测是训练目标的平滑。等价核决定了每个训练点对预测的贡献权重。
def equivalent_kernel(x_train, x_test, basis_func, alpha, beta, **basis_kwargs):
"""计算等价核"""
Phi_train = basis_func(x_train, **basis_kwargs)
Phi_test = basis_func(x_test, **basis_kwargs)
M = Phi_train.shape[1]
S_N_inv = alpha * np.eye(M) + beta * Phi_train.T @ Phi_train
S_N = np.linalg.inv(S_N_inv)
# 等价核矩阵
K = beta * Phi_test @ S_N @ Phi_train.T
return K
# 可视化等价核
np.random.seed(42)
N = 20
x_train = np.sort(np.random.uniform(-3, 3, N))
x_test_point = np.array([0.0]) # 关注 x=0 处的等价核
x_test_all = np.linspace(-3, 3, 200)
centers = np.linspace(-3, 3, 10)
alpha = 2.0
beta = 10.0
K = equivalent_kernel(x_train, x_test_all, gaussian_basis, alpha, beta,
centers=centers, width=0.5)
# 找到最接近 x=0 的测试点
idx = np.argmin(np.abs(x_test_all))
kernel_weights = K[idx, :]
fig, ax = plt.subplots(figsize=(10, 5))
ax.stem(x_train, kernel_weights, linefmt='b-', markerfmt='bo', basefmt='k-')
ax.set_xlabel('Training point x_n')
ax.set_ylabel('Kernel weight k(0, x_n)')
ax.set_title('Equivalent Kernel at x=0')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('equivalent_kernel.png', dpi=150)
plt.close()
print(f"等价核权重之和: {kernel_weights.sum():.4f}")
print(f"最大权重: {kernel_weights.max():.4f} (对应 x={x_train[np.argmax(kernel_weights)]:.3f})")三、实战案例
3.1 多项式回归的偏差-方差权衡
def bias_variance_tradeoff():
"""研究多项式回归的偏差-方差权衡"""
np.random.seed(42)
N_train = 25
N_test = 200
n_trials = 100
x_test = np.linspace(-3, 3, N_test)
true_function = np.sin(x_test)
M_values = range(1, 13)
biases_sq = []
variances = []
irreducible = 0.09 # sigma^2 = 0.3^2
for M in M_values:
predictions = []
for _ in range(n_trials):
x_train = np.sort(np.random.uniform(-3, 3, N_train))
t_train = np.sin(x_train) + np.random.normal(0, 0.3, N_train)
Phi = polynomial_basis(x_train, M)
w = np.linalg.inv(Phi.T @ Phi + 0.01*np.eye(M)) @ Phi.T @ t_train
Phi_test = polynomial_basis(x_test, M)
predictions.append(Phi_test @ w)
predictions = np.array(predictions)
mean_pred = predictions.mean(axis=0)
bias_sq = np.mean((mean_pred - true_function) ** 2)
variance = np.mean(np.var(predictions, axis=0))
biases_sq.append(bias_sq)
variances.append(variance)
print(f"{'M':>3} {'Bias²':>10} {'Variance':>10} {'Total':>10}")
print("-" * 40)
for M, b, v in zip(M_values, biases_sq, variances):
print(f"{M:3d} {b:10.4f} {v:10.4f} {b+v+irreducible:10.4f}")
bias_variance_tradeoff()3.2 正则化路径
def regularization_path():
"""绘制正则化路径：权重随正则化系数的变化"""
from sklearn.linear_model import Ridge, Lasso
np.random.seed(42)
N, D = 100, 15
X = np.random.randn(N, D)
true_w = np.zeros(D)
true_w[:5] = [3, -2, 1.5, 0.5, -1]
t = X @ true_w + 0.5 * np.random.randn(N)
# 正则化路径
alphas = np.logspace(-3, 3, 50)
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 岭回归路径
ridge_coefs = []
for alpha in alphas:
ridge = Ridge(alpha=alpha, fit_intercept=True)
ridge.fit(X, t)
ridge_coefs.append(ridge.coef_)
ridge_coefs = np.array(ridge_coefs)
for j in range(D):
axes[0].plot(np.log10(alphas), ridge_coefs[:, j],
alpha=0.7, linewidth=1.5 if j < 5 else 0.5)
axes[0].set_xlabel('log10(α)')
axes[0].set_ylabel('Coefficient')
axes[0].set_title('Ridge Regression Path')
axes[0].grid(True, alpha=0.3)
axes[0].axhline(y=0, color='k', linestyle='--', alpha=0.3)
# Lasso路径
lasso_coefs = []
for alpha in alphas * 0.01:
lasso = Lasso(alpha=alpha, fit_intercept=True, max_iter=5000)
lasso.fit(X, t)
lasso_coefs.append(lasso.coef_)
lasso_coefs = np.array(lasso_coefs)
for j in range(D):
axes[1].plot(np.log10(alphas * 0.01), lasso_coefs[:, j],
alpha=0.7, linewidth=1.5 if j < 5 else 0.5)
axes[1].set_xlabel('log10(α)')
axes[1].set_ylabel('Coefficient')
axes[1].set_title('Lasso Regression Path')
axes[1].grid(True, alpha=0.3)
axes[1].axhline(y=0, color='k', linestyle='--', alpha=0.3)
plt.tight_layout()
plt.savefig('regularization_path.png', dpi=150)
plt.close()
print("岭回归: 所有权重随α增大平滑收缩到零")
print("Lasso: 权重在某个α值后精确变为零（特征选择）")
regularization_path()3.3 贝叶斯模型比较
def bayesian_model_comparison():
"""使用边缘似然比较不同复杂度的模型"""
np.random.seed(42)
N = 30
x = np.sort(np.random.uniform(-3, 3, N))
t = np.sin(x) + np.random.normal(0, 0.3, N)
M_values = range(1, 10)
log_evidences = []
for M in M_values:
Phi = polynomial_basis(x, M)
# 证据近似
alpha = 0.01
beta = 1.0
for _ in range(50):
S_inv = alpha * np.eye(M) + beta * Phi.T @ Phi
S = np.linalg.inv(S_inv)
m = beta * S @ Phi.T @ t
eigenvalues = np.linalg.eigvalsh(beta * Phi.T @ Phi)
gamma = np.sum(eigenvalues / (alpha + eigenvalues))
alpha = gamma / (m @ m + 1e-10)
beta = (N - gamma) / (np.sum((t - Phi @ m)**2) + 1e-10)
# 对数边缘似然
E_m = beta/2 * np.sum((t - Phi @ m)**2) + alpha/2 * m @ m
log_evd = M/2 * np.log(alpha) + N/2 * np.log(beta) - E_m \
- 0.5 * np.log(np.linalg.det(S_inv)) - N/2 * np.log(2*np.pi)
log_evidences.append(log_evd)
print(f"M={M}: ln(evidence)={log_evd:.2f}, α={alpha:.4f}, β={beta:.4f}")
best_M = M_values[np.argmax(log_evidences)]
print(f"\n贝叶斯模型选择: 最优 M = {best_M}")
bayesian_model_comparison()四、本章小结
第3章系统介绍了线性模型在回归问题中的完整理论框架：
基函数展开：通过将输入映射到特征空间，线性模型可以处理非线性问题。多项式基函数、高斯基函数和Sigmoid基函数是最常用的选择。
正则化：岭回归（L2）平滑收缩权重，Lasso（L1）产生稀疏解并进行特征选择。正则化系数控制模型复杂度与泛化能力的权衡。
贝叶斯线性回归：将参数视为随机变量，提供完整的后验分布和预测不确定性量化。预测方差由数据噪声和模型不确定性两部分组成。
证据近似：通过最大化边缘似然自动确定超参数 $\alpha$ 和 $\beta$，避免了手动调参的繁琐。
等价核：揭示了线性回归预测的本质——训练目标的加权平滑。等价核的形状决定了模型的平滑程度。
线性模型的理论框架为后续章节奠定了基础。第4章将线性模型扩展到分类问题，第5章的神经网络可以理解为多层线性模型的非线性组合，第6章的核方法是无穷维线性模型的对偶形式，第8章的图模型利用线性高斯的性质进行高效推断。
五、思考与练习
思考题
为什么岭回归的解总是存在的（$\mathbf{\Phi}^T \mathbf{\Phi} + \lambda \mathbf{I}$ 总是可逆的），而最小二乘法的解可能不存在？
Lasso为什么能产生稀疏解？请从优化问题的几何角度和次梯度角度分别解释。
贝叶斯线性回归的预测方差在训练数据稀疏的区域会变大。这种性质在实际应用中有什么意义？
证据近似中，为什么边缘似然关于超参数 $\alpha$ 的优化不会过拟合？
等价核与核方法中的核函数有什么联系和区别？
编程练习
弹性网络：实现弹性网络（Elastic Net），结合L1和L2正则化：$E = \frac{1}{2}\|\mathbf{\Phi}\mathbf{w} - \mathbf{t}\|^2 + \lambda_1 \|\mathbf{w}\|_1 + \lambda_2 \|\mathbf{w}\|^2$。比较不同 $\lambda_1/\lambda_2$ 比例下的稀疏性和预测精度。
贝叶斯学习曲线：实现贝叶斯线性回归，随着训练数据量的增加，绘制后验预测区间的变化。验证"数据越多，不确定性越小"的直觉。
高斯基函数宽度的影响：固定高斯基函数的中心位置，改变宽度参数 $s$，观察模型拟合效果的变化。讨论宽度参数与正则化系数的关系。