Regression spline, penalized spline, and smoothing spline

  

1 全局多项式

1.1 一阶多项式与高阶多项式

  首先生成一个待拟合的数据，我们在(-1,1)区间上以三角函数为基础生成100个数据点，代码和图像如下：

# 安装与载入
# install.packages("ggplot2") 
library(ggplot2)
# install.packages("rms")
library(rms)

## Loading required package: Hmisc

## 
## Attaching package: 'Hmisc'

## The following objects are masked from 'package:base':
## 
##     format.pval, units

# install.packages("ctexart")

# 生成数据
set.seed(123)
n <- 100
x <- seq(-1, 1, length.out = n)
y <- sin(1.5 * pi * x) + rnorm(n, sd = 0.5)
ggplot(data = data.frame(x = x, y = y), aes(x, y)) + geom_point()

  先考虑线性回归：\(y_{i}=\beta_{0}+\beta_{1}x_{i}+\text{error}_{i}\)，等式的右边可以看作\(1\)和\(x\)以\(\beta_{0}\)和\(\beta_{1}\)作为权重线性相加，此时可以称这个模型的基（basis）是\(1\)和\(x\)。基函数图像如下图所示。

一阶线性基函数

  如果定义一个\(X\)矩阵并将这两个基函数写入，可以表示为： \[\boldsymbol{X}=\left[\begin{array}{cc} 1 & x_{1} \\ \vdots & \vdots \\ 1 & x_{n} \end{array}\right]\]

  同时， \[ \widehat{\beta}=(X^{T}X)^{-1}X^{T}y\;\;,\;\;\widehat{y}=X(X^{T}X)^{-1}X^{T}y=Hy \]

  我们将散点图使用线性回归进行拟合，代码和对应图像如下：

# 一阶多项式拟合

ggplot(data = data.frame(x = x, y = y), aes(x, y)) + geom_point() +
  geom_smooth(method="lm",formula='y ~ x')

  

  将模型推广到包括二次项，模型形式为：\(y_{i}=\beta_{0}+\beta_{1}x_{i}+\beta_{2}x_{i}^2+\text{error}_{i}\)，等式右边同样可以看作\(1\)、\(x\)、\(x^2\)以\(\beta_{0}\)、\(\beta_{1}\)、\(\beta_{2}\)为权重线性相加，这个模型的基即为\(1\)、\(x\)、\(x^2\)。基函数图像如下图所示：

二阶线性基函数

  如果同样将基函数写入\(X\)矩阵，可以表示为 \[ X={\left[\begin{array}{l l l}{1}&{x_{1}}&{x_{1}^{2}}\\ {\vdots}&{\vdots}&{\vdots}\\ {1}&{x_{n}}&{x_{n}^{2}}\end{array}\right]} \]

  系数和因变量拟合值的表达式同上，只是替换了\(X\)矩阵。

  如果我们将散点图使用二次项回归进行拟合，代码和对应图像如下:

# 二阶多项式拟合

ggplot(data = data.frame(x = x, y = y), aes(x, y)) + geom_point() +
  geom_smooth(method="lm",formula='y ~ poly(x,2)')

  

  将模型推广到包括k次项，模型形式为：\(y_{i}=\beta_{0}+\beta_{1}x_{i}+\beta_{2}x_{i}^2+...+\beta_{k}x_{i}^k+\text{error}_{i}\)，等式右边同样可以看作\(1\)、\(x\)、\(x^2\)…\(x^k\)以\(\beta_{0}\)、\(\beta_{1}\)、\(\beta_{2}\)…\(\beta_{0}\)为权重线性相加，这个模型的基即为\(1\)、\(x\)、\(x^2\)…\(x^k\)。

  如果同样将基函数写入\(X\)矩阵，可以表示为 \[ X=\left[\begin{array}{c c c c}{{1}}&{{x_{1}}}&{{x_{1}^{2}}}&{{\cdots}}&{{x_{1}^{k}}}\\ {{\vdots}}&{{\vdots}}&{{\vdots}}&{{\ddots}}&{{\vdots}}\\ {{\vdots}}&{{x_{n}}}&{{x_{n}^{2}}}&{{\cdots}}&{{x_{n}^{k}}}\end{array}\right] \]

  系数和因变量拟合值的表达式仍然同上，只是替换了\(X\)矩阵。

  如果我们将散点图使用更高阶的多项式进行拟合，代码和对应图像如下：

# 三阶多项式

ggplot(data = data.frame(x = x, y = y), aes(x, y)) + geom_point() +
  geom_smooth(method="lm",formula='y ~ poly(x,3)')

# 八阶多项式

ggplot(data = data.frame(x = x, y = y), aes(x, y)) + geom_point() +
  geom_smooth(method="lm",formula='y ~ poly(x,8)')

  从以上尝试可以看出，使用多项式回归模型时，低阶多项式不能很好地刻画数据的突然变动，阶数的逐渐提高可以更好地遍历数据，但其波动可能并不代表数据中的任何特征，更像一个”一般的高阶多项式”，有时一点小的扰动对拟合函数导数的影响远远大于对拟合函数本身的影响，导致在导数上无法反应真实数据之间应有的关系，经常出现与真实数据趋势相反且导数很大的情况。

  

1.2 正交多项式

  普通多项式具有远离正交性的缺点，非正交性会导致拟合曲线尤其在边界上数值不稳定，因此逐渐产生了多种正交多项式。正交性在向量中表示为点乘（inner product）为零，如果三个向量正交，那么他们之间是线性独立的，并且可以张成一个三维空间，这意味着空间内的每一点都可以用这三个向量的线性组合来表示。在多项式或者基函数中正交性即为任意两个基函数的乘积的积分为零，\(\textstyle{\int_{a}^{b}f_{i}(x)f_{j}(x)d x=0,i\neq j}\)。

  

1.2.1 勒让德多项式（Legendre polynomial）

  勒让德多项式是勒让德微分方程\((1-x^{2}){\frac{\partial^{2}P(x)}{\partial x^{2}}}-2x{\frac{\partial P(x)}{\partial x}}+\lambda P(x)=0,\lambda=n(x+1)\)的解，常用于物理学中球极坐标相关计算，在\([−1,1]\)区间上正交，其具体形式是： \[ \begin{array}{l c r}{P_{0}(x)=1}\\ {P_{1}(x)=x}\\ {P_{2}(x)={\frac{1}{2}}(3x^{2}-1)}\\ {P_{3}(x)={\frac{1}{2}}(5x^{3}-3x)}\\ {P_{4}(x)={\frac{1}{8}}(35x^{4}-30x^{2}+2)}\\ {P_{5}(x)={\frac{1}{8}}(63x^{5}-70x^{3}+15x)}\end{array} \]

  各阶Legendre多项式图像如下：

legendre多项式基函数

  我们分别使用Legendre三阶、八阶多项式拟合数据，代码和对应图像如下：

# legendre多项式

#设定阶数
degree <- 3

# 定义Legendre多项式矩阵
x_norm <- (x - mean(x)) / sd(x) #标准化x
legendre_poly <- matrix(0, nrow = length(x), ncol = degree + 1) #构造矩阵
legendre_poly[, 1] <- rep(1, length(x))
legendre_poly[, 2] <- x_norm
legendre_poly[, 3] <- (3 * x_norm^2 - 1) / 2
legendre_poly[, 4] <- (5 * x_norm^3 - 3 * x_norm) / 2

# 拟合
fit <- lm(y ~ legendre_poly - 1)

# 画出散点图和拟合图
ggplot(data = data.frame(x = x, y = y),aes(x,y)) +
  geom_point() +
  geom_line(mapping = aes(x, predict(fit)))

#设定阶数
degree <- 8

# 定义Legendre多项式矩阵
x_norm <- (x - mean(x)) / sd(x) #标准化x
legendre_poly <- matrix(0, nrow = length(x), ncol = degree + 1) #构造矩阵
legendre_poly[, 1] <- rep(1, length(x))
legendre_poly[, 2] <- x_norm
legendre_poly[, 3] <- (3 * x_norm^2 - 1) / 2
legendre_poly[, 4] <- (5 * x_norm^3 - 3 * x_norm) / 2
legendre_poly[, 5] <- (35 * x_norm^4 - 30 * x_norm^2 + 3) / 8
legendre_poly[, 6] <- (63 * x_norm^5 - 70 * x_norm^3 + 15 * x_norm) / 8
legendre_poly[, 7] <- (231 * x_norm^6 - 315 * x_norm^4 + 105 * x_norm^2 - 5) / 16
legendre_poly[, 8] <- (429 * x_norm^7 - 693 * x_norm^5 + 315 * x_norm^3 - 35 * x_norm) / 16
legendre_poly[, 9] <- (6435 * x_norm^8 - 12012 * x_norm^6 + 6930 * x_norm^4 - 1260 * x_norm^2 + 35) / 128

# 拟合
fit <- lm(y ~ legendre_poly - 1)

# 画出散点图和拟合图
ggplot(data = data.frame(x = x, y = y),aes(x,y)) +
  geom_point() +
  geom_line(mapping = aes(x, predict(fit)))

  

1.2.2 埃米尔特多项式（Hermite polynomial）

  埃米尔特多项式是埃米尔特方程\(\frac{\partial H_{n}^{2}}{\partial^{2}H_{n}}\!-\!2x H_{n}\!+\!n H_{n}\!=\!0\)的解，常作为量子谐振子哈密尔顿解的一部分出现，在\((-\infty,\infty)\)区间上正交，其具体形式是： \[ \begin{array}{l}{{H_{0}\bigl(x\bigr)=1}}\\ {{H_{1}\bigl(x\bigr)=2x}}\\{{H_{2}\bigl(x\bigr)=4x^{2}-2x+2}}\\ {{H_{3}\bigl(x\bigr)=8x^{3}-12x}}\\ {{H_{4}\bigl(x\bigr)=16x^{4}-48x^{2}+12}}\\ {{H_{5}\bigl(x\bigr)=32x^{5}-160x^{3}+120x}}\end{array} \]

Hermite多项式基函数

  使用埃米尔特三阶、八阶多项式拟合数据的代码和对应图像为：

# Hermite多项式

# 定义Hermite多项式
hermite <- function(n, x) {
  if (n == 0) {
    return(rep(1, length(x)))
  } else if (n == 1) {
    return(2 * x)
  } else {
    h0 <- rep(1, length(x))
    h1 <- 2 * x
    for (i in 2:(n - 1)) {
      hi <- 2 * x * h1 - 2 * i * h0
      h0 <- h1
      h1 <- hi
    }
    return(h1)
  }
}

# 定义阶数
degree <- 3

# 计算Hermite多项式矩阵
hermite_matrix <- matrix(nrow = length(x), ncol = degree + 1)
for (i in 0:degree) {
  hermite_matrix[, i + 1] <- hermite(i, x)
}

# 拟合
fit <- lm(y ~ hermite_matrix - 1)

# 画出散点图和拟合图
ggplot(data = data.frame(x = x, y = y), aes(x, y)) + 
  geom_point() + 
  geom_line(mapping = aes(x, predict(fit)))

# 定义阶数
degree <- 8

# 计算Hermite多项式矩阵
hermite_matrix <- matrix(nrow = length(x), ncol = degree + 1)
for (i in 0:degree) {
  hermite_matrix[, i + 1] <- hermite(i, x)
}

# 拟合
fit <- lm(y ~ hermite_matrix - 1)

# 画出散点图和拟合图
ggplot(data = data.frame(x = x, y = y), aes(x, y)) + 
  geom_point() + 
  geom_line(mapping = aes(x, predict(fit)))

  

1.2.3 三角多项式（Trigonometric polynomials）

  三角多项式是三角函数\(sin(x)\)和\(cos(x)\)的有限线性组合，多被用于周期性数据的拟合，表达式为\(f(t)=\sum_{n=0}^{N}\Bigl(a_{n}\mathrm{cos}{(}n t)+b_{n}\mathrm{sin}(n t){\Bigr)}\)，另一种表达形式是\(f(t)=\sum_{n=-N}^{N}(c_{n}e^{i n t})\)。前一种形式通过\({{cos\left(n t\right)=(e^{i n t}+e^{-i n t})/2}}\)和\({{sin\left(n t\right)=(e^{i n t}-e^{-i n t})/2}}\)可以转换为后一种形式，后一种表达形式更常使用，因为\(\Bigl\{e^{ i n t}\Bigl\}_{-\infty,\infty}=\left\{\ldots,e^{-2i t},e^{-i t},1,e^{i t},e^{2i t}\ldots\right\}\)，此时如果要构造正交的基底，对系数\(c_{n}\)的计算会更方便。

  三角多项式拟合数据代码和对应图像为：

# 三角多项式

# 设定频数：选用几个谐波来构建多项式
n_harm <- 4

# 定义三角多项式矩阵
X <- matrix(0, nrow = length(x), ncol = 2 * n_harm)
for (j in 1:n_harm) {
  X[, 2*j - 1] <- sin(j * x)
  X[, 2*j] <- cos(j * x)
}

# 拟合
fit <- lm(y ~ X - 1)

# 画出散点图和拟合图
ggplot(data = data.frame(x = x, y = y), aes(x, y)) + 
  geom_point() + 
  geom_line(mapping = aes(x, predict(fit)))

# 三角多项式

# 设定频数：选用几个谐波来构建多项式
n_harm <- 8

# 定义三角多项式矩阵
X <- matrix(0, nrow = length(x), ncol = 2 * n_harm)
for (j in 1:n_harm) {
  X[, 2*j - 1] <- sin(j * x)
  X[, 2*j] <- cos(j * x)
}

# 拟合
fit <- lm(y ~ X - 1)

# 画出散点图和拟合图
ggplot(data = data.frame(x = x, y = y), aes(x, y)) + 
  geom_point() + 
  geom_line(mapping = aes(x, predict(fit)))

  

1.2.4 伯恩施坦多项式 (Bernstein polynomials)

  伯恩施坦多项式首次被Bernstein用于Weierstrass近似定理的构造性证明，在\([0,1]\)区间上正交，与贝塞尔曲线紧密相关。

  n阶伯恩施坦多项式的基函数定义为：

\[ b_{\nu,n}\Big(x\Big)=C_{n}^{v}x^{\nu}\bigl(1-x\bigr)^{n-\nu},\nu=0,...,n,C_{n}^{v}=\frac{n!}{\nu!\Big(n-\nu\Big)!} \]

  具体写出前几项有：

\[ \begin{array}{l l}{{b_{0,0}(x)=1}}&{{}}&{{}}\\ {{b_{0,1}(x)=1-x}}&{{b_{1,1}(x)=x}}\\ {{b_{0,2}(x)=(1-x)^{2}}}&{{b_{1,2}(x)=2x(1-x)}}&{{b_{2,2}(x)=x^{2}}}\\ {{b_{0,3}(x)=(1-x)^{3}}}&{{b_{1,3}(x)=3x(1-x)^{2}}}&{{b_{2,3}(x)=3x^{2}(1-x)}}&{{b_{3,3}(x)=x^{3}}}\end{array} \]

  因此拟合曲线就写为\(B_{n}(x)=\sum_{\nu\equiv0}^{n}\beta_{\nu}b_{\nu,n}(x)\)。

  伯恩施坦多项式的基函数有以下重要性质：

  （1）对称性，\(b_{\nu,n}{(}1-x{)}=\,b_{n-\nu,n}{(}x{)}\)；

  （2）导数可以写成两个低阶基的加总，\(b_{\nu,n}^{\prime}\!\left(x\right)=n(b_{\nu-1,n-1}(x)\!-\!b_{\nu,n-1}\!\left(x\right)\!)\)；

  （3）单位分解（partition of unity），\(\sum_{\nu=0}^{n}b_{\nu,n}(x)=1\)。

  三阶、四阶Bernstein多项式基函数图像如下：

三阶Bernstein多项式基函数

四阶Bernstein多项式基函数

  我们使用伯恩施坦三阶、八阶多项式拟合数据，相关代码和拟合图像为

# Bernstein多项式

# 定义Bernstein多项式
bernsteinPoly <- function(k, n, x) {
  choose(n, k) * x^k * (1 - x)^(n - k)
}

# 定义阶数
degree <- 3

# 生成Bernstein多项式矩阵
X <- matrix(0, nrow = length(x), ncol = degree + 1)
for (i in 0:degree) {
  X[, i + 1] <- bernsteinPoly(i, degree, x)
}

# 拟合
fit <- lm(y ~ X)

# 画出散点图和拟合图
ggplot(data = data.frame(x = x, y = y), aes(x, y)) + 
  geom_point() + 
  geom_line(mapping = aes(x, predict(fit)))

# 定义阶数
degree <- 8

# 生成Bernstein多项式矩阵
X <- matrix(0, nrow = length(x), ncol = degree + 1)
for (i in 0:degree) {
  X[, i + 1] <- bernsteinPoly(i, degree, x)
}

# 拟合
fit <- lm(y ~ X)

# 画出散点图和拟合图
ggplot(data = data.frame(x = x, y = y), aes(x, y)) + 
  geom_point() + 
  geom_line(mapping = aes(x, predict(fit)))

  贝塞尔曲线就是以伯恩施坦多项式作为基底、控制点作为系数画出的曲线，而B样条可以理解为贝塞尔曲线的拼接，从而克服贝塞尔曲线改变控制点会影响整条曲线的不足。有关贝塞尔曲线和B样条的内容会在第2节和第3节详细说明。

  

2 局部多项式

2.1 低阶截断多项式

  为适应在不同区间的潜在趋势变化很大（甚至反转）的数据结构，同时缓解使用多项式拟合会出现的边界数值不稳定问题，接下来引入局部多项式。

  首先最自然的是”折线”类型的基函数，称为截断基函数（Truncated basis function）。一阶的截断基函数形如\(\left(x-0.6\right)_{+}\),具体含义是当\(x<0.6\)时取\(0\)，当\(x>0.6\)时取\(x-0.6\)。从进行模型改良的比较基准的角度说，一般将截断基函数与多项式基函数保持同阶，即在一阶多项式模型（线性模型）中加入截断项时一般用一阶截断项，以此类推。以一阶截断基为例，加入线性模型后表示为\(y_{i}=\beta_{0}+\beta_{1}x_{i}+\beta_{11}\left(x-0.6\right)_{+}+e r r o r_{i}\)。等式右边同样可以看作\(1\)、\(x\)、\(\left(x-0.6\right)_{+}\)以\(\beta_{0}\)、\(\beta_{1}\)、\(\beta_{11}\)为权重线性相加，这个模型的基即为\(1\)、\(𝑥\)、\(\left(x-0.6\right)_{+}\)，基函数图像如下图所示：

一阶截断多项式基函数

  如果同样将基函数写入\(X\)矩阵，可以表示为 \[X={\left[\begin{array}{l l l}{1}&{x_{1}}&{\left(x_{1}-0.6\right)_{+}}\\ {\vdots}&{\vdots}&{\vdots}\\ {1}&{x_{n}}&{\left(x_{n}-0.6\right)_{+}}\end{array}\right]}\]

  系数和因变量拟合值的表达式仍然相同，只是替换了\(X\)矩阵。

  

  带截断基的模型有两个简单的推广方向，其一是增加结点，主要针对数据趋势变动较多的区间；其二是增加阶数，主要为了达到更好的拟合效果。

  首先在增加结点方面，同样以一阶为例，假如在\((0.5,1)\)区间内以\(0.05\)为间隔均匀铺入结点可得基函数\(\left(x-0.5\right)_{+}\)、\(\left(x-0.55\right)_{+}\)…\(\left(x-0.95\right)_{+}\)，则基函数图像如下图所示

一阶截断多项式基函数增加结点

  方法之二是增加阶数。线性样条的拟合结果在结点较密时看起来像一条曲线，但实际应是一整条连续的分段折线。要想使拟合函数避免分段线性，或者想使用增加结点以外的方式使曲线变得平滑，更好地拟合图像中的波峰和波谷，可以引入二次项构造二阶样条。

  将线性样条中每一个截断项换为二次项并加入一般二次项，可得模型\(y_{i}\!=\!\beta_{0}\!+\!\beta_{1}x_{i}\!+\!\beta_{2}x_{i}^{2}\!+\!\sum_{k=1}^{K}\!\beta_{2k}\Big(x\!-\!\kappa_{k}\Big)_{+}^{2}+e r r o r_{i}\)，其中基函数为\(1\)、\(x\)、\(x^2\)、\(\Big(x\!-\!\kappa_{1}\Big)_{+}^{2}\)…\(\Big(x\!-\!\kappa_{k}\Big)_{+}^{2}\)。换为二次项后，基函数在结点处有连续一阶导数，直观表示为基函数没有尖锐转角。此时系数和因变量拟合值的表达式仍然同上，只是替换了\(X\)矩阵。

  在实际中，结点较少时二次样条拟合能比线性样条更有效地平滑峰值和谷值，不再存在折角，但是在使用线性样条时如果选取足够多的结点并设置惩罚项，两种拟合效果的差异就可以忽略不计，惩罚项将在第三节详细介绍。

  如果要得到更平滑的拟合，或者对拟合函数的导数有更高的要求，那可以继续提高样条的阶数，但现实中一般不使用超过三阶的样条。\(p\)阶截断幂模型表示为\(y_{i}=\beta_{0}+\beta_{1}x_{i}+\beta_{2}x_{i}^{2}+\dots+\beta_{p}x^{p}+\sum_{k=1}^{p}\beta_{p k}{\Big(}x-\kappa_{k}{\Big)}_{+}^{p}+e r r o r_{i}\)，如果同样将基函数写入\(𝑿\)矩阵，可以表示为

\[ X={\left[\begin{array}{l l l l l l}{1}&{x_{1}}&{\cdots}&{\left(x_{1}-\kappa_{1}\right)_{+}^{p}}&{\cdots}&{\left(x_{1}-\kappa_{k}\right)_{+}^{p}}\\ {\vdots}&{\vdots}&{\vdots}&{\vdots}&{\vdots}&{\vdots}\\ {1}&{x_{n}}&{\cdots}&{\left(x_{n}-\kappa_{1}\right)_{+}^{p}}&{\cdots}&{\left(x_{n}-\kappa_{k}\right)_{+}^{p}}\end{array}\right]} \]   共有\(1+𝑘+𝑝\)列。系数和因变量拟合值的表达式仍然同上，只是替换了\(𝑿\)矩阵。

  

  我们使用三阶截断幂基拟合数据，选取均匀分布的6个结点，相关的代码和拟合结果为：

# 截断幂基

knots <- quantile(x, probs = seq(0, 1, 0.2))

# 三阶截断幂
b=cbind(1,x,x^2,x^3,outer(x,knots,function(p,q) (p-q)^3*(p>q)))
fit5<-lm(y ~ b, data = data.frame(x = x, y = y))
ggplot(data = data.frame(x = x, y = y),aes(x = x, y = y)) + geom_point() +
  geom_line(mapping = aes(x, predict(fit5)))

# 截断幂基

knots <- quantile(x, probs = seq(0, 1, 0.2))

# 八阶截断幂
b=cbind(1,x,x^2,x^3,x^4,x^5,x^6,x^7,x^8,outer(x,knots,function(p,q) (p-q)^8*(p>q)))
fit3<-lm(y ~ b, data = data.frame(x = x, y = y))
ggplot(data = data.frame(x = x, y = y),aes(x = x, y = y)) + geom_point() +
  geom_line(mapping = aes(x, predict(fit3)))

  如上所述的分段线性函数就可以称为结点在\(\kappa_{1}\)…\(\kappa_{k}\)处的样条函数，“样条”得名于分段函数对数据的灵活适应。值得重视的是，截断基的使用极大依赖结点的选择。从截断基的图像直观理解，结点将截断基分段的两部分通过这个点连接起来；从拟合数据的角度理解，结点应该更多地出现在数据变化较剧烈或趋势有明显改变的区间。如果在整个区间内选择较密的一组结点将提升拟合效果，但容易过拟合，意味着会将较小、较不重要的信息也纳入模型进行处理。

  下图分别为均匀选取11个结点、均匀选取21个结点的三阶截断幂样条对数据的拟合结果。

# 增加结点
knots <- quantile(x, probs = seq(0, 1, 0.1))

b=cbind(1,x,x^2,x^3,outer(x,knots,function(p,q) (p-q)^3*(p>q)))
fit5<-lm(y ~ b, data = data.frame(x = x, y = y))
ggplot(data = data.frame(x = x, y = y),aes(x = x, y = y)) + geom_point() +
  geom_line(mapping = aes(x, predict(fit5)))

# 增加结点
knots <- quantile(x, probs = seq(0, 1, 0.05))

b=cbind(1,x,x^2,x^3,outer(x,knots,function(p,q) (p-q)^3*(p>q)))
fit5<-lm(y ~ b, data = data.frame(x = x, y = y))
ggplot(data = data.frame(x = x, y = y),aes(x = x, y = y)) + geom_point() +
  geom_line(mapping = aes(x, predict(fit5)))

  

2.2 B样条

  使用截断幂基进行样条的说明比较易于理解，且如果仔细选择结点、使用惩罚系数压缩使用，可以运用到实践中，但由于截断幂基仍然存在着远离正交性的缺点，表现为在结点数量过多、惩罚参数很小时会导致数值不稳定，因此在实践中建议使用具有更稳定数值性质的等效基。

  其中最常见的选择是来源于贝塞尔曲线（Bezier Curve）的B样条。

  

  首先简要介绍一下贝塞尔曲线。贝塞尔曲线是伯恩施坦多项式的延伸，最开始由法国工程师皮埃尔·贝塞尔用于汽车设计，后主要用于二维图形应用程序中的数学曲线绘制，PS中的钢笔工具就运用了这一方法。曲线由起始点、终止点和控制点组成，通过调整控制点来调整贝塞尔曲线的形状，这些点称为de Boor点。直观来说贝塞尔曲线就是一个点从起始点到终止点、以控制点为条件划过的轨迹。

  一阶贝塞尔曲线没有控制点，若起始点与终止点分别为\(P_0\)和\(P_1\)，方程为\(B_{1}(t)=P_{0}+(P_{1}-P_{0})t=(1-t)P_{0}+tP_{1}\)，\(t\in[0,1]\)。

  二阶贝塞尔曲线有一个控制点，若起始点、控制点、终止点分别为\(P_0\)、\(P_1\)、\(P_2\)，则分别在\(P_0P_1\)、\(P_1P_2\)两条连线上再设定两个控制点， \[\begin{array}{l c r}{{P_{0}^{\prime}=\left(1-t\right)P_{0}+t P_{1}}}\\ {{P_{1}^{\prime}=\left(1-t\right)P_{1}+t P_{2}}}\end{array}\]

  进而得到二阶贝塞尔曲线 \[\begin{array}{c}{{B_{2}\Big(t\Big)=\left(1-t\right)P_{0}^{\prime}+tP_{1}^{\prime}}}\\ {{=\left(1-t\right)\left(\left(1-t\right)P_{0}+tP_{1}\right)+t\left(\left(1-t\right)P_{1}+tP_{2}\right)}}\\{{=\left(1-t\right)^2P_{0}+2t\left(1-t\right)P_{1}+t^2P_{2},t\in[0,1]}}\end{array}\]

  继续推导可得n阶贝塞尔曲线\[B\left(t\right)=\sum_{i=0}^{n}C_{h}^{i}P_{i}t^{i}\left(1-t\right)^{n-i},\ \,t\in\left[0,1\right]\]

  从推导过程可以看出\(n\)阶贝塞尔曲线会经过\(n-1\)次再选点完成降阶，从函数形式可以看出\(n\)阶贝塞尔曲线每一项的系数是二项式\((t+(1-t))^n\)的展开公式。

  贝塞尔曲线图像如下所示：

一阶贝塞尔曲线

二阶贝塞尔曲线

三阶贝塞尔曲线

  贝塞尔曲线有几个重要的特点：

  （1）贝塞尔曲线的阶数是总控制点数减1，且各项系数之和为0；

  （2）移动一个控制点会影响整条曲线的形状，但画出的图像始终在包含了所有控制点的最小凸多边形中，这个凸多边形称为de Boor多边形，这一性质对路径规划来说很重要；

  （3）起始点处的一阶导数等于起始点与第一个控制点连线的斜率，终止点处的一阶导数等于最后一个控制点与终止点连线的斜率，因此两条贝塞尔曲线进行拼接时，只要离拼接点最近的两个控制点与拼接点三点共线即可满足拼接后是导数连续的。

  

  为克服贝塞尔曲线”牵一发而动全身”并且控制点和阶数绑定的缺点，B样条采用贝塞尔曲线的拼接，从而可以指定阶次，并且做到移动控制点只改变曲线的部分形状。B样条基函数代替了贝塞尔曲线里的二项式展开项，基函数表达式为 \[ N_{i,0}(u)=\left\{ \begin{aligned} 1&,if u_i\leq u\leq u_{i+1} \\ 0&,otherwise \end{aligned} \right.\\ N_{i,p}(u)=\frac{u-u_{i}}{u_{i+p}-u_{i}}*N_{i,p-1}(u)+\frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}*N_{i+1,p-1}(u) \]   上述公式称为Cox-de Boor递归公式，其中\(i\)表示基函数所在的结点区间，\(p\)表示基函数的阶数。在\([0,1]\)上具有七个不均匀结点的一阶、二阶、三阶B样条基的图像展示如下。

B样条基函数

  具体来看B样条基函数的计算方式。零阶B样条基函数是阶梯函数，在各自的区间上取1，其他区间取0，如果用零阶B样条基函数直接拟合数据，最终得到的拟合函数是一个个离散控制点。一阶B样条曲线就是把这些控制点连起来的曲线。更高阶的样条曲线我们使用三角计算格式进行\(N_{(i,p)}(u)\)的计算说明，第一列是结点分割的区间，第二列是零阶B样条基函数，即阶梯函数，两两线性插值得到一阶B样条基函数，继续插值得到二阶基函数，以此类推。B样条基的三角计算格式如下图：

B样条基函数的三角计算式1

  从B样条基函数图上可以看出，阶数决定了每一个基函数在\(x\)轴横跨的范围，例如一阶B样条基函数横跨一个结点（边缘结点除外），且导数在这个结点处转向，二阶B样条基函数横跨两个结点（区间边缘除外），二阶导数在这两个结点上转向。同时，根据三角计算式我们也可以发现这样的规律，某个基函数的非零定义域是它向左下和左上画一个三角形后竖直边上对应的区间。因此，基函数\(N_{i,p}(u)\)在\(\left[u_{i},u_{i+p+1}\right)\)这\(p+1\)个小区间上非零，即第i个结点仅在\(\left[u_{i},u_{i+p+1}\right)\)区间内被使用，这个性质被称为局部支撑（local support）。

  在相反的方向上，给定一个结点区间，要知道这个区间最终的拟合函数会由哪些基函数构成，需要从这个区间出发向右上和右下画边长为\(p+1\)的等边三角形后竖直边对应的基函数，即从\(N_{i-p,p}(u)\)到\(N_{i,p}(u)\)。\(p\)阶样条在一个区间上至多有\(p+1\)个基函数，边缘区间上的基函数数量会减少。

  以上两个问题的图示说明如下图：

B样条基函数的三角计算式2

B样条基函数的三角计算式3

  在导数上，首先因为继承了贝塞尔曲线的性质，B样条的导数还是一个B样条。B样条公式为：\[B(u)\!=\!\sum_{i=0}^{n}P_{i}N_{i,p}(u)\]

  基函数求导等于低阶基函数求导的线性相加，因此\[\frac{d}{d u}C(u)\!=\!C^{\prime}(u)\!=\!\sum_{i=0}^{n-1}N_{i+1,p-1}(u)Q_{i}\]

  总结一下\(p\)阶B样条的性质：

  （1）\(p\)阶B样条由\(p+1\)个区间多项式块组成，每个区间多项式块的阶数是\(p\)，多项式块在\(p\)个内结点处连接，在连接点处高达\(p-1\)阶导数是连续的；

  （2）B样条在\(p+2\)个区间上是正的，其他区间都为零；

  （3）除了在边界上，它与边上的\(2p\)个多项式块重叠；

  （4）在给定\(x\)处，有\(p+1\)个B样条是非零的。

  

  B样条与截断幂基对比：

  首先，B样条因为相比截断幂基更具正交性而在数值上更稳定；其次，截断幂基构造的拟合曲线在每个结点处改变函数式（加上新的一项），但B样条基通常在结点处改变导数性质。

  如果两组基能够张成（span）同样的函数集，即它们的所有线性组合是相同的函数集，就称这两个基函数等价。同阶的截断幂基和B样条是等价的，如果将B样条基写入\(X\)矩阵的列，那么拟合效果将与结点相同、阶数相同的截断幂基函数相同。在数学上，两组基等价意味着一组基将列写入得到的\(X\)矩阵乘以一个可逆方阵，可以得到另一组基将列写入得到的\(X\)矩阵。假设\(X_T\)是截断幂基对应的\(X\)矩阵，\(X_B\)是相同阶数、相同结点位置的B样条基对应的\(X\)矩阵，则可以表示为\(X_B=X_T L_p\)。

  

  自然样条：

  对三次样条的一个常见修正方法是使用自然三次样条基，即在三次样条的基础上对边界结点施加线性限制，这一限制是通过要求样条在边界结点满足二阶导和三阶导都为零实现的。对于\(k\)阶的样条，在满足结点间有连续\(k-1\)阶导的前提下，要求在最边缘的区间将阶数降为\((k-1)/2\)。对于截断幂基可以构造自然样条变体，此时只需要与结点数相同个数的基函数即可。

  三阶自然样条基函数图像如下图所示：

自然样条基函数

3 B样条及其变体

3.1 回归样条

  前文所述的分段多项式回归即为样条回归。通过将整体划分成多个区间，并使用分段多项式来拟合，避免了使用全局多项式时模型容易变得复杂的不足，也克服了全局多项式拟合情况易受单个样本点变化影响的缺陷。使用分段多项式时，我们希望不同区间之间的拟合曲线连接较为平滑，因此需要要求区间两侧多项式有阶数-1阶连续导数，同时我们希望分段多项式具有较好的正交性，这能带来更好的数值稳定性质。因此，最常用的基底为B样条基。

  使用来自\(splines\)包的\(bs()\)函数，\(bs()\)函数即表示使用B样条拟合，完整代码为：\[bs(x, df = NULL, knots = NULL, degree = 3, intercept = FALSE, Boundary.knots = range(x))\]

  \(df\)用于指定自由度，指定后函数会等分位距选择\(df-degree(-1)\)个结点，\(knots\)可以根据数据散点图观察结果手动指定，也可以根据数据范围均匀分布，还可以使用\(quantile()\)函数进行分位数指定；\(degree\)用于指定阶数，没有指定时默认值为3；\(intercept\)用于指定是否包含截距，默认不包含；\(Boundary.knots\)用于设定边界结点。

  使用B样条和6个均匀结点对数据进行拟合的代码和图像为：

library(ggplot2)
library(splines)

# 设置 B 样条的结点
knots <- quantile(x, probs = seq(0, 1, 0.2))

# 计算 B 样条拟合曲线
fit <- lm(y ~ bs(x, degree=3, knots = knots))

# 绘制数据散点图和拟合曲线
ggplot(data = data.frame(x = x, y = y),aes(x = x, y = y)) + geom_point() +
  geom_line(mapping = aes(x, predict(fit))) 

  

3.2 惩罚B样条

  B样条的结点数量和位置的最优比较复杂，在平滑样条（smoothing spline）中，Reinsch (1967)提出了基于二阶导平方的积分的惩罚项，O’Sullivan(1986,1988)将其推广到B样条中，提倡使用较多的结点，并为了防止过拟合加入设置在二阶导数上的惩罚。

  如果选择\(K\)个内部结点并使用三阶B样条，则共有\(K+8\)个结点（其中前四个重合在区间左端，最后四个重合在区间右端），这些结点定义了\(K+4\)个B样条基函数，\(B\)矩阵大小即为\(n×(K+4)\)，惩罚矩阵\(D\)的大小即为\((K+4)×(K+4)\)， \[D_{i j}=\int_{a}^{b}B_{i}^{\prime\prime}(x)\,B_{j}^{\prime\prime}(x)d x\]

  于是拟合函数为\[\hat{f}(x;\lambda)=B(B^{T}B+\lambda D)^{-1}B^{T}y\]

  D矩阵可简化为\[D=\left(\tilde{B}^{\prime\prime}\right)^{T}d i a g(w)\tilde{B}^{\prime\prime}\]

  其中\(\tilde{B}^{\prime\prime}\)为\(3(K+7)×(K+4)\)的矩阵，第\((i,j)\)个元素\(B_{j}^{\prime\prime}({\tilde{x}}_{i})\)，\({\tilde{x}}_{i}\)是第\({\tilde{x}}\)的第i个元素，\(w\)是\(3(K+7)×1\)的向量。 \[{\overline{{x}}}=(\kappa_{1},\frac{\kappa_{1}+\kappa_{2}}{2},\kappa_{2},\kappa_{2},\kappa_{3},\kappa_{3},\ldots,\kappa_{K+7},{\frac{\kappa_{K+7}+\kappa_{K+8}}{2}},\kappa_{K+8})\] \[w=\left(\frac{1}{6}\left(\Delta\kappa\right)_{1},\frac{4}{6}\left(\Delta\kappa\right)_{1},\frac{1}{6}\left(\Delta\kappa\right)_{1}...\frac{1}{6}\left(\Delta\kappa\right)_{K+7},\frac{4}{6}\left(\Delta\kappa\right)_{K+7},\frac{1}{6}\left(\Delta\kappa\right)_{K+7}\right)\]

  

3.3 P样条：B样条的另一种惩罚项

  Eilers and Marx(1996)提出使用相对较多的等距结点和相邻B样条系数的差值惩罚，并称之为P样条。P样条没有边界效应（拟合曲线在边界的拓展通常伴随着向零的弯曲），是广义线性回归模型的直接拓展。这一方法实际上简化并推广了O’Sullivan的思想，相较于使用拟合曲线的高阶导数的平方的积分构造惩罚项，选择了简单地对相邻B样条的系数本身进行差分惩罚，可以很容易在广义回归模型中加入任意阶的惩罚。

  首先，设\(B_j(x;p)\)是给定等距结点下\(x\)处的第\(j\)个\(p\)阶\(B\)样条，拟合曲线是每个数据点的线性相加，\({\widehat{y}}=\sum_{j=1}^{n}{\widehat{a}}_{j}B_{j}(x;p)\)，则B样条的一阶导数是\[h\sum_{j}a_{j}B_{j}^{\prime}(x;p)\!=\!\sum_{j}a_{j}B_{j}\Big(x;p\!-\!1\Big)\!-\!\sum_{j}a_{j+1}B_{j+1}\Big(x;p\!-\!1\Big)\!-\!\sum_{j}\Delta a_{j+1}B_{j}\Big(x;p\!-\!1\Big)\]

  其中h是结点之间的距离，\(\Delta a_{j+1}=a_{j+1}-a_{j}\)。通过归纳法可以得到二阶导数为\[h^{2}\sum_{j}a_{j}B_{j}^{\prime\prime}(x;p)=\sum_{j}\Delta^{2}a_{j}B_{j}(x;p-2)\]

  其中，\(\Delta^{2}a_{j}=\Delta\Delta a_{j}=a_{j}-2a_{j-1}+a_{j-2}\)。

  

  接下来介绍惩罚项的具体形式和最小化目标。如果有m个数据点和n个B样条，则最小化目标为\[S=\sum_{i=1}^{m}\{y_{i}-\sum_{j=1}^{n}a_{j}B_{j}(x_{i})\}^{2}\]

  如果进行直观对比，O’Sullivan建议使用较多结点但是限制曲线的变化程度，提出了基于二阶导数的惩罚项，最小化目标为 \[S=\sum_{i=1}^{m}\{y_{i}-\sum_{j=1}^{n}a_{j}B_{j}(x_{i})\}^{2}+\lambda\int_{x m i n}^{x m a x}\left\{\sum_{j=1}^{n}a_{j}B_{j}^{\prime\prime}(x)\right\}^{2}d x\]

  Eilers and Marx提出使用相邻B样条的（高阶）有限差分来进行惩罚，最小化目标为\[S=\sum_{i=1}^{m}\{y_{i}-\sum_{j=1}^{n}a_{j}B_{j}(x_{i})\}^{2}+\lambda\sum_{j=k+1}^{n}(\Delta^{k}a_{j})^{2}\]

  这种处理可以将问题的维数降到B样条的个数\(n\)，而不是样本个数\(m\)。差分惩罚是对导数的平方积分的一个很好的离散逼近，有了这个惩罚项之后数据的矩是守恒的。

  

  接下来将说明Eilers and Marx使用二阶差分的惩罚与O’Sullivan使用二阶导的惩罚之间的关联性。O’Sullivan在三阶B样条中使用\[h^{2}P=\lambda\int_{x m i n}^{x m a x}\{\sum_{j}a_{j}B_{j}^{\prime\prime}(x;3)\}^{2}d x\]

  B样条的导数性质可将导数算出来\[h^{2}P=\lambda\int_{x m i n}^{x m a x}\left(\sum_{j}\Delta^{2}a_{j}B_{j}\Big(x;1\Big)\right)^{2}d x=\lambda\int_{x m i n}^{x m a x}\sum_{j}\sum_{k}\Delta^{2}a_{j}\Delta^{2}a_{k}B_{j}\Big(x;1\Big)B_{k}\Big(x;1\Big)d x\]

  大部分一次样条的乘积都会消失，因为当\(𝑗\)等于\(𝑘−1\)、\(𝑘\)、\(𝑘+1\)时乘积才不为零，因此 \[ h^{2}P=\lambda\int_{x m i n}^{x m a x}\left\{\Delta^{2}a_{j}B_{j}\Big(x;1\Big)\right\}^{2}+2\sum_{j}\Delta^{2}a_{j}\Delta^{2}a_{j-1}B_{j}\Big(x;1\Big)B_{j-1}\Big(x;1\Big)\Big\}d x\\ =\lambda\sum_{j}\Bigl(\Delta^{2}a_{j}\Bigr)^{2}\int_{x m i n}^{x m a x}\Biggl(B_{j}\Bigl(x;1\Bigr)\Biggl)^{2}\nonumber++\,2\lambda\sum_{j}\Delta^{2}a_{j}\Delta^{2}a_{j-1}\int_{x m i n}^{x m a x}B_{j}\biggl(x;1\biggr)\,B_{j-1}\biggl(x;1\biggr)d x\\=\lambda\{c_{1}\sum_{j}\Bigl(\Delta^{2}a_{j}\Bigr)^{2}+c_{2}\sum_{j}\Delta^{2}a_{j}\Delta^{2}a_{j-1}\} \]

  其中\(c_1\)和\(c_2\)对于给定结点是常数，\(c_{1}=\int_{x m i n}^{x m a x}\left(b_{j}{\Big(}x;1\Big)\right\}^{2}d x\)，\(c_{2}=\int_{x m i n}^{x m a x}B_{j}{\biggl(}x;1{\biggr)}B_{j-1}{\biggl(}x;1{\biggr)}dx\)。最终得到的惩罚项里，第一项是二阶差分惩罚项，第二项是相邻二阶差分的交乘项，如果B样条在同一点重叠得过多就会使上述计算变得复杂，如果只使用差分惩罚就可以避免这个问题。

  

  文章建议使用信息准则选择\(λ\)， \[A I C\Bigl(\lambda\Bigr)=d e v(y;a,\lambda)+2d i m\left(a,\lambda\right)\]

  \(d i m\left(a,\lambda\right)\)是参数向量的有效维度，用平滑矩阵的迹作为近似；\(d e v(y;a,\lambda)\)可以使用残差平方和进行调整来估计，还要使用\(λ=0\)时的残差方差，可得 \[A I C\Bigl(\lambda\Bigr)=\sum_{i=1}^{m}\frac{\left(y_{i}-\widehat{\mu}_{i}\right)^{2}}{\widehat{\sigma}_{0}^{2}}+2t r\Bigl({H}\Bigr)-2m\ln\widehat{\sigma}_{0}-m\ln2\pi\]

  为避免使用方差，可以使用交叉验证和和广义交叉验证法。交叉验证法的最小化目标是 \[C V\Bigl(\lambda\Bigr)=\sum_{i=1}^{m}\left\{\frac{y_{i}-\widehat{y}_{i}}{1-h_{i i}}\right\}^{2}\]

  \(h_{ii}\)是平滑矩阵的对角线元素。

  广义交叉验证法的最小化目标是\[G C V\Bigl(\lambda\Bigr)=\sum_{i=1}^{m}\left\{{\frac{y_{i}-\widehat{y}_{i}}{m-{\sum_{i=1}^{m}h_{i i}}}}\right\}^{2}\]

  这两种方法算出的最优λ用于计算\({\widehat{\sigma}}_{0}\)。

  

3.4 B样条代码实现

  无论是带有二阶导平方的积分惩罚项的B样条还是P样条，我们都可以使用\(mgcv\)包中的\(gam\)命令对数据进行拟合： \[gam(formula, family=gaussian(), data=list(), weights=NULL, subset=NULL, na.action, offset=NULL, method="GCV.Cp", optimizer=c("outer","newton"), control=list(), scale=0, select=FALSE, knots=NULL, sp=NULL, min.sp=NULL, H=NULL, gamma=1, fit=TRUE, paraPen=NULL, G=NULL, in.out, drop.unused.levels=TRUE, drop.intercept=NULL, nei=NULL, discrete=FALSE,...)\]

  \(formula\)指需要将数据拟合到什么形式，可以填写平滑项\(s\)、\(te\)、\(ti\)、\(t2\)。

  其中，\(s\)用于帮助使用基于样条曲线的平滑来建立模型， \[s(..., k=-1, fx=FALSE, bs="tp", m=NA, by=NA, xt=NULL, id=NULL, sp=NULL, pc=NULL)\]

  \(k\)是使用的平滑项的基的维度，默认值是\(10\)和\(m[1]\)中较大值；\(fx\)取\(TRUE\)表明是固定自由度的回归样条，取\(False\)表明是带惩罚项的回归样条；\(bs\)指定了平滑项具体使用什么基底，比如”\(cr\)“表示具有二阶导数平方的积分惩罚项的三次样条，”\(bs\)“表示具有导数平方的积分惩罚项的B样条基，”\(ps\)“表示P样条基；\(m\)里输入了样条的阶数和惩罚的阶数，如果只输入一个数则表示样条阶数，如果多于两个数则表示加入了多种惩罚项。

  其他平滑项代表其他形式，\(te\)用于在每个边界基上设置一个与导数有关的惩罚；\(ti\)用于排除与边际平滑的\(main effects\)有关的基函数，并加上一些相互作用；\(t2\)用于张量基（\(tensor product\)）构造，导致更多惩罚。

  \(family\)用于指定数据的分布类型，\(binomial(link = "logit")，gaussian(link = "identity")，Gamma(link = "inverse")，inverse.gaussian(link = "1/mu^2")，poisson(link = "log")，quasi(link = "identity", variance = "constant")，quasibinomial(link = "logit")，quasipoisson(link = "log")\)；

  \(data\)用于指定数据；   * \(weights\)用于指定数据对似然的先验权重；

  \(method\)用于指定平滑参数估计方法，\(GCV.Cp\)指用\(GCV\)用于未知标度参数，\(Cp\)或\(UBER\)或\(AIC\)用已知标度；\(REML\)指用\(REML\)方法估计，\(ML\)指用\(ML\)方法标度；

  \(optimizer\)用于指定优化平滑参数估计标准的数值优化方法，与\(method\)有关；

  \(select\)和\(gamma\)用于增加额外的惩罚，如果\(select\)取值为\(TRUE\)且\(gamma\)取1以上的值则增加了额外的惩罚。

  因此，如果使用\(gam\)函数分别进行带二阶导平方的积分惩罚项的B样条和P样条拟合，前者应该使用\(bs=“bs”\)，通常\(s(x,bs="bs",m=c(3,2))\)为三阶B样条带二阶导数惩罚项；后者应该使用\(bs=“ps”\)，通常\(s(x,bs="ps",m=c(3,2))\)为三阶P样条带二阶差分惩罚项。具体拟合代码为：

library(mgcv)

## Loading required package: nlme

## This is mgcv 1.8-42. For overview type 'help("mgcv-package")'.

fit <- gam(y ~ s(x,bs="bs",m=c(3,2)), family = "gaussian", 
  data = data.frame(x = x, y = y), method="REML")

ggplot(data = data.frame(x = x, y = y),aes(x = x, y = y)) + geom_point() +
  geom_line(mapping = aes(x, predict(fit))) 

summary(fit)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## y ~ s(x, bs = "bs", m = c(3, 2))
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.04520    0.04607   0.981    0.329
## 
## Approximate significance of smooth terms:
##        edf Ref.df     F p-value    
## s(x) 6.222  7.106 30.64  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.687   Deviance explained = 70.6%
## -REML = 74.671  Scale est. = 0.21225   n = 100

library(mgcv)

fit <- gam(y ~ s(x,bs="ps",m=c(3,2)), family = "gaussian",
  data = data.frame(x = x, y = y), method="GCV.Cp")

ggplot(data = data.frame(x = x, y = y),aes(x = x, y = y)) + geom_point() +
  geom_line(mapping = aes(x, predict(fit))) 

summary(fit)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## y ~ s(x, bs = "ps", m = c(3, 2))
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.04520    0.04602   0.982    0.328
## 
## Approximate significance of smooth terms:
##        edf Ref.df    F p-value    
## s(x) 5.449  6.148 35.5  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.687   Deviance explained = 70.5%
## GCV = 0.22634  Scale est. = 0.21174   n = 100

  在平滑参数选择上，通过指定\(method\)为\(REML\)或者\(GCV.Cp\)，模型自动选择最优的平滑参数，我们可以使用\(summary(fit)\)查看拟合效果。以前一种情况的结果为例，描述中显示了平滑项的近似显著性，其中\(edf\)取值为6.222，拟合优度为0.687，\(REML\)值为74.671，即使用\(REML\)方法最大化似然后的最大取值为74.671。后一种情况显示\(GCV\)=0.226，为使用\(GCV\)方法最小化\(GCV(λ)\)后的最小取值，拟合优度和前一种方法很接近。从图像也可以看出，两种设置惩罚项的方法在选取不同的确定平滑参数方法后的表现仍然比较接近。

  

4 平滑参数选择

4.1 均方误差：偏差和方差的权衡

  线性回归与带惩罚项的截断幂基样条虽然本质上是不同的方法，但他们也有共同点。

  线性回归模型中最小化残差平方和之后可得拟合值为\({\hat{y}}=X{(}X^TX{)}^{-1}X^Ty=Hy\)。\(H\)被称为帽子矩阵；而带惩罚项的p阶截断幂基样条可得拟合值为\({\hat{{y}}}=X(X^{T}X+\lambda^{2p}D)^{-1}X^{T}y=S_{\lambda}y\)。\(S_λ\)可以被称为平滑矩阵。以上两个式子都可以提炼为\({\hat{{y}}}=Ly\)。即拟合相当于对\(y\)进行线性组合，满足与\(y\)无关的\(L\)被称为线性平滑类。通常情况下\(L\)不会完全地与\(y\)无关，但在此可以将基于数据的平滑参数认为是固定的，从而进行偏差和方差的讨论。

  设定一个通用的平滑模型\(y_i=f(x_i )+ε_i\), 是对\(f(x)\)的估计。那么均方误差表示为 \[M S E \{{\hat{{f}}}(x)\big\}=E[{(\widehat{f}}(x)-f(x))^2]\\=[E (\widehat{f}(x)) -f(x)]^2+E[{\widehat{f}}(x)-E(\hat f(x))]^2\\={({b i a s(\widehat{f}(x))})^{2}}+v a r \{ \hat{f}(x)\}\]

  由于现在我们对整条拟合曲线更感兴趣，因此计算平均积分平方误差 \[M S E \{ \hat f(\cdot)\}{=}\int_{\mathcal{X}}M S E\left\{\hat{{{f}}}(x)\right\}d x \]

  但为计算过程避免对\(χ\)的依赖，选择计算平均求和平方误差 \[M S S E\left\{\hat{f}(\cdot)\right\}=E\sum_{i=1}^{n}\Biggl\{\hat{f}\left(x_{i}\right)-f\left(x_{i}\right)\Biggr\}^{2}\]

  同时MSSE可以用矩阵简化表示，定义\(\mathbf{ \hat{f} }=[\hat f(x_1),…, \hat f(x_n)]\)，则 \[M S S E\left\{ \mathbf {\hat{f}}\right\}=E\;\vert\vert \hat{{{f}}}-f\;\vert\vert^2\]

  在线性平滑类的情况下，\(\mathbf{\hat{{f}}}=\mathbf L \mathbf y\)，此时平均求和平方误差可以计算： \[M S S E \{ \mathbf{\hat{f}} \} =\sum_{i=1}^{n}\Bigg\{E \hat{f}\Big(x_{i}\Big)-f\left(x_{i}\right)\Bigg\}^{2}+v a r\{\hat{f}(x_{i})\}\\=\sum_{i=1}^{n}\left\{E\left(\mathbf L \mathbf y\right)_{i}-f_{i}\right\}^2+v a r{\{}{\bigl(}\mathbf L \mathbf y{)}_{i}{\}}\\=\sum_{i=1}^{n} \{E\left(\mathbf L \mathbf y\right)_{i}-f_{i} \}^2+\{C o v\left(\mathbf L \mathbf y\right)\}_{i i}\\={\parallel}\left(\mathbf L-\mathbf I\right)\mathbf f\left|\left|^2+t r\right\{C o v\left(\mathbf L \mathbf y\right)\right\}\\=\parallel\left(\mathbf L-\mathbf I\right)\mathbf f\parallel^{2}+t r\Biggl\{\mathbf L C o v\bigl(\mathbf y\bigr)\mathbf L^{T}\Biggr\}\]

  如果假设同方差，那么 \[M S S E \{ \mathbf{\hat{f}} \}=\parallel\left(\mathbf L-\mathbf I\right)\mathbf f\parallel^{2}+\sigma_{\varepsilon}^{2}t r\{\mathbf L \mathbf L^T\}\]

  对于最小化\(MSSE\)这个目标，前一项表示偏差的平方的贡献，后一项表示方差的贡献。\(λ\)表示进行惩罚的程度，偏差总和随λ的增大而增大，方差总和随λ的增大而减小。直观来看，则如果\(λ\)是对结点数的惩罚，较大的\(λ\)将导致权重被压缩到0附近；如果\(λ\)是对导数的惩罚，较大的\(λ\)将导致图像的转折变得较为平缓，因此\(MSSE\)将存在一个理想的最小值，其对应的\(λ\)就是最优的，但需要\(f\)和\(σ_ε\)已知。如下图所示：

偏差方差权衡

  

4.2 秩与自由度

4.2.1 平滑矩阵的秩

  经典的平滑样条使用\(n\)个基函数，线性惩罚样条和三阶径向基使用\(K\)个结点因此使用\(K+2\)个基函数，远少于\(n\)个。后面的两种方法在样本量比较大的时候会节省大量计算量，那么这会有对应的成本吗？接下来将说明后两种方法丢弃了那些不重要的结点，不会影响最终拟合的结果。

  仍然考虑以\({\hat{\mathbf f}}=\mathbf L \mathbf y\)表示的广义线性平滑，\(L\)是对称且半正定的\(n×n\)矩阵，设定\(λ_1≥⋯≥λ_n\)是\(L\)的经过排序的特征值，\(v_1,…,v_n\)是对应的特征向量，\(Lv_i=λ_i v_i\)。如果对于一个\(n\)为20的均匀观测，将有20个结点的平滑样条和有9个结点的带惩罚项的三阶径向基样条的特征值分别画出来，如下图所示。

特征值

  相对应的20个结点的平滑样条的特征值图像如下图所示：

特征向量

  从特征值图像中可以看出，平滑样条的拟合主要使用了前十个特征值，三阶径向基样条的拟合只使用了前十一个特征值。从特征向量图像中可以看出，三阶径向基样条没有使用的信息是振荡得最激烈的、精度较差的一部分信息，因此这些减少了结点的平滑方法自带了省略最不重要信息的功能。   回到平滑矩阵的秩（rank）的含义，\(L\)的秩是\(L\)中线性无关的列的个数，也是\(L\)中非零特征值的个数，因此我们可以使用平滑矩阵的秩来衡量某个平滑方法的计算复杂程度。使用少于\(n\)个基函数的平滑方法称为低秩平滑，使用接近于\(n\)个基函数的平滑方法称为满秩平滑。对于数量更大的样本和涉及多种方法的平滑模型，使用低秩平滑方法节省的计算成本是巨大的。

  

4.2.2 平滑矩阵的自由度

  首先考虑针对结点数量的惩罚项，直观来说，不同的\(λ\)选择会对拟合图像产生不同影响，在结点数量相同时，较小的\(λ\)导致拟合曲线较不平滑，较大的\(λ\)导致拟合曲线趋向于基准模型，但大小适中的\(λ\)的取值将如何影响拟合曲线还不清楚。   我们将平滑矩阵的迹\(tr(S_λ)\)称为与平滑参数λ相对应的自由度，由于参数回归中帽子矩阵\(H\)的迹是待拟合参数的数量，也是自由度，而平滑矩阵\(S_λ\)是\(H\)的一般形式，因此\(tr(S_λ)\)也是自由度的一般推广，可以粗略理解为参数的等效数量，并且可以与多项式对应，一个自由度为\(ν\)的平滑估计量将数据拟合到与\((ν-1)\)阶多项式大致相同的程度。\(tr(S_λ)\)可以缩写为\(df_fit\)。

  对于一个有\(K\)个结点的\(p\)阶惩罚样条，我们有 \[d f_{f i t}=t r\{(\mathbf X^{T}\mathbf X+\lambda^{2p}\mathbf D)^{-1}\mathbf X^{T}\mathbf X\}\]

  在平滑系数\(λ\)为\(0\)时，\(tr(S_0)=p+1+K\)，而\(λ\)趋向无穷大时，\(tr(S_λ)→p+1\)，因此随\(λ\)的取值不同，\(df_fit\)在\(p+1\)至\(p+1+K\)的范围内变动。不同的\(λ\)值对应\(df_fit\)的取值之间的差异大小决定了图像外观的差异大小。

  

4.2.3 残差的自由度

  参数回归模型中残差的自由度为\(n-p\)，这一自由度用于校正使用\(RSS\)估计\(σ^2\)时的偏差。非参数估计中，定义残差自由度为 \[d f_{r e s}=n-2t r(\mathbf S_{\lambda})+t r\bigl(\mathbf {S_{\lambda}}\mathbf{ S_{\lambda}}{}^{T}\bigr)\]

  同样用于估计\(𝜎^2\)。

  设定真实模型为\(\mathbf y=\mathbf f+\mathbf{\varepsilon}\)，同样根据估计有\({\mathbf {\hat{f}}_{\lambda}}=\mathbf S_{\lambda}\mathbf y\)，于是 \[ E\left(R S S\right)=E\parallel{\hat{\mathbf f}}_{\lambda}-\mathbf y\parallel^{2}\\=E\parallel(\mathbf S_{\lambda}-\mathbf I)\mathbf y\parallel^{2}\\=E\{\mathbf y^{T}(\mathbf S_{\lambda}-\mathbf I)^{T}(\mathbf S_{\lambda}-\mathbf I)\mathbf y\}\\=\mathbf f^{T}(\mathbf S_{\lambda}-\mathbf I)^{T}(\mathbf S_{\lambda}-\mathbf I)\mathbf f+\sigma^{2}t r\{(\mathbf S_{\lambda}-\mathbf I)^{T}(\mathbf S_{\lambda}-\mathbf I)\}\\=\parallel\left(\mathbf S_{\lambda}-\mathbf I\right)\mathbf f\,\parallel^{2}+\sigma^{2}\bigl\{t r\bigl(\mathbf S_{\lambda}\mathbf S_{\lambda}{}^{T}\bigr)-2t r(\mathbf S_{\lambda})+n\bigr\}\\=\|\;(\mathbf S_{\lambda}-\mathbf I)\mathbf f\;\|^{2}+d f_{r e s}\sigma^{2} \]

  假设偏差项\(\|\,(\mathbf S_{\lambda}-\mathbf I)\mathbf f\,\|^{2}\)可以忽略不计那么\(R S S/d f_{r e s}\)是\(𝜎^2\)的无偏估计。

  且\(𝑑𝑓_{𝑓𝑖𝑡}\)和\(𝑛−𝑑𝑓_{𝑟𝑒𝑠}\)在一定程度上可替代，参数回归模型中二者相等，非参数模型中，\(𝑛−𝑑𝑓_{𝑟es}\)对于适中的\(𝜆\)取值的变动更敏感。如下图所示：

平滑矩阵自由度与残差自由度

  

4.3 经验法则

  以下部分主要针对结点数量选择，因为结点位置通常设定为等距分布或等分位数分布，位置对拟合结果影响较小。

  主要有两种经验法则，其一是在数据区间均匀分布结点，不考虑实际数据分布，只关注覆盖区间长度，其二是根据数据的分布来均匀分布结点，一般选择每两个结点有4-5个观测值，同时保持总结点数在20~40个之间。

  一个常用的具体设置方法：第\(k\)个结点的位置\(κ_k\)为样本的\((k+1)/(K+2)\)分位，结点总个数\(K\)选取样本中不重复\(x_i\)个数的四分之一与35之间的较小值。

  

4.4 CV

  考虑残差平方和 \[R S S=\textstyle\sum_{i=1}^{n}\{y_{i}-{\hat{y}}_{i}\}^{2}=|| \mathbf y-{\hat{\mathbf y}}||^{2}\]

  由于最小化残差平方和使用自身信息时会对样本点本身赋予很高权重，因此将导致在每个样本点处得到的残差平方和最小，拟合曲线会极度不光滑，相当于没有施加惩罚，因此在对\(𝑦_𝑖\)估计时将\(𝑦_𝑖\)去掉，即可以使用交叉验证法绕开这个问题。使用\({\hat{f}}(x;\lambda)\)表示在点\(x\)处的非参数估计，平滑参数为\(𝜆\)，于是残差平方和可以重写为 \[{R S S(\lambda)}={\textstyle\sum_{i=1}^{n}}{\{y_{i}-{\hat{f}}(x;\lambda)\}}^{2}\]

  交叉验证准则即为最小化 \[C V(\lambda)=\sum_{i=1}^{n}\{y_{i}-{\hat{f}}_{-i}(x;\lambda)\}^{2}\]

  下图为一个具体例子中不同\(λ\)对应的\(CV\)和\(RSS\)图像，可以看出\(λ\)趋向于0时\(RSS\)取最小值，但\(λ\)取到一个适当的中间值时\(CV\)取最小值，为平滑参数的选择提供了比较合适的参考。

CV与RSS

  由于使用\(G V(\lambda)=\textstyle\sum_{i=1}^{n}\{y_{i}-{\hat{f}}_{-i}(x;\lambda)\}^{2}\)进行计算会存在较大的计算成本，先引入一个快速算法。设\(𝑺_𝜆\)是与\(𝑓\)对应的平滑矩阵，则拟合值向量为 \[\mathbf S_{\lambda} \mathbf y=\left[\begin{array}{c}{{\hat{f}(x_{1};\lambda)}}\\ {{\vdots}}\\ {{\hat{f}(x_{n};\lambda)}}\end{array}\right]\]

  因此 \[\hat{f}(x_{i};\lambda)=\Sigma_{j=1}^{n}S_{\lambda,i j}y_{j}\]

  对于大多数平滑方法来说，\[\hat{f}_{-i}(x;\lambda)=\frac{\Sigma_{j\ne i}S_{\lambda,i j}y_{j}}{\Sigma_{j\ne i}s_{\lambda,i j}}, $\textstyle\sum_{j=1}^{n}S_{\lambda,i j}=1\]因此\[\hat{f}_{-i}(x;\lambda)=\frac{\sum_{j\neq i}S_{\lambda,i j}y_{j}}{\sum_{j\neq i}S_{\lambda,i j}}=\frac{\sum_{j\neq i}S_{\lambda,i j}y_{j}}{1-S_{\lambda,i i}}\]

  即用这个式子带入交叉验证的计算， \[C V(\lambda)=\sum_{i=1}^{n}\left\{y_{i}-\frac{\sum_{j\neq i}S_{\lambda,i j}y_{j}}{1-S_{\lambda,i i}}\right\}^{2}\\=\sum_{i=1}^{n}\left({\frac{y_{i}-{\hat{f}}(x_{i};\lambda)}{1-S_{\lambda,i}}}\right)^{2}\\=\sum_{i=1}^{n}\left({\frac{\{(\mathbf I-\mathbf S_{\lambda})\mathbf y\}_{i}}{1-S_{\lambda,ii}}}\right)^{2}\\=\sum_{i=1}^{n}\left\{{\frac{y_{i}-{\hat{y}}_{i}}{1-S_{\lambda,i i}}}\right\}\]

  此时CV可以仅用普通残差和平滑矩阵对角线元素来进行计算。

  

4.5 GCV

  在CV的快速算法还没出现时，通常对CV进行简化处理，即用\(\textstyle{\frac{1}{n}}\sum_{i=1}^{n}S_{\lambda,i i}={\frac{1}{n}}\,t r(S_{\lambda})\)代替\(𝑆_{𝜆,𝑖𝑖}\)，则广义交叉验证准则就是最小化\[G C V(\lambda)=\sum_{i=1}^{n}\bigg\{{{\frac{\{(\mathbf I-\mathbf S_{\lambda})\mathbf y\}_{i}}{1-{\frac{1}{n}}t r(\mathbf S_{\lambda})}}}\bigg\}^2=\frac{R S S(\lambda)}{\left\{1-\frac{1}{n}t r(\mathbf S_{\lambda})\right\}^{2}} \]

  广义交叉验证只是交叉验证的”近似”推广方法，在平滑样条中是由Craven and Wahba (1979)提出的。

  

4.6 信息准则

  信息准则的构造都基于对拟合优度和模型复杂程度的权衡。

  Akaike’s information criterion（AIC）表达式为 \[A I C(\lambda)=\log\{R S S(\lambda)\}+2d f_{f i t}(\lambda)/n\]

  修正AIC表达式为\[A I C_{C}(\lambda)=\log\{R S S(\lambda)\}+{\frac{2\{d f_{f i t}(\lambda)+1\}}{n-d f_{f i t}(\lambda)-2}}\]

  

4.7 似然法

  似然法的使用与混合模型紧密相关。

  首先简要介绍一下混合模型，混合模型是回归模型的拓展，允许引入随机效应\(u_i\)，允许不同个体（组别）有不同的截距项，\(u_{i}{\sim}N(0,\sigma_{u}^{2})\)，模型形式为\(y_{i j}=\alpha+u_{i}+\beta x_{i j}+\varepsilon_{i j}\)，拟合值可以写成\(\hat{\mathbf y}=\mathbf H_{0}\mathbf y+\mathbf H_{u}\mathbf y+\mathbf H_{x}\mathbf y\)，因此每个分量的自由度为 \[d f_{0}=t r(\mathbf H_{0})=1$、$d f_{x}=t r(\mathbf H_{x})=1\] \[d f_{u}\,=\,t r(\mathbf H_{u})\,=\,\frac{(m-1)n_{i}\sigma_{u}^{2}}{n_{i}\sigma_{u}^{2}+\sigma_{\varepsilon}^{2}}=\,\frac{(m-1)n_{i}}{n_{i}+\sigma_{\varepsilon}^{2}/\sigma_{u}^{2}}\]

  说明参数的等效数量取决于方差比值。这体现了这个随机截距模型是介于两个固定效应模型之间的，即只存在固定截距项的模型（删去\(u_i\)）和将每个个体的\(u_i\)赋予一个常数值（也变为固定效应）。

  类似线性模型，混合模型也可以推广到矩阵形式，\(\mathbf y=\mathbf X \mathbf{\beta}+\mathbf Z \mathbf w+\mathbf{{\varepsilon}}\)，其中\(E\left[\begin{array}{c}{{\mathbf u}}\\ {{\mathbf\varepsilon}}\end{array}\right]=\left[\begin{array}{c}{{\mathbf 0}}\\ {{\mathbf 0}}\end{array}\right]\)，\(Cov=\left[\begin{array}{c}{{\mathbf u}}\\ {{\mathbf\varepsilon}}\end{array}\right]=\left[\begin{array}{c}{{\mathbf G}}&{{\mathbf 0}}\\ {{\mathbf 0}}&{{\mathbf R}}\end{array}\right]\)，\(\mathbf G=\sigma_{u}^{2}\mathbf I\)，\(\mathbf R=\sigma_{\varepsilon}^{2}\mathbf I\)。

  要得到β的估计可以将模型重写为\(\mathbf y=\mathbf X\beta+\varepsilon^{*}\)，\(\varepsilon^{*}=\mathbf Z \mathbf u+\varepsilon\)，可得\(C o v\Big(\varepsilon^{*}\Big)=\mathbf V=\mathbf Z \mathbf G \mathbf Z^ T+\mathbf R\)，因此\(𝜷\)的\(GLS\)估计量为\({{\tilde{\beta}}}=(\mathbf X^{T}\mathbf V^{-1}\mathbf X)^{-1}\mathbf X^{T}\mathbf V^{-1}\mathbf y\)，随机效应的估计为\({{\tilde{u}}}=\mathbf G \mathbf Z^{T}\mathbf V^{-1}(\mathbf y-\mathbf X\beta)\)，可将\({{\tilde{\beta}}}\)代替\(𝜷\)求得。

  Henderson’s justification通过假设分布\(\mathbf y|\mathbf u\!\sim\!N(\mathbf X\beta+\mathbf Z \mathbf u,\mathbf R),\mathbf u\!\sim\!N(\mathbf 0,\mathbf G)\)，并通过最大化似然建立准则\((\mathbf y-\mathbf X\beta-\mathbf Z \mathbf u)^{T}\mathbf R^{-1}(\mathbf y-\mathbf X\beta-\mathbf Z \mathbf u)+\mathbf u^{T}\mathbf G^{-1}\mathbf u\)，即要得到最优线性无偏预测相当于一个带有惩罚项的广义最小二乘。

  

  关于协方差矩阵的估计主要使用极大似然法（ML）和限制性极大似然法（REML）。

  使用ML进行方差估计，其原理是对于给定的数据点，找到参数使得这部分数据点在此分布下发生的概率最大，那么这些数据点就最有可能来自这个分布，ML法对\(V\)的估计基于\(y~N(Xβ,V)\)，因此对每个样本点来说有密度函数 \[f_{i}=\frac{1}{\sigma{\sqrt{2\pi}}}e^{-\frac{(y_{i}-\alpha-\beta x_{i})^{2}}{2\sigma^{2}}}\]

  因此似然函数为 \[{ L}=\displaystyle\prod_{i=1}^{n}\frac{1}{\sigma\sqrt{2\pi}}\,e^{-\frac{(y_{i}-\alpha-\beta x_{i})^{2}}{2\sigma^{2}}}\]

  对数化后可得 \[\ln L=-\frac{n}{2}\ln{\sigma^{2}}-\frac{n}{2}\ln{2\pi}-\frac{1}{2\sigma^{2}}{\sum_{i=1}^{n}}(y_{i}-\alpha-\beta x_{i})^{2}\]

  似然函数写为矩阵形式为 \[l(\beta,\mathbf V)=-\frac{1}{2}[n l o g(2\pi)+\log\vert \mathbf V\vert+(\mathbf y-\mathbf X\beta)^{T}\mathbf V^{-1}(\mathbf y-\mathbf X\beta)\]将\({{\tilde{\beta}}}\)

  带入后可得 \[l_{p}(\mathbf V)=-\frac{1}{2}[n l o g(2\pi)+\log|\mathbf V|+\left(\mathbf y-\mathbf X\tilde{\beta}\right)^{T}\mathbf V^{-1}(\mathbf y-\mathbf X\tilde{\beta})\\=-\frac12[n l o g(2\pi)+\log|\mathbf V|+\mathbf y^{T}\mathbf V^{-1}\{\mathbf I-\mathbf X(\mathbf X^{T}\mathbf V^{-1}\mathbf X)^{-1}\mathbf X^{T}\mathbf V^{-1}\}\mathbf y\}\]

  最大化上式可以获得\(V\)的估计。同时求偏导后求期望可得\(E(\hatσ^2)=σ^2-σ^2/n\)，因此在高斯分布下极大似然法对方差的估计是有偏的，因此一般使用REML估计\(\hat σ^2\)并选择惩罚程度\(λ\)。

  REML相对于ML的主要优势在于考虑了固定效应（\(β\)）的自由度，纠正了ML在方差估计中的偏误，将ML的基础\(\mathbf y\sim N\left(\mathbf X\beta,\mathbf V\right)\)变换为\(A^{\prime}Y\sim N(0,\sigma^{2}A^{\prime}A)\)，从而转变为对\(A'Y\)进行极大似然估计，以校正偏误。REML的似然函数为 \[l_{R}(\mathbf V)=l_{p}(\mathbf V)-\log|\mathbf X^{T}\mathbf V^{-1}\mathbf X|\]

  对于惩罚样条，可将\(\mathbf y=\mathbf X\beta+\mathbf Z \mathbf u+\varepsilon\)中的\(𝒁\)看作样条项，\(𝒖\)即为样条项系数，那么惩罚样条的最小化目标可以写为（除以\(𝜎_𝜀^2\)后） \[\textstyle{\frac{1}{\sigma^{2}}}\parallel \mathbf y-\mathbf X\beta-\mathbf Z \mathbf u\parallel^{2}+{\frac{\lambda^{2}}{\sigma^{2}}}\parallel \mathbf u\parallel^{2}\]

  如果将𝒖看作一组协方差为\(\sigma_{u}^{2}\mathbf I\)且\(\sigma_{u}^{2}=\sigma_{\varepsilon}^{2}/\lambda^{2}\)的随机系数，就与之前提到的Henderson’s justification的表达形式相似，以上即可得到惩罚样条的混合模型表示。因此如果最小化了REML得到了\((\hat{\sigma}_{\varepsilon,R E M L}^{2},\hat{\sigma}_{u,R E M L}^{2})\)，可以得到平滑参数为\(\hat{\lambda}_{R E M L}=\left(\hat{\sigma}_{\varepsilon,R E M L}^{2}/\hat{\sigma}_{u,R E M L}^{2}\right)^{1/2p}\)。

  基于最大似然思想的平滑参数选择方法依赖于惩罚样条的混合模型表示，因此基于经典模型选择思想的更通用的方法在现实中使用得更多。

5 平滑样条

  平滑样条的思路与前两大类样条不太相同。通过使用所有的样本点作为结点，因而避免了结点选择问题，并通过添加惩罚项来压缩系数以避免模型过于复杂，因此平滑样条考虑的是以下问题：在所有具有连续二阶导的函数中找到一个使下列带惩罚项的残差平方和最小的函数\[\sum_{i=1}^{n}{\bigl(}y_{i}-f(x_{i}){\bigr)}^{2}+\lambda\int f^{\prime\prime}(x)^{2}d x\\=\sum_{i=1}^{n}{\bigl(}y_{i}-f(x_{i}){\bigr)}^{2}+\lambda\int\left(\sum_{i}\beta_{j}f_{j}{}^{\prime\prime}(x)\right)^{2}d x\\=\sum_{i=1}^{n}{\bigl(}y_{i}-f(x_{i}){\bigr)}^{2}+\lambda\sum_{i,j}\beta_{i}\beta_{j}\int f_{i}^{\prime\prime}(x)f_{j}^{\prime\prime}(x)d x\\=\|\mathbf y-\mathbf F\beta\ \|_{2}^{2}+\lambda\beta^{T}\Omega\beta\]

  这个最小化目标定义在无限维的函数空间上，使用λ在拟合优度和模型复杂度之间建立了权衡，\(λ=0\)意味着没有惩罚，\(λ=∞\)意味着模型被压缩到线性模型，因为不允许任何二阶导数的存在。这个问题有一个唯一的、有限维的解，是结点位于每个样本点处的三阶自然样条。因此平滑样条的使用思路与三阶自然样条类似，同时由于三阶自然样条与三阶B样条、三阶截断幂样条都等价，他们的使用思路也都类似。

  

  平滑样条可以通过\(stats\)包中的\(smooth.spline\)函数进行使用： \[smooth.spline(x, y = NULL, w = NULL, df, spar = NULL, lambda = NULL, cv = FALSE, all.knots = FALSE, nknots = .nknots.smspl, keep.data = TRUE, df.offset = 0, penalty = 1, control.spar = list(), tol = 1e-6 * IQR(x), keep.stuff = FALSE)\]

   - \(w\)表示可对\(x\)设定权重；

   - \(df\)表示所需的等效自由度，即平滑矩阵的迹，必须小于不重复的\(x\)的个数；

   - \(spar\)和\(lambda\)用于设定平滑参数；

   - \(cv\)取值为\(TRUE\)则表明使用经典的去一法选择平滑参数，取值为\(FALSE\)则表明使用\(GCV\)选择平滑参数；

   - \(all.knots\)取值为\(TRUE\)则所有\(x\)都被用作结点，取值为\(FALSE\)则使用\(nknots\)设定结点。

  具体拟合代码为：

library(splines)

# 1
fit1 <- smooth.spline(x, y, df=4)

ggplot(data = data.frame(x = x, y = y), aes(x, y)) + 
  geom_point() + geom_line(mapping = aes(x, fit1$y))

# 2
fit2 <- smooth.spline(x, y, cv=TRUE, all.knots=TRUE)

ggplot(data = data.frame(x = x, y = y), aes(x, y)) + 
  geom_point() + geom_line(mapping = aes(x, fit2$y))

# 3
fit3 <- smooth.spline(x, y, cv=FALSE, all.knots=TRUE)

ggplot(data = data.frame(x = x, y = y), aes(x, y)) + 
  geom_point() + geom_line(mapping = aes(x, fit3$y))

fit1

## Call:
## smooth.spline(x = x, y = y, df = 4)
## 
## Smoothing Parameter  spar= 0.9174261  lambda= 0.01987269 (11 iterations)
## Equivalent Degrees of Freedom (Df): 3.999631
## Penalized Criterion (RSS): 29.35932
## GCV: 0.3185665

fit2

## Call:
## smooth.spline(x = x, y = y, cv = TRUE, all.knots = TRUE)
## 
## Smoothing Parameter  spar= 0.8392485  lambda= 0.0008675094 (12 iterations)
## Equivalent Degrees of Freedom (Df): 7.561031
## Penalized Criterion (RSS): 19.42816
## PRESS(l.o.o. CV): 0.2300711

fit3

## Call:
## smooth.spline(x = x, y = y, cv = FALSE, all.knots = TRUE)
## 
## Smoothing Parameter  spar= 0.8344394  lambda= 0.000800811 (13 iterations)
## Equivalent Degrees of Freedom (Df): 7.693516
## Penalized Criterion (RSS): 19.37156
## GCV: 0.2273527

  同样我们可以输入\(fit\)命令查看拟合效果描述，以上分别是上述三种稍有差异的拟合的效果，可以看到当我们指定自由度为4时，等效自由度其实为3.999，\(GCV\)为0.318，而输入使用所有样本点都作为结点时，第二种方法显示等效自由度为7.561，\(CV\)值为0.230，第三种方法显示等效自由度为7.693，\(GCV\)值为0.227，后两种方法的等效自由度非常接近，而第一种方法的\(GCV\)值高于第三种方法，意味着后两种方法的拟合效果优于第一种。

  

6 多维情形

6.1 径向基

  径向基（Radial Basis Function）最初由Michael J. D. Powell于1977年提出，用于求解偏微分方程和数值逼近问题，后被David Broomhead和David Lowe应用于机器学习，后常用于计算机图像、神经网络等领域。

  任意一个满足\(\Phi{(x)}=\Phi(\lVert x \rVert)\)特性的函数\(\Phi\)都叫做径向基函数，因此具体形式较为丰富。径向基函数的取值只取决于与固定点之间的距离，这个距离通常指欧式距离。由于按点移位且在每个方向上径向对称，径向基函数最大的优点之一就是适用于高维，数据除了需要位于不同的点之外没有受到任何限制，且在数据变得密集时径向基仍然能保持高精度和快速收敛。

  一阶径向基函数以绝对值为基础进行定义，图像如下图所示。\(p\)是奇数时的径向基函数与截断幂基函数和B样条基函数等价。

径向基

  使用径向基函数获得的拟合函数形式为： \[{\hat{y}}(x)=\textstyle{\sum_{i=1}^{n}w_{i}\varphi(\lVert x-x_{i}\rVert)}\]

  即将以不同的结点为中心的基函数以\(w_i\)为权重相加。另有一些常用变体形式，例如高斯形式\(\varphi(x)=e^{-(\varepsilon\|x-x_{i}\|)^{2}}\)

  如果将径向基函数推广到高维，即将\(x\)换为多维输出向量，\(x_i\)换为多维中心向量。高维数据中数据点之间的距离变得更为稀疏，但过大的宽度参数\(ε\)又会使模型欠拟合，因此选择合适的宽度参数是非常关键的，可以使用交叉验证方法进行宽度选择。另外为防止过拟合问题，也可以加入惩罚项进行限制。

  中心向量或者中心点的选取可以选取使用\(K-means\)、层次聚类等聚类算法对输入数据进行聚类，从每个类别中选取一个代表点作为中心结点，这个代表点可以是类别的中心点，也可以是类别中距离其他样本最远的点，聚类方法在保留数据特征的同时减少了结点数量。中心结点还可以使用基于核函数的谱分解方法选择，在定义好核函数之后计算数据的核矩阵，即将样本点对应的特征向量带入核函数计算相似度矩阵，接着对和矩阵进行谱分解得到特征值和特征向量，选择前若干个最大的特征值对应的特征向量作为中心结点，这种方法确定的中心结点更加精确，但计算量较大。其余方法还有基于格点采样、基于最大边缘采样等。

  

6.2 张量积

  张量积样条函数为一维样条函数提供了一个向高维推广的简单明了的方式。

  张量积即为外积，假设我们有一个n维向量V，和一个m维向量W，张量积\(V⊗W\)定义了\(m*n\)维的向量空间。张量积在量子物理和深度学习的张量乘法中有广泛的应用，仅关注数据拟合时，对应于前文的一维样条适用于将散点拟合成曲线，张量积适用于将多维数据拟合成超平面，例如使用张量积可以定义一个光滑的函数来表示相互作用的协变量，且这样的基底对所有协变量的重缩放（rescaling）都不变。

  例如对与两个平滑函数 \[f_{x}(x)=\sum_{i=1}^{I}\alpha_{i}a_{i}(x)\] \[f_{z}(z)=\sum_{l=1}^{L}\delta_{l}d_{l}(z)\]

  如果要求\(f(x)\)随\(z\)平滑地变化，则可以通过允许\(f(x)\)的参数随\(z\)平滑变化来实现， \[\alpha_{i}(z)=\sum_{l=1}^{L}\delta_{i l}d_{l}(z)\]

  现在可以在\(x,z\)上定义一个光滑函数 \[f_{x z}(x,z)=\sum_{i=1}^{I}\sum_{i=1}^{L}\delta_{i l}d_{i}(z)a_{i}(x)\]

  可以改写成以下形式\(f_{x z}(x,z)=a(x)^{T}\delta d(z)\)

  以上的例子展示了两个光滑的一维函数\(f(x)\)和\(f(z)\)的张量积可以张成一个二维子空间。

  

7 参考资料

1. Ruppert, D., Wand, M. P., & Carroll, R. J. (2003). Semiparametric regression (No. 12). Cambridge university press.
2. Hastie, T., Tibshirani, R., Friedman, J. H., & Friedman, J. H. (2009). The elements of statistical learning: data mining, inference, and prediction (Vol. 2, pp. 1-758). New York: springer.
3. https://bookdown.org/ssjackson300/Machine-Learning-Lecture-Notes/
4. 正交多项式：https://www.sydney.edu.au/science/chemistry/\~mjtj/CHEM3117/Resources/poly_etc.pdf
5. 三角多项式： https://people.math.carleton.ca/~ckfong/C4.pdf
6. 伯恩施坦多项式：https://web.mit.edu/hyperbook/Patrikalakis-Maekawa-Cho/node8.html
7. B样条的惩罚项:

• O’Sullivan, F. (1986). A statistical perspective on ill-posed inverse problems. Statistical science, 502-518.
• O’Sullivan, F. (1988). Fast computation of fully automated log-density and log-hazard estimators. SIAM Journal on scientific and statistical computing, 9(2), 363-379.
• Wand, M. P., & Ormerod, J. (2008). On semiparametric regression with O’Sullivan penalized splines. Australian & New Zealand Journal of Statistics, 50(2), 179-198.
• https://search.r-project.org/CRAN/refmans/mgcv/html/smooth.construct.bs.smooth.spec.html

8. P样条:

• Eilers, P. H., & Marx, B. D. (1996). Flexible smoothing with B-splines and penalties. Statistical science, 11(2), 89-121.
• Perperoglou, A., Sauerbrei, W., Abrahamowicz, M., & Schmid, M. (2019). A review of spline function procedures in R. BMC medical research methodology, 19(1), 1-16.
• https://search.r-project.org/CRAN/refmans/mgcv/html/smooth.construct.ps.smooth.spec.html

9. 径向基: http://www.scholarpedia.org/article/Radial_basis_function