\[ \begin{align}\begin{aligned}\newcommand{\ba}{\boldsymbol{a}} \newcommand{\bb}{\boldsymbol{b}} \newcommand{\be}{\boldsymbol{e}} \newcommand{\bq}{\boldsymbol{q}} \newcommand{\bk}{\boldsymbol{k}} \newcommand{\bw}{\boldsymbol{w}} \newcommand{\bx}{\boldsymbol{x}} \newcommand{\by}{\boldsymbol{y}} \newcommand{\bz}{\boldsymbol{z}} \newcommand{\bd}{\boldsymbol{d}} \newcommand{\bv}{\boldsymbol{v}} \newcommand{\bs}{\boldsymbol{s}}\\\newcommand{\btheta}{\boldsymbol{\theta}} \newcommand{\bbeta}{\boldsymbol{\beta}} \newcommand{\bgamma}{\boldsymbol{\gamma}} \newcommand{\bsigma}{\boldsymbol{\sigma}} \newcommand{\md}{\mbox{d}} \newcommand{\bmu}{\boldsymbol{\mu}} \newcommand{\bone}{\boldsymbol{1}} \newcommand{\bzero}{\boldsymbol{0}} \newcommand{\bepsilon}{\boldsymbol{\epsilon}} \newcommand{\bphi}{\boldsymbol{\phi}} \newcommand{\bh}{\boldsymbol{h}} \newcommand{\bc}{\boldsymbol{c}} \newcommand{\br}{\boldsymbol{r}} \newcommand{\bQ}{\boldsymbol{Q}} \newcommand{\bK}{\boldsymbol{K}} \newcommand{\bV}{\boldsymbol{V}} \newcommand{\bSigma}{\boldsymbol{\Sigma}} \newcommand{\bg}{\boldsymbol{g}} \newcommand{\bxi}{\boldsymbol{\xi}} \newcommand{\bvarepsilon}{\boldsymbol{\varepsilon}} \newcommand{\bdelta}{\boldsymbol{\delta}} \newcommand{\bq}{\boldsymbol{q}} \newcommand{\bk}{\boldsymbol{k}} \newcommand{\bJ}{\boldsymbol{J}} \newcommand{\bp}{\boldsymbol{p}} \newcommand{\bi}{\boldsymbol{i}} \newcommand{\bo}{\boldsymbol{o}} \newcommand{\bE}{\boldsymbol{E}} \newcommand{\bH}{\boldsymbol{H}} \newcommand{\bL}{\boldsymbol{L}} \newcommand{\bu}{\boldsymbol{u}} \newcommand{\bLambda}{\boldsymbol{\Lambda}} \newcommand{\trans}{^{\rm\scriptsize T}} \newcommand{\var}{\mathrm{var}}\\\newcommand{\bA}{\boldsymbol{A}} \newcommand{\bB}{\boldsymbol{B}} \newcommand{\bC}{\boldsymbol{C}} \newcommand{\bD}{\boldsymbol{D}} \newcommand{\bG}{\boldsymbol{G}} \newcommand{\bI}{\boldsymbol{I}} \newcommand{\bM}{\boldsymbol{M}} \newcommand{\bP}{\boldsymbol{P}} \newcommand{\bS}{\boldsymbol{S}} \newcommand{\bU}{\boldsymbol{U}} \newcommand{\bW}{\boldsymbol{W}} \newcommand{\bX}{\boldsymbol{X}} \newcommand{\bY}{\boldsymbol{Y}} \newcommand{\bZ}{\boldsymbol{Z}} \newcommand{\cotp}{\textcolor[RGB]{48,209,88}{TP}} \newcommand{\cotn}{\textcolor[RGB]{100,210,255}{TN}} \newcommand{\cofp}{\textcolor[RGB]{94,92,230}{FP}} \newcommand{\cofn}{\textcolor[RGB]{191,90,242}{FN}}\\\newcommand{\numcotp}{\textcolor[RGB]{48,209,88}{50}} \newcommand{\numcotn}{\textcolor[RGB]{100,210,255}{30}} \newcommand{\numcofp}{\textcolor[RGB]{94,92,230}{10}} \newcommand{\numcofn}{\textcolor[RGB]{191,90,242}{10}} \DeclareMathOperator*{\argmin}{arg\,min}\end{aligned}\end{align} \]

一元变量的线性回归

学习目标与记号

1. 准确解释可复现实验、线性回归与最小二乘估计；

2. 根据公式和张量维度分析蒙特卡洛模拟，并识别常见实现错误；

3. 用 Python 实现或验证向量化，通过受控实验解释结果。

　　本节沿用统一记号：普通小写字母表示标量，粗体小写字母表示向量，粗体大写字母表示矩阵或高阶张量；样本或时间编号写作下标；转置写作 \(\trans\)⁠。除非另有说明，批量样本按行存放。正文与练习中的程序都应同时检查数值结果和数组维度。本节练习的参考答案见 Python 线性回归实验与向量化答案⁠。

　　在本节中，我们将依托一元变量的线性回归模型，介绍如何对一个统计模型进行蒙特卡洛模拟，并对模型参数的估计量进行统计推断。在这个过程中，我们将分析步骤拆分成若干独立的函数，并利用 NumPy 的基本命令实现编程。

随机种子的设定

　　在讨论一元线性回归之前，先介绍 随机种子 （random seed）。模拟研究按照预先设定的真实模型产生样本；固定伪随机数生成器的种子，可以使相同程序重复生成同一组随机数，从而便于复现实验结果。若不固定种子，每次运行得到的样本和数值结果通常会有所不同。需要注意，固定随机种子只能保证给定软件环境下的随机过程可复现，并不能证明模型设定或研究结论正确。

　　首先，我们通过两个简单的例子展示 随机种子 的设定，以及其在随机数生成中的作用。下面这段代码中，我们利用 np.random.normal() 产生两个随机数。由于我们并没有预先设定随机种子，所产生的两个随机数并不相同。

import numpy as np #加载库

a = np.random.normal(1)
b = np.random.normal(1)
print(a,b)

0.4852213282828558 0.7160725818075863

　　第一段程序没有重置随机状态，a 和 b 是随机序列中连续抽取的两个标量，具体数值每次运行可能改变，二者也通常不同。接下来利用 np.random.seed() 在每次抽样前都把生成器重置到同一状态；预期第二段程序打印两个相同的数。

np.random.seed(1234)
a = np.random.normal(1)
np.random.seed(1234)
b = np.random.normal(1)
print(a,b)

1.471435163732493 1.471435163732493

　　两次都以种子 1234 开始并执行同一条正态随机数命令，因此 a 与 b 相同。这说明随机种子有助于复现给定软件环境中的随机过程，但频繁重置全局生成器可能意外重复样本；较大的程序通常更适合创建并传递独立的随机数生成器。设定随机种子对模拟数据分析结果的复现性非常重要。

备注

　　在使用 Python 进行模拟实验时，我们可能会用到其他的库产生随机数，例如 random 库。对于 random 库而言，我们则需要使用 random.seed() 设定随机种子。

生成训练集

真实模型假设

　　在本节中，我们考虑以下线性回归模型：

(1)\[y_i=b_0+w_0x_i+\epsilon_i, \qquad i=1,\ldots,n,\]

其中，真实参数 \(b_0=w_0=1\)⁠，特征 \(x_i\sim\mathcal{N}(2,2^2)\)⁠，误差项 \(\epsilon_i\sim\mathcal{N}(0,1)\)⁠，且 \(x_i\) 与 \(\epsilon_i\) 相互独立；\(n\) 为样本量，\(y_i\) 为第 \(i\) 个样本的标签；\(\mathcal{N}(\mu,\sigma^2)\) 表示均值为 \(\mu\)⁠、方差为 \(\sigma^2\) 的正态分布，下标 0 表示用于生成数据的真实参数。

备注

　　真实模型与分析模型。 真实模型描述数据实际的产生过程，分析模型则是研究者为了估计参数、解释变量之间的关系或预测结果而建立的模型。在模拟实验中，研究者预先确定真实模型，因此可以比较估计结果与已知参数，检查分析方法能否正确恢复数据中的规律。在真实数据分析中，真实模型通常无法完全获知，分析模型只能根据研究目的、专业知识和已有数据提出，因而往往是真实模型的近似。即使分析模型能够较好地拟合已有数据，也不能据此断定它就是真实模型；还需要结合残差、新数据上的预测表现以及不同模型设定下的结果，判断所得结论是否可靠。

备注

　　蒙特卡洛模拟和真实数据分析在目标上具有一致性，均用于揭示规律、预测结果或评估不确定性，且都依赖统计方法和概率理论。然而，两者的核心差异在于数据来源和应用逻辑。蒙特卡洛模拟基于 预先设定的真实模型生成人工数据⁠，通过大量随机抽样模拟复杂系统或极端场景，其优势在于灵活控制参数，但结果高度依赖真实模型的设定。相比之下，真实数据分析直接处理观测或实验数据，反映现实世界的真实分布；其对应的 真实模型往往较为复杂，且无法预先获知⁠。此外，真实数据受限于数据质量、噪声或获取成本，且可能无法覆盖所有潜在情况。蒙特卡洛模拟适用于理论验证、风险建模或数据稀缺时的假设测试，而真实数据分析更适用于实证研究或实际问题的直接分析。两者互补性强：模拟可辅助设计实验或验证方法，真实数据则能校准模型假设，提高模拟的可靠性。选择何种方法取决于研究目的、数据可得性及对现实逼近的需求。

训练集的生成

　　接下来，我们将编写用于产生训练集的 Python 函数，这里将用到 np.random.normal() 函数。函数输入样本量 n 和随机种子 rn，预期返回两个大小均为 \(n\times1\) 的数组：特征 x 与按真实模型生成的标签 y。

def train_data_generation(n, rn):
    # n: 样本量
    # rn: 随机种子

    np.random.seed(rn)  # 设置随机种子为 rn
    x = np.random.normal(2, 2, (n,1))   #从均值为2、标准差为2的正态分布中生成大小为nx1的特征向量
    epsilon = np.random.normal(0, 1, (n,1))  #从均值为0、标准差为1的正态分布中生成大小为nx1的随机误差
    y = 1 + x + epsilon   #通过(1)式生成y

    return x, y #打包输出相关结果

　　上面一段代码定义了一个名为 train_data_generation 的函数。该函数的输入由样本量 n 以及随机种子 rn 两部分组成，其输出为一个由 x 和 y 组成的元组。在函数内部，我们利用注释解释了函数的输入以及相应的代码。我们希望各位同学在将来自己编写函数时，也能养成用简单的注释解释代码的习惯。此外，在函数编写时，我们用的是有意义的字符串表示相应的结果。例如，我们用 epsilon 表示误差项，用 x 以及 y 表示特征和标签。我们用到的符号与 真实模型 中的符号相同。我们也希望同学们将来在编写代码时，也能够用具有意义的字符串，以便让代码具有较强的可读性。

分析模型

　　当分析数据时，我们往往不知道真实模型的形式。在这种情况下，我们往往借助数据可视化或者专家经验，给出分析模型的形式。当然，在蒙特卡洛模拟分析阶段，我们是知道真实模型的形式的。下面先生成包含 \(1\,000\) 个样本的 x 与 y，再用 matplotlib.pyplot 绘制散点图；横轴为特征、纵轴为标签，预期图像呈现带有随机波动的上升直线趋势。

import matplotlib.pyplot as plt  #我们通常将matplotlib.pyplot简单记为plt
x, y = train_data_generation(1000, 100)

plt.scatter(x,y)
plt.xlabel("Feature")
plt.ylabel("Label")
plt.show()

　　通过观察，我们发现标签 \(y\) 和特征 \(x\) 之间呈现线性关系。因此，我们选择如下分析模型对训练集进行拟合：

(2)\[\mathbb{E}(Y\mid X=x)=b+wx,\]

其中，\(b\) 和 \(w\) 为模型参数。

备注

　　建模前进行必要的数据可视化，有助于发现变量关系和异常观测。式 (2) 给出的分析模型与式 (1) 给出的真实模型形式一致；分析真实数据时无法直接观察真实模型，因此必须通过诊断和样本外评估检查模型设定是否合理。

备注

　　由于真实模型未知，线性、同方差或独立性等假设都可能失效。可通过残差图、样本外误差和针对性敏感性分析检查这些假设。神经网络通常比线性模型更灵活，但一个固定架构仍是有限维参数模型；灵活性不能自动消除测量误差、遗漏变量、分布漂移或因果解释问题。

损失函数

　　式 (2) 给出的分析模型由参数 \(b\) 和 \(w\) 确定。回归损失用于度量预测值与观测值之间的差异。针对线性回归等问题，下面的均方误差（mean squared error，MSE）是常见的损失函数，也常称为代价函数（cost function）：

(3)\[(\widehat{b},\widehat{w}) =\operatorname*{arg\,min}_{(b,w)} \frac{1}{n}\sum_{i=1}^n(y_i-b-wx_i)^2,\]

其中，\(\widehat{b}\) 和 \(\widehat{w}\) 是在 MSE 意义下的最优参数估计量。

备注

　　损失函数与代价函数。 按照一种常见且较严格的区分，损失函数（loss function）衡量模型对单个样本的预测误差，例如第 \(i\) 个样本的损失可以写成 \(\ell_i(\btheta)\)⁠；代价函数（cost function）则对训练集中所有样本的损失求和或取平均，并且可以加入正则化项，用于确定模型参数。不同教材和软件库并不总是严格区分这两个名称，常把整个训练目标也称为损失函数。因此，阅读公式或代码时，应注意它作用于单个样本还是一组样本，以及其中是否包含正则化项。式 (3) 对全部样本的平方误差取平均，按照上述区分属于代价函数，也常称为均方误差损失。

　　MSE 对较大误差施加二次惩罚，因此对异常值较敏感；均方根误差（RMSE）与观测标签的量纲一致。平均绝对误差（MAE）对异常值更稳健，但在误差为零处不可导。Huber 损失在小误差区域采用二次形式、在大误差区域采用线性形式，在光滑性与稳健性之间取得折中。损失函数应根据噪声分布、任务目标以及对高估和低估的不同代价来选择。

参数估计

　　记 \(\tilde{\bx}_i=(1,x_i)\trans\)⁠，\(\bX=[\tilde{\bx}_1\trans;\ldots;\tilde{\bx}_n\trans]\in\mathbb{R}^{n\times2}\)⁠，\(\by=(y_1,\ldots,y_n)\trans\)⁠。若 \(\bX\) 列满秩，由 正规方程 可知，最小化式 (3) 的参数估计量为

(4)\[\begin{split}\begin{pmatrix}\widehat{b}\\\widehat{w}\end{pmatrix} =\left(\sum_{i=1}^n\tilde{\bx}_i\cdot\tilde{\bx}_i\trans\right)^{-1}\cdot \sum_{i=1}^n\tilde{\bx}_iy_i =(\bX\trans\cdot\bX)^{-1}\cdot\bX\trans\cdot\by.\end{split}\]

　　在式 (1) 给出的模型满足条件均值假设且设计矩阵满列秩时，式 (4) 中的 \(\widehat{b}\) 与 \(\widehat{w}\) 均为无偏估计量。

　　第一个等式对应于基于循环运算的“求和法”，而第二个等式对应于向量化。作为本次的编程训练，我们将实现两种方法，并比较两种方法的计算效率。

基于 for 循环的求和法

　　下面按照式 (4) 的求和形式累积 \(\bX\trans\cdot\bX\) 和 \(\bX\trans\cdot\by\)⁠，再用 np.linalg.solve 解正规方程。函数输入按行存放的 \(n\times1\) 特征 x 和标签 y，输出长度为 2 的数组，依次给出截距与斜率估计。这样保留循环教学目的，同时避免显式计算逆矩阵。

def estimation_summation(x,y):
    # x: 长度为n的特征向量
    # y: 长度为n的标签向量

    n = len(x) #获取样本量n的大小
    aug_x = np.concatenate((np.ones_like(x),x), axis = 1) #将截距项1加入到特征中，得到增广特征矩阵，记为aug_x
    xx = np.zeros((2,2)) #用2X2的零矩阵初始化'xx'
    xy = np.zeros((2,1)) #用2X1的零矩阵初始化'xy'

    for i in range(n): #根据（3）的第一个方程计算'xx'和'xy'的求和
        example_i = aug_x[i,:].reshape((2,1))
        xx += example_i @ example_i.transpose()
        xy += example_i * y[i]

    par_est = np.linalg.solve(xx, xy)  # 解正规方程，不显式求逆

    return par_est.flatten()

　　aug_x 在特征左侧补上一列 1 以表示截距；循环逐个样本累加两个交叉乘积，最后一次性求解参数。返回值展平后形如 [b_hat, w_hat]。这种写法便于观察求和公式，但 Python 循环在大样本下通常较慢；若两列设计变量线性相关，solve 也无法得到唯一解。

向量化

　　为公平比较循环与向量化，下面两个版本都形成正规方程并使用同一个 solve 求解器，区别只在交叉乘积是循环累加还是矩阵乘法。下面的函数同样输入 \(n\times1\) 的 x、y，输出 [b_hat, w_hat]；随后再给出实际分析中更推荐的 lstsq 版本。

def estimation_vectorization(x,y):
    # x: 长度为n的特征向量
    # y: 长度为n的标签向量

    aug_x = np.concatenate((np.ones_like(x),x), axis = 1)  #将截距项1加入到特征中，得到增广特征矩阵，记为aug_x

    xx = aug_x.transpose() @ aug_x
    xy = aug_x.transpose() @ y
    par_est = np.linalg.solve(xx, xy)

    return par_est.flatten()

　　这里的两次矩阵乘法一次完成循环版中的全部累加，因此在相同输入上应得到相同的两个参数估计，只可能存在很小的浮点舍入差异。它通常比显式循环快，但正规方程会放大设计矩阵的条件数，数值稳定性不如直接最小二乘求解。

　　下面，我们将调用以上两个函数，验证两种方法的估计值是否一致。

print(estimation_summation(x,y))
print(estimation_vectorization(x,y))
print('以上两个结果应在浮点误差范围内一致')
print('两个参数的真实值为1和1。')

[1.01322026 0.99776934]
[1.01322026 0.99776934]
以上两个结果应在浮点误差范围内一致
两个参数的真实值为1和1。

　　前两行输出分别来自循环版和向量化版；二者相近说明两种程序实现了同一组正规方程，而估计值接近 [1, 1] 则与数据生成时的真实参数一致。这只是一次随机训练集上的核对，不能代替对多组样本和异常输入的测试。

　　为减少正规方程带来的数值误差，下面定义 estimation_stable。它接收同样的 x 与 y，直接对增广设计矩阵做最小二乘求解，返回截距和斜率；若设计矩阵不满列秩，则明确报错。

def estimation_stable(x, y):
    """使用 QR/SVD 类最小二乘求解，避免形成 X.T @ X。"""
    aug_x = np.concatenate((np.ones_like(x), x), axis=1)
    par_est, _, rank, singular_values = np.linalg.lstsq(
        aug_x, y, rcond=None
    )
    if rank < aug_x.shape[1]:
        raise np.linalg.LinAlgError("设计矩阵不满列秩")
    return par_est.flatten()

　　np.linalg.lstsq 通常通过 QR 或 SVD 类算法求解，不需要形成 \(\bX\trans\cdot\bX\)⁠；rank 用于确认两个参数能够被唯一识别。该函数比前两个版本更适合实际数值计算，但仍假定 x 与 y 为形状相容的二维列数组，且这里只返回点估计，不负责模型诊断。

备注

　　尽管两个方法的估计值相同，但其计算效率差异较大。我们将在 计算效率 一节中对两种方法进行比较。

备注

　　向量化在数据处理和算法实现中具有重要意义。通过将数据表示为向量或矩阵形式，向量化允许使用线性代数中的矩阵运算，这些运算通常 比传统的循环运算更快⁠，特别是在处理大型数据集时。这种表示方式 不仅提高了计算效率，还简化了代码实现⁠，因为可以使用统一的数学运算来处理数据，减少了需要编写的代码量。此外，向量化还便于利用现代计算机体系结构的并行计算能力，如 GPU，从而进一步加速数据处理和算法执行。因此，向量化是数据科学和机器学习领域中一种基础且关键的数据处理手段。我们在接下来的课程中，将详细介绍该方法在不同深度学习框架下的应用。

　　在传统的统计模型分析中，统计推断是一个重要的内容。但该部分内容涉及到较多的数理统计相关内容，基础较为薄弱的同学可以跳过下面的内容。

统计推断

　　通常，我们会进行 1,000 次蒙特卡洛模拟：根据式 (1) 给出的真实模型，独立生成 1,000 个大小为 \(n\) 的训练集，并得到参数估计 \((\widehat{b},\widehat{w})\)⁠。汇总这些估计即可研究偏差、方差、均方误差与置信区间覆盖率。循环版和向量化版实现的是同一个估计量，后续只使用更稳健且高效的向量化版本。

　　我们首先产生 1,000 个训练集，并将对应的参数估计值保存在一个 \(1000\times2\) 的矩阵中，其中每一行代表基于对应训练集得到的参数估计值。

n = 1000 #样本量
par_result = np.zeros((1000, 2))  # 保存 1000 次参数估计

for i in range(1000):
    x, y = train_data_generation(n=n, rn=i) # the random seed is set to be i
    par_result[i,:] = estimation_stable(x, y)

　　循环用不同种子生成 \(1\,000\) 份独立训练集，并把每次得到的 [b_hat, w_hat] 存入 par_result 的一行。因此最终数组大小为 \(1000\times2\)⁠；它描述估计量在重复抽样下的变化，而不是一份训练集中的 \(1\,000\) 个模型参数。

　　接下来我们计算 1,000 次蒙特卡洛模拟估计值的偏差、方差以及均方误差。此处，我们将用到上节所复习的 np.mean() 和 np.var() 两个函数。

bias_est = np.mean(par_result,axis = 0)-1 #由于参数真值均为1，故我们减掉1，以得到估计量的蒙特卡洛偏差
var_est = np.var(par_result,axis = 0)
mse_est = bias_est**2 + var_est
print('参数b_0和w_0的估计偏差为：')
print(bias_est)
print("*"*6)
print('参数b_0和w_0的估计方差为：')
print(var_est)
print("*"*6)
print('参数b_0和w_0的MSE为：')
print(mse_est)

参数b_0和w_0的估计偏差为：
[5.59456735e-05 2.83427856e-04]
******
参数b_0和w_0的估计方差为：
[0.00200172 0.0002554 ]
******
参数b_0和w_0的MSE为：
[0.00200173 0.00025548]

　　三组长度为 2 的输出依次对应截距和斜率。偏差衡量重复估计的平均值与真值 1 的差，方差衡量估计值的波动，bias_est**2 + var_est 给出按当前总体方差约定计算的蒙特卡洛 MSE。结果只反映本节指定的数据生成模型、样本量和 \(1\,000\) 次重复，不能直接推广到其他噪声或特征分布。

　　在式 (1) 给出的模型设定下，给定设计矩阵后，估计量 \((\widehat{b},\widehat{w})\) 的条件协方差为

(5)\[\begin{split}\operatorname{Cov}\!\left[ \begin{pmatrix}\widehat{b}\\\widehat{w}\end{pmatrix} \mathrel{\bigg|}\bX\right] =\sigma^2(\bX\trans\cdot\bX)^{-1},\end{split}\]

其中，\(\sigma^2=1\) 是误差项方差。在常用的正则条件下，\(n^{-1}\bX\trans\cdot\bX\) 收敛到正定矩阵，因此估计量的协方差矩阵与 \(n^{-1}\) 同阶。实际推断时，可用残差方差估计 \(\sigma^2\)⁠。令

\[\widehat{\epsilon}_i=y_i-\widehat{b}-\widehat{w}x_i.\]

我们利用 \(\{\widehat{\epsilon}_i:i=1,\ldots,n\}\) 的残差平方和除以自由度 \(n-2\)⁠，得到 \(\widehat{\sigma}^2\)⁠。即 \(\widehat{\sigma}^2 = (n-2)^{-1}\sum_{i=1}^n\widehat{\epsilon}_i^2\)⁠。

　　下面的 se_est 输入一份 \(n\times1\) 的 x 与 y，先拟合截距和斜率，再返回长度为 2 的标准误数组，分别量化两个参数估计在该样本下的不确定性。

def se_est(x,y):
    # x: 长度为n的特征向量
    # y: 长度为n的标签向量

    est_par = estimation_stable(x, y)  # 获取 (b_0, w_0) 的估计量
    error_est = y - est_par[0] - est_par[1] * x  #获取误差项
    dof = len(x) - 2  # 两个待估参数：截距与斜率
    if dof <= 0:
        raise ValueError("至少需要 3 个样本来估计残差方差")
    sigma2_est = np.sum(error_est**2) / dof

    aug_x = np.concatenate((np.ones_like(x),x), axis = 1) #用截距项增广特征向量
    xx = aug_x.transpose() @ aug_x  #通过 X^TX 获取 'xx'，其中 X 是大小为 nX2 的增广矩阵
    xx_inv = np.linalg.solve(xx, np.eye(xx.shape[0]))
    sd_par = np.sqrt(sigma2_est * np.diag(xx_inv))

    return sd_par

　　程序用残差平方和除以 \(n-2\) 估计误差方差，再从 \(\widehat{\sigma}^2(\bX\trans\bX)^{-1}\) 的对角线取得两个方差并开平方。它要求至少 3 个样本且增广设计矩阵满列秩；标准误公式还依赖线性模型、独立同方差误差等假设。

　　接下来，将函数给出的标准误估计与 \(1\,000\) 次蒙特卡洛模拟中估计量的经验标准差进行比较，以检查程序是否正确。

n = 1000
sd_result = np.zeros((1000,2))

for i in range(1000):
    x, y = train_data_generation(n=n, rn=i) # the random seed is set to be i
    sd_result[i,:] = se_est(x,y)
print('函数计算结果为:')
print(np.mean(sd_result,axis = 0))
print('蒙特卡洛模拟结果为：')
print(np.std(par_result, axis=0, ddof=1))

函数计算结果为:
[0.04477234 0.01583325]
蒙特卡洛模拟结果为：
[0.04476301 0.01598918]

　　第一组输出是 \(1\,000\) 份样本内标准误的平均，第二组是 \(1\,000\) 个参数估计的经验标准差；两者接近表明标准误公式与重复抽样波动相符。它们不必逐项完全相等，因为都含有有限次数模拟造成的随机误差。

　　如果我们对方差的估计结果非常接近于蒙特卡洛模拟结果，这就在某种程度上说明了我们所编写的方差估计函数没有问题。这种利用理论结果和蒙特卡洛模拟结果检验程序的方法，在统计分析中较为常用。基于此，我们可以进一步构造双边置信区间，并检验所置信区间的覆盖率。这一步分析对于传统统计模型来说较为重要。我们可以证明，在式 (1) 给出的模型设定下，当样本量 \(n\) 趋于无穷大时，以下渐近结果成立：

\[\begin{split}\widehat{\sigma}^{-1}(\bX\trans\cdot\bX)^{1/2}\cdot \left[ \begin{pmatrix}\widehat{b}\\\widehat{w}\end{pmatrix} -\begin{pmatrix}b_0\\w_0\end{pmatrix} \right] \xrightarrow{d}\mathcal{N}(\bzero,\bI_2),\end{split}\]

其中，\((\bX\trans\cdot\bX)^{1/2}\) 表示矩阵“平方根”，\(\bI_2\) 为 \(2\times2\) 单位矩阵。由此可分别构造 \(b_0\) 和 \(w_0\) 的双侧 \(95\%\) 置信区间：

\[\begin{split}\begin{aligned} &\left(\widehat{b}-q_{0.975}\widehat{\sigma}_b, \widehat{b}+q_{0.975}\widehat{\sigma}_b\right),\\ &\left(\widehat{w}-q_{0.975}\widehat{\sigma}_w, \widehat{w}+q_{0.975}\widehat{\sigma}_w\right), \end{aligned}\end{split}\]

其中，\(q_{0.975}\) 为标准正态分布的 \(97.5\%\) 分位数，\(\widehat{\sigma}_b\) 和 \(\widehat{\sigma}_w\) 分别为两个估计量的标准误估计。对每个模拟样本检查置信区间是否覆盖真实参数，即可估计覆盖率。

　　下面的 cr_indicator 输入一份训练集和显著性水平 alpha，为截距、斜率各构造一个正态近似双侧区间，并返回两个布尔值，表示相应区间是否覆盖真实参数 1。

def cr_indicator(x,y,alpha):
    # x: 长度为n的特征向量
    # y: 长度为n的标签向量
    # alpha: 显著性水平。例如，alpha=0.05对应95%的双侧置信区间。

    quan_normal = scipy.stats.norm.ppf(1-alpha/2) #获取标准正态分布的第 (1-alpha/2) 分位数。
    est_par = estimation_stable(x, y)  # 获取参数估计值
    est_se = se_est(x,y) #获取标准误估计值。

    # 检查 b_0 与 w_0 的真实值 1 是否落在相应置信区间内
    ind_b = est_par[0] - quan_normal * est_se[0] < 1 < est_par[0] + quan_normal * est_se[0]
    ind_w = est_par[1] - quan_normal * est_se[1] < 1 < est_par[1] + quan_normal * est_se[1]

    return np.array([ind_b, ind_w])

　　ppf 给出区间所需的临界值，随后程序用“估计值 \(\pm\) 临界值 \(\times\) 标准误”形成区间并检查真值。返回的 [ind_b, ind_w] 只记录是否覆盖，不返回区间端点；该实现还把真值 1 写死，只适用于本节模拟设定。

　　下面，我们计算对应置信水平下的覆盖率。

import scipy.stats # 为了得到标准正态分布的分位数
cr_est = np.zeros_like(par_result, dtype=bool)
for i in range(1000):
    x, y = train_data_generation(n=n, rn=i) # the random seed is set to be i
    cr_est[i,:]  = cr_indicator(x, y, alpha = 0.05)
print('在置信水平为0.95时，双侧置信区间的覆盖率为:')
print(np.mean(cr_est, axis = 0))

在置信水平为0.95时，双侧置信区间的覆盖率为:
[0.955 0.952]

　　cr_est 的每行对应一次模拟、两列对应两个参数；对布尔值按列取平均，就是两个 \(95\%\) 区间的经验覆盖率。有限的 \(1\,000\) 次模拟会带来蒙特卡洛误差，因此输出无需恰好等于 0.95。

　　若估计覆盖率在蒙特卡洛误差范围内接近 0.95，说明结果符合渐近理论的相关结论，但不能单凭这一项检查证明整个程序正确。此外，小样本下若依赖高斯线性模型的精确推断，因此在统计推断时，使用 \(t\) 分位数其实更为稳妥。

计算效率

　　尽管基于 for 循环的求和法与向量化方法在数值误差范围内给出相同估计，但求和法通常较慢。下面以 \(n=1\,000\) 为例比较运行时间。针对深度学习模型，正式比较不同模型程序的运行时间时，应先运行几次，使程序完成必要的初始化；随后重复计时多次，并报告运行时间的中位数以及四分位距，以说明程序通常的运行速度及其波动情况。

import time # 用来记录算法的计算时间
T1 = time.time()
for i in range(1000):
    x, y = train_data_generation(n=n, rn=i) # the random seed is set to be i
    par_est = estimation_summation(x,y)
T2 = time.time()
print('基于 for 循环的求和法耗时约 %4.4f 秒' % (T2 - T1))

T1 = time.time()
for i in range(1000):
    x, y = train_data_generation(n=n, rn=i) # the random seed is set to be i
    par_est = estimation_vectorization(x,y)
T2 = time.time()
print('向量化方法耗时约 %4.4f 秒' % (T2 - T1))

基于 for 循环的求和法耗时约 1.3826 秒
向量化方法耗时约 0.0307 秒

　　两行输出是在相同的 \(1\,000\) 份数据规模下完成相同参数估计所用的总时间；通常向量化版本更快。这里每段代码只计时一次，而且计时还包含随机数据生成，结果会受机器负载和软件环境影响；正式基准应把待比较部分隔离，并按上文所述预热、重复计时和报告波动范围。

核心推导与实现核验

核心关系

\[\widehat{\btheta}=(\bX\trans\cdot\bX)^{-1}\cdot\bX\trans\cdot\by.\]

　　推导路径。 对平方损失求梯度得 \(2\bX\trans\cdot(\bX\cdot\btheta-\by)\)⁠，令其为 \(\bzero\) 得到正规方程；满列秩时可写出解析解。数值实现应使用 solve 或 lstsq，避免显式求逆。

数值梯度检验

　　对参数向量 \(\btheta=(\theta_1,\ldots,\theta_p)\trans\)⁠，梯度 \(\nabla_{\btheta}\mathcal{J}(\btheta)\) 由损失函数对各个参数的一阶偏导数组成。公式推导或后向传播程序会给出一个梯度，下面将其称为 解析梯度⁠。为了检查这个梯度是否计算正确，可以在不使用原推导过程的情况下，直接观察参数发生很小变化时损失函数怎样变化，由此得到 数值梯度⁠，再比较两种结果。

　　记 \(\be_j\) 为第 \(j\) 个分量等于 1、其余分量等于 0 的向量。选择一个很小的正数 \(h\)⁠，只改变第 \(j\) 个参数而保持其他参数不变，则

\[g_{\mathrm{num},j} =\frac{\mathcal{J}(\btheta+h\be_j) -\mathcal{J}(\btheta-h\be_j)}{2h} \approx \frac{\partial\mathcal{J}(\btheta)}{\partial\theta_j},\]

其中，\(g_{\mathrm{num},j}\) 是第 \(j\) 个参数对应的数值梯度。这个方法称为 中心差分 （central difference）：它分别在当前参数值的两侧计算一次损失，再用两次损失之差除以两侧参数值之差 \(2h\)⁠，从而近似当前位置的变化率。对所有参数逐个进行这项计算，便得到数值梯度向量 \(\bg_{\mathrm{num}}\)⁠。

　　可以按照下面的步骤检验梯度：

1. 选取一个规模很小、可以重复计算的例子，固定数据和参数，并使用 float64 保存数值；

2. 通过公式推导或后向传播程序计算解析梯度 \(\bg\)⁠；

3. 每次只把一个参数增加 \(h\) 或减少 \(h\)⁠，分别重新计算损失，并用中心差分得到 \(\bg_{\mathrm{num}}\)⁠；

4. 比较两个梯度的每个分量，同时计算最大绝对误差和整体相对误差，例如

   \[E_{\mathrm{abs}} =\max_j\left|g_j-g_{\mathrm{num},j}\right|, \qquad E_{\mathrm{rel}} =\frac{\lVert\bg-\bg_{\mathrm{num}}\rVert_2} {\max\!\left\{10^{-12}, \lVert\bg\rVert_2+ \lVert\bg_{\mathrm{num}}\rVert_2\right\}}.\]

5. 尝试几个不同的 \(h\)⁠，例如 \(10^{-4}\)⁠、\(10^{-5}\) 和 \(10^{-6}\)⁠。在光滑的小型 float64 例子中，相对误差达到 \(10^{-6}\) 或更小通常能够为梯度实现提供较强的正确性证据，但这个数值只是参考，不能脱离函数、数据尺度和计算精度机械套用。

　　下面的程序展示了上述步骤。loss_fn 接收参数向量并返回一个标量损失；gradient_fn 表示需要接受检验的梯度程序。计算数值梯度时，每次改变一个参数，计算完成后立即恢复该参数的原值。

import numpy as np

def numerical_gradient(loss_fn, theta, h=1e-5):
    theta = np.asarray(theta, dtype=np.float64).copy()
    grad = np.zeros_like(theta)

    for j in range(theta.size):
        original_value = theta[j]

        theta[j] = original_value + h
        loss_plus = loss_fn(theta)

        theta[j] = original_value - h
        loss_minus = loss_fn(theta)

        theta[j] = original_value
        grad[j] = (loss_plus - loss_minus) / (2 * h)

    return grad

def check_gradient(loss_fn, gradient_fn, theta, h=1e-5):
    analytic_grad = np.asarray(
        gradient_fn(theta), dtype=np.float64
    )
    numeric_grad = numerical_gradient(loss_fn, theta, h=h)

    max_abs_error = np.max(
        np.abs(analytic_grad - numeric_grad)
    )
    relative_error = np.linalg.norm(
        analytic_grad - numeric_grad
    ) / max(
        1e-12,
        np.linalg.norm(analytic_grad)
        + np.linalg.norm(numeric_grad),
    )
    return max_abs_error, relative_error

# 示例：检验平方损失的梯度
X = np.array([[1.0, -1.0], [1.0, 0.5], [1.0, 2.0]])
y = np.array([0.0, 1.0, 3.0])
theta = np.array([0.2, 0.8])

def loss_fn(theta):
    error = X @ theta - y
    return np.mean(error**2)

def gradient_fn(theta):
    return 2 * X.T @ (X @ theta - y) / y.size

max_abs_error, relative_error = check_gradient(
    loss_fn, gradient_fn, theta
)
print("最大绝对误差：", max_abs_error)
print("相对误差：", relative_error)

　　如果两种梯度的符号、大小或个别分量明显不同，应依次检查损失函数是否完全相同、求导符号是否写反、对样本求和与取平均是否混用，以及数组维度或广播是否导致了错误计算。修改程序后，应更换几组数据和参数并重复检验，不能只让一个例子通过。

自动微分对照

　　如果模型由 PyTorch、JAX 等框架支持的运算构成，还可以让框架通过 自动微分 （automatic differentiation）计算梯度，再与手工推导或自行编写的后向传播结果逐项比较。对照时必须使用完全相同的数据、参数、损失定义以及求和或取平均方式，并关闭 Dropout 等随机变化；否则，梯度差异可能来自计算设置不同，而不一定来自求导错误。自动微分能够提供第二种实现作为参照，但它仍然建立在同一模型和损失定义上，不能检查模型定义本身是否符合研究问题。

手算结果与可信实现对照

　　对规模很小的例子，可以先手工计算中间量和最终结果，再与程序逐项比较；若已有经过广泛使用的程序库，还可以让自行实现和可信库处理完全相同的输入，并比较输出与梯度。使用这种方法时，应先统一数据类型、参数、数组轴、边界处理、是否取平均等计算规则，并设置与数值精度相符的误差范围。手算或可信库对照能够发现公式翻译、数组维度和边界处理中的错误，但“两个程序结果一致”也不能单独证明它们采用的计算规则一定适合当前任务。

使用数值梯度检验时需要注意

　　\(h\) 过大时，中心差分不能准确描述当前位置附近的变化；\(h\) 过小时，两次非常接近的损失相减会放大计算机舍入误差。因此，应尝试多个 \(h\)⁠，而不是只看一次结果。还要固定程序中的随机变化，并避开损失函数不可导的位置。数值梯度需要为每个参数至少计算两次损失，计算量很大，所以它适合检查小型例子或抽查少量参数，不适合代替训练时的后向传播。

关键条件

　　在线性模型条件均值正确、设计矩阵满列秩且误差满足 \(\mathbb{E}(\bepsilon\mid\bX)=\bzero\) 时，给定 \(\bX\) 有 \(\mathbb{E}(\widehat{\btheta}\mid\bX)=\btheta_0\)⁠，再用重期望即可得到无偏性。

数据规模

　　若有 \(n\) 个样本、每个样本有 \(p\) 个特征，加入一列常数作为截距后，设计矩阵 \(\bX\) 有 \(n\) 行、\(p+1\) 列；标签向量 \(\by\) 有 \(n\) 个数，待估计参数 \(\widehat{\btheta}\) 有 \(p+1\) 个数。

常见误区

　　把随机种子设在函数内部会使每次调用产生完全相同的数据；训练集与测试集若复用同一随机数，或者在计算标准化所需的均值和标准差、缺失值填充值等信息时使用了测试集，也会使训练过程提前获得测试集的信息，造成数据依赖或数据泄漏。

动手检查

　　验证正规方程残差 X.T @ (X @ w_hat - y) 接近 0，并与 np.linalg.lstsq 对照；重复模拟再检查估计偏差随重复数稳定。

数值稳定性与规模

　　最小二乘不要显式计算逆矩阵；使用 np.linalg.lstsq 或 QR/SVD。共线性严重时同时报告条件数，并考虑中心化、标准化或正则化。

本节小结

1. 随机种子保证可复现，但不应混淆独立重复。

2. 正规方程来自平方损失的一阶条件。

3. 中心差分 可以在小型例子中检验解析梯度或后向传播程序。

4. 向量化、稳定线性代数等计算设计与对应模型的理论正确性同等重要。

综合练习

　　题目依次覆盖记号、概念、推导、证明、维度、实现、数值核验和受控实验。程序题应固定随机种子、写出维度断言并报告运行环境。全部参考答案见 Python 线性回归实验与向量化答案⁠。

1. 公式推导。 从正文定义的损失函数出发，推导 \(\widehat{\btheta}=(\bX\trans\cdot\bX)^{-1}\cdot\bX\trans\cdot\by\)⁠。

2. 证明题。 在线性模型条件均值正确且设计矩阵满列秩时，证明最小二乘估计量条件无偏。写清假设、关键等式以及结论的适用范围。

3. 维度检查。 针对“Python 线性回归实验与向量化”，按正文的批量约定写出核心关系中输入、参数、中间量和输出的维度；逐项验证乘法、求和、转置或广播是否合理。

4. 核心编程。 实现稳定的最小二乘估计、重复模拟及计时函数，并返回带标准误的汇总表。函数应检查输入、避免不必要的隐式广播，并返回便于核验的中间量。

5. 数值稳定性。 针对“最小二乘估计”及 \(\widehat{\btheta}=(\bX\trans\cdot\bX)^{-1}\cdot\bX\trans\cdot\by\) 检查真正可能出现的溢出、下溢、除零、病态或消减问题；不要机械罗列与本节无关的风险，并给出稳定改写。

6. 受控实验。 改变样本量、噪声方差和共线性程度，比较估计偏差、方差、覆盖率与运行时间。除研究因素外保持数据划分、随机种子集合、训练预算和评价代码一致。

7. 可视化。 为“Python 线性回归实验与向量化”的受控实验选择能直接检验假设的横轴、纵轴和分组变量；显示均值与波动范围，并写出一个图中不能支持的因果结论。

8. 综合任务。 以 \(\widehat{\btheta}=(\bX\trans\cdot\bX)^{-1}\cdot\bX\trans\cdot\by\) 为理论核心，实现“实现稳定的最小二乘估计、重复模拟及计时函数，并返回带标准误的汇总表”，再完成“改变样本量、噪声方差和共线性程度，比较估计偏差、方差、覆盖率与运行时间”；写清数据、伪代码、评价指标和结论边界。

Make live