5.2 分段多项式和样条函数

从本节至第 5.7 节，我们都假设 $X$ 是一维变量。将 $X$ 的定义域划分为相接的（contiguous）区间，在每个区间中用单独的多项式来代表 $f$，则得到了一个分段多项式函数 $f(X)$。图 5.1 展示了两个简单的分段多项式函数。

左上图为分段常数，它的三个基函数为：

由于它们在三个不相交的区域上取值为一，所以模型 $f(X)=\sum_{m=1}^3\beta_mh_m(X)$ 的最小二乘估计为 $\hat{\beta}(m)=\bar{Y}_m$，即在区域 $m$ 上 $Y$ 的均值。

右上图是一个分段线性的拟合。模型需要增加三个基函数：

除非是特殊的场景，左下图通常是更好的选择，它同样是分段线性的，但限制了函数在两个节点处的连续性。连续性的限制相当于对参数的线性约束条件；例如限制 $f(\xi_1^-)=f(\xi_1^+)$ 相当于 $\beta_1+\xi_1\beta_4=\beta_2+\xi_1\beta_5$。图中有两个节点便有两个限制，所以有 2 个参数可从线性约束中推导出，即模型整体有 4 个可自由变化的参数。

一个更直接的方法是使用已经满足了约束条件的基函数：

其中 $t_+$ 的含义为 t 的正部¹。图 5.1 的右下展示了函数 $h_3$ 的曲线。

我们通常想要更平滑的函数形式，这可以通过增加局部多项式的级数来实现。图 5.2 展示了在相同样本数据上拟合的四个分段三次多项式函数，它们在节点处的连续程度依次增加。

右下图的函数是连续的，而且在节点处有连续的一阶和二阶导数。这就是 三次样条（cubic spline） 函数。如果再多加一阶（三阶）的（导数）连续性约束，就会得到全局的三次多项式函数。可证明（练习 5.1）节点为 $\xi_1$ 和 $\xi_2$ 的三次样条可以用以下的基函数来表达：

这样就得到了 6 个基函数，对应着 6 个维度的（线性）函数空间。从分段三次函数和节点约束得到的参数个数可以佐证这个维度个数：3 个区域上各有 4 个参数，共 12 个参数；2 个节点处各有 3 个限制²，共 6 个限制；汇总后共有 $3\times4-2\times3=6$ 个自由参数。

更一般性地，包含了节点 $\xi_j,j=1,\dots,K$ 的 $M$ 阶（order）样条函数为一个分段 $M-1$ 次（degree）多项式函数³，其小于等于 $M-2$ 阶的导数都连续。三次样条的 $M$ 为 4。图 5.1 中的分段常数函数实际上为 1 阶样条，而连续的分段线性函数实际上为 2 阶样条。一般情况下的截断幂（truncated-power）基函数组为：

三次样条号称是肉眼无法看出节点不连续性的最低阶数样条。除非对连续的导数有特殊要求，否则基本没有选择高于三次的样条函数的必要。在实践中，最常用的阶数为 $M=1,2,4$。

固定节点的样条也被称为 回归样条（regression spline）。模型中要选择样条的阶数、节点的个数、以及节点的位置。一个简单的方法是将基函数的个数或自由度作为样条模型的参数，再从观测样本 $x_i$ 的范围决定节点的位置。例如，R 中的 bs(x, df=7) 语句⁴会生成一个三次样条函数的基函数在 $N$ 个点 $x$ 上的取值矩阵，其中包含了 $7-3=4$ 个内部节点⁵，分布在 $x$ 的（20%，40%，60%，80%）分位数上。或者可以更直接地定义样条函数：bs(x, degree=1, knots=c(0.2, 0.4, 0.6)) 生成含有 3 个内部节点的线性样条的基函数，返回值为一个 $N\times4$ 的矩阵。

指定了阶数和节点位置后，样条函数的空间为一个向量空间，因此存在多种等价的可张成这个空间的基（可类比普通的多项式空间）。截断幂（truncated-power）基函数虽然定义简洁，但大数的幂运算可能导致严重的取整问题。本章附录中介绍的 B-样条基函数，在节点数量 $K$ 很大时仍然可以快速地计算结果。

5.2.1 自然三次样条

多项式函数的拟合在边界附近的表现通常很不稳定，而且外推（extrapolation）可能会非常不可靠。这个问题在样条函数中更为严重。在边界节点外侧区域的多项式拟合要比在该区域中全局多项式拟合更加不可靠。最小二乘拟合样条函数的逐点方差可以清楚地说明这个问题（关于方差的计算参考下一小节中的例子）。图 5.3 比较了几个不同模型的逐点方差。在边界点附近，方差很明显地增加，而其中三次样条的幅度最大。

自然三次样条（natural cubic spline） 增加了额外的约束条件，即函数在边界节点外侧区域为线性。这减少了 4 个自由度（每个边界减少了 2 个），从而可在内部区域增加更多的节点。图 5.3 展示了方差的取舍。在边界附近的偏差会增大，然而函数在边界附近是线性的假设通常是合理的（因为这个区域本身提供的信息就不多）。

含有 $K$ 个节点的自然三次样条可被 $K$ 个基函数表述。可以从三次样条的基函数出发，再通过边界约束推导出新的基函数。例如，从上述的截幂基函数出发，可得（练习 5.4）：

其中

这些基函数在 $X\geq\xi_K$ 时的二阶和三阶导数均为 0。

5.5.2 实例：南非心脏病数据（续）

第 4.4.2 节用线性的对数几率模型拟合南非心脏病数据。以下用自然样条来拟合非线性的模型。模型的函数形式为：

其中每个 $\theta_j$ 为自然样条基函数 $h_j$ 向量对应的系数向量。⁶

模型中每一项中包括 4 个自然样条基函数。例如，若 $X_1$ 代表变量 $\text{sbp}$，则 $h_1(X_1)$ 为四个基函数组成的向量。这意味着有三个（而不是两个，因为排除了每个 $h$ 内部的常数项）内部节点和两个取值在数据极值的边界节点。内部节点位于变量取值范围的均匀分位点（uniform quantile）上。

由于 $\text{famhist}$ 为二分类因子，将其编码为 0/1 变量或哑变量，在模型拟合中只对应一个系数。

可以将所有 $p$ 个基函数向量（以及一个常数项）拼成一个大的向量 $h(X)$，模型的形式就变成了 $h(X)^T\theta$；将每项中的参数个数加起来，这个模型共有 $df=1+\sum_{j=1}^pdf_j$ 个参数。每个基函数在 $N$ 个样本点上计算取值，得到 $N\times df$ 维度的基函数矩阵 $\mathbf{H}$。这时的模型与之前的线性对数几率模型是一样的，第 4.4.1 节中的计算方法仍然适用。

模型的计算是一个后向逐步剔除的过程，以基函数组结构（$h_j(X)$）的方式保留或剔除变量，而不会只单独地排除某一个系数。筛选标准为 AIC 统计量（见第 7.5 节），将最终模型中的任意一组变量排除，均会导致 AIC 的上升（见表 5.1）。

图 5.4 展示了逐步回归选择的最终模型的曲线。其中的曲线为每个变量 $X_j$ 的 $\hat{f}_j(X_j)=h_j(X_j)^T\hat{\theta}_j$。协方差矩阵 $\operatorname{Cov}(\hat{\theta})=\mathbf{\Sigma}$ 的估计量为 $\hat{\mathbf{\Sigma}}=(\mathbf{H}^T\mathbf{W}\mathbf{H})^{-1}$，其中的 $\mathbf{W}$ 为对数几率回归中的对角权重矩阵。因此 $v_j(X_j)=\operatorname{Var}[\hat{f}_j(X_j)]=h_j(X_j)^T\hat{\mathbf{\Sigma}}_{jj}h_j(X_j)$ 为 $\hat{f}_j$ 的逐点方差函数，其中 $\operatorname{Cov}(\hat{\theta}_j)=\hat{\mathbf{\Sigma}}_{jj}$ 为 $\hat{\mathbf{\Sigma}}$ 中对应的子矩阵。图中每个曲线附近的阴影区域为 $\hat{f}_j(X_j)\pm 2 v_j(X_j)$。

在剔除变量方面，AIC 统计量要比似然差检验（偏差检验）统计量更宽容一些。变量 $\text{sbp}$ 和 $\text{obesity}$ 都被包含在了这个模型中，但没有包含在线性模型中。从图中可大致看出原因，这两个变量对模型的影响为非线性的。两者的影响看似出乎意料，但回溯数据的特性可能可以解释。这些测量值是在患者有过心脏病发作一段时间后获取的，而常常这些患者在患病后改善饮食和剩湖哦习惯，因此低水平的 obesity 和 sbp 反而体现出了患病风险的增加。表 5.1 给出了最终模型的概述。

表 5.1：经过分步剔除自然样条项后，最终的对数几率回归模型。其中的“LRT”列为该行的变量是否可排除模型的似然比检验统计量，实际上也是排除这个变量后与整体模型（第一行，none）的偏差（deviance）的变化。

5.2.3 实例：音素（phoneme）识别

在这个例子中，样条函数被用来降低而不是增加模型的灵活度；这在更广泛的函数建模领域经常使用。图 5.5 中的上图展示了两个音素“aa”和“ao”在 256 个频率上的各 15 个对数周期图（log-periodogram）。目标是利用这些数据来识别一个发音。之所以选择这两个音素是因为两者很难区分。

输入特征变量 $x$ 为一个长度为 256 的向量，可视为一个函数 $X(f)$ 在频率网格 $f$ 上的取值向量。在现实世界中，连续的模拟信号是频率的一个函数，这里是对这个函数一个取样。

图 5.5 中下图的灰色曲线线性对数几率回归模型系数的最大似然拟合值，训练样本为从 695 个“aa”和 1022 个 “ao”中抽取的 1000 个数据。图中将系数绘制成频率的函数曲线，而实际上可以从这个模型的连续形式来理解：

上式在样本中以下式来近似：

在两类音素的对数周期图中有差异的频率区间，会对应取值比较大的系数，这样的系数可以分辨出两个音素。

灰色曲线非常不平滑。由于输入的信号变量有比较强的正自相关，使得系数存在比较强的负自相关。另外，相较输入变量维度来说，样本数量不大，每个系数只对应大概 4 个观测值。

这个应用场景中，一个很自然的正则化约束是要求系数为频率的一个平滑的函数。图 5.5 的下图中红色曲线为用这些数据拟合的平滑系数曲线。平滑系数曲线不仅让两个音素的对比更容易理解，同时也使得分类结果更准确：

平滑的红色曲线是使用自然三次样条方法得出的。将系数函数表达为一个样条拓展的形式 $\beta(f)=\sum_{m=1}^Mh_m(f)\theta_m$。实际操作中，可写为 $\beta=\mathbf{H}\theta$，其中 $\mathbf{H}$ 为 $p\times M$ 的自然三次样条的基函数矩阵，基函数的定义域为频率的取值范围。此例中使用了 12 个基函数，节点等距地分布在代表频率的整数 $1,2,\dots,256$ 之间。由于 $x^T\beta=x^T\mathbf{H}\theta$，可以对输入特征变量进行转换 $x^*=\mathbf{H}^Tx$，然后用 $x^*$ 拟合线性对数几率回归系数 $\theta$。因而红色曲线即为 $\hat{\beta}(f)=h(f)^T\hat{\theta}$。

本节练习

练习 5.1

证明式 5.3 中的截断幂（truncated power）基函数是能够表达两个节点的三次样条函数的一组基函数。

练习 5.3

编写一个程序来复现第 145 页的图 5.3。

练习 5.4

考虑 $K$ 个内部节点的三次样条的截断幂表达式。定义：

证明自然三次样条中的边界条件（第 5.2.1 节）就相当于以下对系数的线性约束：

所以可推导出式 5.4 和 5.5 的形式。