线性模型

线性模型（linear model）用若干已知量的线性组合表示输出。在统计建模中，“线性”通常指对未知参数线性：输入可以经过平方、对数等已知变换，曲线也未必是一条直线。本篇以线性回归为主，区分模型、拟合方法与由模型得到的推断；一般线性方程组的结构见线性代数。

四次水位读数，三个不同的对象

沿用最小二乘法中的合成教学数据：在 $t=0,1,2,3$ 分钟读取水位，得到 $h=1,2,2,4$ 厘米。想回答的是：在这个时间范围内，平均水位怎样随时间变化？若再观测一次，其读数大约在哪里？

先提出模型 $$H_i=b+a{t}_i+{\epsilon }_i.$$ 这里 $H_i$ 是第 $i$ 次观测的随机水位，$h_i$ 是已经拿到的数值；$b$ 是未知截距，单位厘米；$a$ 是未知斜率，单位厘米／分钟。误差 ${\epsilon }_i$ 表示模型直线之外的变化，也用厘米计。

假设在给定观测时刻后误差均值为零，即 $E[{\epsilon }_i\mid t_i]=0$，那么 $$E[H_i\mid t_i]=b+a{t}_i.$$ 这是一条关于平均响应的直线，不要求每次读数都在线上。数据拟合后得到的 $\hat{h}(t)=0.9+0.9t$ 则是这条未知平均直线的一次估计。模型中的真参数、由样本算出的估计值、下一次带随机误差的读数，是三个不同对象。Penn State STAT 501：Simple Linear Regression

实心点确定拟合线，空心点用于后续检验。模型直线描述平均水位，观测点还包含误差。

“对参数线性”为什么允许弯曲

把模型写成 $$H_i={\beta }_0{\varphi }_0(t_i)+\cdots +{\beta }_{p-1}{\varphi }_{p-1}(t_i)+{\epsilon }_i,$$ 其中 ${\varphi }_j$ 是事先选择的已知函数，${\beta }_j$ 是待估计的系数。每个系数只乘一个已知数，再相加，这就是这里的线性。

例如 ${\varphi }_0(t)=1,{\varphi }_1(t)=t,{\varphi }_2(t)=t^2$ 给出 $$E[H\mid t]={\beta }_0+{\beta }_{1}t+{\beta }_2{t}^2.$$ 它可以是一条抛物线，却仍是线性参数模型。三种系数的单位分别是厘米、厘米／分钟、厘米／分钟平方；每项相加前具有同样的水位单位。若模型含 $e^{-\beta t}$，未知参数进入指数，便不再是上述线性参数形式。对参数线性与对输入线性需要区分，NIST 的线性回归章节也采用这一约定。NIST：Linear Least Squares Regression

这也解释了与线性代数用语的差别：$t\mapsto b+at$ 在 $b\ne 0$ 时是仿射映射；把 $(b,a)$ 看作输入，而把全部预测看作输出时，这个参数到预测的映射是线性的。

从一行观测组成设计矩阵

每次观测写一行，每个模型项写一列。直线模型的设计矩阵为 $$X=\left(\begin{array}{cc}1& 0\\ 1& 1\\ 1& 2\\ 1& 3\end{array}\right),\beta =\left(\begin{array}{c}b\\ a\end{array}\right),H=X\beta +\epsilon .$$ 四行是四次观测，第一列的四个一让每次都含同一个截距。若加入平方项，只需另加列 $(0,1,4,9{)}^{𝖳}$。矩阵尺寸为 $n\times p$，$n$ 是观测数，$p$ 是包括截距在内的参数数目。

设计矩阵还告诉我们参数能否区分。若两列完全相同，改变一列的系数、反向改变另一列系数，预测不会变化。即使数据很多，也无法从这种设计中识别两个独立贡献。参数唯一可识别要求 $\mathrm{rank}X=p$，即列满秩；问题在于独立信息，不只是行数多。

最小二乘怎样估计这个模型

模型选好之后，还要规定如何估计参数。普通最小二乘法选择 $$\hat{\beta }={\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}_{\beta }\sum _{i=1}^n(h_i-(X\beta {)}_i{)}^2.$$ 残差是观测减拟合值。令目标对每个系数的导数为零，得到 $$X^{𝖳}X\hat{\beta }=X^{𝖳}h.$$ 对当前数据，方程是 $$4\hat{b}+6\hat{a}=9,6\hat{b}+14\hat{a}=18.$$ 第二式减去第一式的 $3/2$ 倍，得到 $5\hat{a}=4.5$，故 $\hat{a}=0.9$；代回得到 $\hat{b}=0.9$。这一步消去了两式中的 $6\hat{b}$，只留下斜率。

四个拟合值为 $0.9,1.8,2.7,3.6$，残差为 $0.1,0.2,-0.7,0.4$，平方和为 $0.7$ 平方厘米。目标的最小性、投影解释与稳定的 QR 解法见最小二乘法；实际计算通常不显式求逆矩阵。

这一步只是一项确定的优化计算，还没有要求误差正态。把“能算出一条拟合线”进一步解释成“参数估计无偏”或“区间有指定覆盖率”，才需要下一节的统计假设。

哪些假设支持哪些结论

以下把整个设计矩阵 $X$ 固定下来，或者条件在给定的 $X$ 上。若真实平均响应为 $X\beta$，并且 $E[\epsilon \mid X]=0$，列满秩时 $$\hat{\beta }=\beta +(X^{𝖳}X{)}^{-1}{X}^{𝖳}\epsilon .$$ 这是把 $H=X\beta +\epsilon$ 代入估计公式后直接得到的。取条件期望，第二项为零，所以 $E[\hat{\beta }\mid X]=\beta$：重复按同一设计采样，估计量的平均值等于真参数。

若还假设误差的条件协方差矩阵为 ${\sigma }^2{I}_n$，即每次误差方差相同、不同误差互不相关，则 $$\mathrm{Cov}(\hat{\beta }\mid X)={\sigma }^2(X^{𝖳}X{)}^{-1}.$$ 原因是对矩阵 $A=(X^{𝖳}X{)}^{-1}{X}^{𝖳}$，有 $\mathrm{Cov}(A\epsilon )=A({\sigma }^2{I}_n)A^{𝖳}$；相乘后中间的 $X^{𝖳}X$ 抵消。

互不相关比独立弱；上述两个矩的结论不需要正态性。若进一步假定误差独立且服从 $N(0,{\sigma }^2)$，则能得到小样本下精确的 t 区间。若只满足零均值而方差不齐，点估计仍可能无偏，但刚才那条协方差公式要改变。漏掉一个与输入有关的系统性变化时，零均值假设本身也可能失败。

在 $n>p$ 和上述同方差、互不相关条件下，常用 $$s^2=\frac{\sum _i{r}_i^2}{n-p}$$ 估计误差方差。拟合用掉了 $p$ 个独立方向，残差所在的正交补只剩 $n-p$ 维；每个正交误差方向期望贡献 ${\sigma }^2$，故残差平方和的期望为 $(n-p){\sigma }^2$。本例 $n=4,p=2$，所以 $s^2=0.7/2=0.35$。它与训练均方误差 $0.7/4$ 的用途不同。

平均水位的区间与下一次读数的区间

选择训练范围内的时刻 $t_0=2.5$ 分钟，拟合水位为 ${\hat{h}}_0=0.9+0.9\times 2.5=3.15$ 厘米。先在独立正态、同方差且平均函数正确的模型下进行计算。

对带截距的一元直线模型，令 $$q(t_0)=\frac{1}{n}+\frac{(t_0-\overline{t}{)}^2}{\sum _i(t_i-\overline{t}{)}^2}.$$ 它来自把预测 $(1,t_0)\hat{\beta }$ 代入上一节的系数协方差。当前 $\overline{t}=1.5$、离差平方和为5，所以 $q(2.5)=1/4+1/5=0.45$。平均响应估计的标准误为 $${\mathrm{S}\mathrm{E}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}=s\sqrt{q(t_0)}=\sqrt{0.35\times 0.45}\approx 0.3969\mathrm{c}\mathrm{m}.$$ 预测一个独立的新读数，还要加入该次新误差的方差 ${\sigma }^2$，所以预测误差的估计标准差为 $${\mathrm{S}\mathrm{E}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{d}}=s\sqrt{1+q(t_0)}=\sqrt{0.35\times 1.45}\approx 0.7124\mathrm{c}\mathrm{m}.$$ 平方根里多出的1，正是新观测本身的随机变化，而不是另一次参数拟合。

为什么这里使用 t 分布？记真实平均响应为 ${\mu }_0=b+a{t}_0$。在上述正态模型下，拟合均值的标准化误差是标准正态变量；它与残差方差估计独立，且 $(n-p)s^2/{\sigma }^2$ 服从自由度为 $n-p$ 的卡方分布。因此 $$\frac{{\hat{h}}_0-{\mu }_0}{s\sqrt{q(t_0)}}\sim t_{n-p}.$$ 独立的新误差再加进来时，分子方差变成 ${\sigma }^2[1+q(t_0)]$，同样得到用于预测的 t 比值。这与学生t分布中“正态误差除以独立估计的标准差”的构造一致。

自由度为 $n-p=2$，双侧95%区间使用学生t分布分位数 $t_{0.975,2}\approx 4.30265$。分别乘上两个标准误： $$\text{平均响应：}3.15\pm 4.30265\times 0.3969,$$ $$\text{新读数：}3.15\pm 4.30265\times \mathrm{0.7124.}$$ 用未舍入值计算，前者约为 $[1.4424,4.8576]$ 厘米，后者约为 $[0.0848,6.2152]$ 厘米。四次观测只剩两个残差自由度，区间很宽，这也体现了小样本能提供的信息有限。

两条区间的中心相同，估计对象不同。区间宽度均按未舍入的标准误和t分位数计算。

两种95%都描述在所列模型下重复采样的覆盖率：一个覆盖固定时刻的真实平均响应，另一个覆盖该时刻新的随机读数。不是整条曲线同时具有95%的覆盖率，也不是所有未来读数必然落在区间内。预测区间的额外方差与覆盖解释见 NIST：预测单次新响应。

残差怎样指出直线遗漏的结构

另取一组明确构造的数据：$x=-3,-2,-1,0,1,2,3$，$y=x^2$。这里没有随机噪声，目的是检查模型形状。若坚持拟合直线，由对称性 $\overline{x}=0$、$\sum x_i{y}_i=\sum x_i^3=0$，斜率为零；截距是 $$\overline{y}=\frac{9+4+1+0+1+4+9}{7}=4.$$ 残差依次为 $5,0,-3,-4,-3,0,5$。它们的和为零，却清楚地在两端偏正、中间偏负。

残差总和为零是含截距最小二乘的代数性质，不能据此认定模型形状正确。上下图共用横轴位置。

加入平方项后，系数 $({\beta }_0,{\beta }_1,{\beta }_2)=(0,0,1)$ 恰好表示全部数据。这个例子同时说明：有系统形状的残差提示遗漏项；加入曲线项仍可能属于线性参数模型。真实样本有随机波动，选择更复杂模型还需要样本外验证，不能只因训练误差更小就不断加参数。NIST：模型验证与残差分析

拟合优度、外推与解释

含截距的普通最小二乘满足总离差平方和等于回归解释部分加残差平方和，因此可定义 $$R^2=1-\frac{\sum _i{r}_i^2}{\sum _i(h_i-\overline{h}{)}^2}.$$ 分母非零时才有这个表达。本例分母为4.75，所以 $R^2=1-0.7/4.75=81/95\approx 0.8526$。它描述这组训练资料中的平方和比例，不是预测正确率，也不表示斜率具有因果含义。无截距拟合或样本外评价也不自动保有训练时 $0\le R^2\le 1$ 的性质。

从零到三分钟的资料推出二点五分钟，是在已有范围内估计；推到一百分钟则是外推，设备容量、进出水规律等都可能变化。此前留出的四、五分钟两点只能示范检验流程，不能保证更远的预测。若要把“时间增加一单位”解释成干预效果，还需控制混杂或说明实验设计，回归公式本身不给出因果保证。

线性模型与差分方程模型也回答不同问题：前者可描述某时刻的平均响应与输入关系，后者明确规定状态怎样从一步更新到下一步。若把滞后状态用作回归输入，还必须处理时间依赖与误差条件，不能将独立样本公式直接照搬。

参考资料与后续

• NIST/SEMATECH：Linear Least Squares Regression，线性参数模型与形式。
• Penn State STAT 501，Lesson 1，总体回归、误差条件、方差估计与拟合优度。
• NIST：预测新响应及其不确定性。
• NIST：How can I tell if a model fits my data?，残差与模型检验。
• Boyd、Vandenberghe：Introduction to Applied Linear Algebra，设计矩阵、最小二乘与模型拟合。
• 先修：线性代数、最小二乘法、期望；后续：正态分布、学生t分布、统计推断。