b样条曲面绘制 opengl_怎样用R构建样条回归模型

1 前言

在之前的文章中(非线性回归之多项式回归(R语言版))我们介绍了在R中如何构建多项式回归模型。但是使用多项式回归模型有一个限制就是，所有的原始数据都会用到。而在实际当中，有可能会遇到自变量在某些点之前和之后发生突变的情况，如政策的变化、战争、革命和地震等突发事件的发生往往会引起函数关系的变化。这时候，若只是简简单单的用一个直线函数或多项式函数来拟合往往不准确。那我们是是否能够构建一个函数，使得函数在不同的区间段内能够很好的拟合这个区间的数据，然后曲线在整个区间上是连续的呢？答案是有的。根据函数是否处处可微分为样条函数和阶跃函数。这里介绍的是样条函数。利用该函数构造样条回归模型。

2 预备数学知识

2.1 引入

多项式函数的形式如下：

写成求和的形式：

多项式函数的形式可以认为是由多个单项式构成{

}。现在，我们可以把每一项都换成一个关于自变量x的函数

:

那么样条函数就是

其中

其中d为节点的个数，q表示多项式的最高次数。

可能说的有点抽象。简单来说，就是将之前的多项式中的每一项换成一个新的多项式。而节点个数，就是突变的点个数(分段点)，比如，像下图所示，节点就为2个（d = 2）

函数的值，y不仅与x有关，还与要拟合的多个多项式的最高次数q有关，准确来说，应该写成

表示在最高次数q确定的情况下，拟合的值。

2.2 需要注意的地方

a.由于函数是连续的，因此要满足:

以上图为例，第一段直线在1处的值和第二段之前在1处的值要相等。不仅如此，还要满足一阶导数和二阶导数相等，这样就可以保证曲线光滑连续。

b. 对于曲线拟合，d一般选择3，通常称为立方回归样条曲线(贝塞尔曲线)。

c.如何计算自由度。自由度就是回归模型中待求参数的个数，由于有d个节点,那么就将一条曲线分成了d+1段，然后每一段为一个q次多项式，q次多项式就有q+1个参数。那么待求参数的个数是不是就是(d+1)(q+1)了？答案并不是原因是整个函数都是连续的，所以根据2,2 中的a可以得到每段之间的关系，这样的话，原来两段的函数有(d+1)(q+1)个参数，现在由于有等式的存在，只需要(d+1)*（q+1）-mq个参数的值，化简就为m+q+1个参数。

3 R语言实现

这里选用的数据集是R中自带的数据集。为MASS包中mcycle数据集。它测试了一系列模拟的交通车事故中，头部的加速度，以此来评估头盔的性能。times为撞击时间(ms)，accel为加速度（g）。首先导入数据，并绘制散点图

> library(MASS)
> x <- mcycle$times
> y <- mcycle$accel
> plot(x,y,xlab = "times / ms",ylab = "acceleration / g")

散点图如下：

看图，可以设置x = 15、20、32、40为节点，这样就有4个节点，然后使用d = 3的样条曲线拟合（也可以设置为4，常用d = 3）

在R中，需要splines包中的bs函数来绘制样条曲线。bs函数的用法如下：

bs(x, df = NULL, knots = NULL, degree = 3, intercept = FALSE,
   Boundary.knots = range(x))

参数解释：

x：自变量，这里为x

df: 自由度（回归函数中的参数个数）,默认为0，由于我们需要截距,2.2 节中c中提到的公式减去1，因此自由度为4+3 = 7

knots：节点的位置，这里为c(15,20,32,40)

degree：q值，默认为3

其他的参数保持默认即可

然后搭配lm 函数，即可画出样条函数曲线

> bspl <- lm(y~bs(x,df = 7,knots = c(15,18,32,40))
> plot(x,y)
> lines(x,fitted(bspl),lwd = 2, col = 2)

最后的拟合图如下：

可以看出，拟合效果很好。

好了，关于样条回归模型的R语言构建就到这里，大家如果有什么问题，欢迎在评论区留言(*^▽^*)

​