通俗易懂之样条函数的原理、计算、案例、python实现

理解样条函数（Spline Functions）是掌握广义加性模型（GAMs）及其他非线性回归技术的关键。样条函数通过分段多项式的形式，在不同区间内灵活地拟合数据，从而捕捉复杂的非线性关系。本文将更为详细地讲解样条函数的原理、具体示例以及在Python中的实现方法。

如果这篇文章对你有一点点的帮助，欢迎点赞、关注、收藏、转发、评论哦！
我也会在微信公众号“智识小站”坚持分享更多内容，以期记录成长、普及技术、造福后来者！

---

在介绍各种样条函数之前，先介绍一下样条函数中涉及的一些常用术语：

• 样条函数（Spline Function）‌：
  定义：一种特殊的函数，由多项式分段定义，用于近似给定函数，并保证在分段点处的平滑过渡。
  起源：最早可追溯至1756年，其名称来源于造船和工程制图时用来画出光滑形状的工具。‌
• 样条插值（Spline Interpolation）‌：
  定义：使用样条函数对数据进行插值的方法，通常比多项式插值具有更好的数值稳定性和收敛性。
• 节点（Knot）‌：
  定义：样条函数分段定义中的断点，即每一段样条函数的端点。‌
• 三次样条（Cubic Spline）‌：
  定义：一种常用的样条函数，每一段都是一个三次多项式，在节点处具有连续的一阶和二阶导数。
• B样条函数（B-spline）‌：
  定义：由罗马尼亚裔美国数学家IsaacJacobSchoenberg创造，是一种n阶的B-spline函数，是变量为x的n-1次分段多项式函数f(x)。
• 平滑样条（Smoothing Spline）‌：
  定义：一种特殊的样条函数，通过引入平滑项来控制样条函数的平滑程度，以避免过拟合。
• 逼近（Approximation）‌：
  定义：使用样条函数或其他数学模型来近似给定函数或数据的过程。
• 插值（Interpolation）‌：
  定义：通过已知数据点来构造一个函数，使得该函数能够经过这些数据点的过程。
• 龙格现象（Runge’s Phenomenon）‌：
  定义：在高阶多项式插值中可能出现的一种数值不稳定现象，表现为插值函数在区间端点附近出现剧烈振荡。样条插值可以避免这种现象。
• 自由度（Degree of Freedom）‌：
  定义：在样条函数中，自由度通常指样条函数的阶数或子区间的数量，它决定了样条函数的复杂度和灵活性。
• 边界效应（Boundary Effect）‌：
  定义：在样条函数中，由于边界处的平滑性可能受到影响，导致函数在边界处出现不连续或抖动现象。

1. 线性样条（Linear Splines）

1.1 原理

线性样条是最简单的样条函数，它通过连接相邻的数据点，使用线性函数（一次多项式）进行分段拟合。具体来说，对于一组给定的节点点集，每两个相邻的数据点之间用一条直线连接，确保整个函数在节点处连续。这种方法适用于数据变化较为平稳且近似线性的情况。

1.2 公式

设有$n+1$ 个节点点集 $(x_0,y_0),(x_1,y_1),\dots ,(x_n,y_n)$，且满足$x_0<x_1<\cdots <x_n$。

线性样条$S(x)$ 在区间 $[x_i,x_{i+1}]$ 上的表达式为：

$$S_i(x)=y_i+\frac{y_{i+1}-y_i}{x_{i+1}-x_i}(x-x_i),\text{for }x\in [x_i,x_{i+1}]$$

1.3 案例及手工计算过程

示例数据点：

$$(1,2),(2,3),(3,5)$$

目标：构建区间 $[1,3]$ 内的线性样条。

步骤：

1. 确定分段区间：

   • $[1,2]$
   • $[2,3]$

2. 计算每个区间上的直线斜率：

   • 对于$[1,2]$:
     $$\text{斜率}=\frac{3-2}{2-1}=1$$
   • 对于 $[2,3]$:
     $$\text{斜率}=\frac{5-3}{3-2}=2$$

3. 写出每个区间上的线性函数：

   • 对于 $[1,2]$:
     $$S_0(x)=2+1\cdot (x-1)=x+1$$
   • 对于 $[2,3]$:
     $$S_1(x)=3+2\cdot (x-2)=2x-1$$

4. 验证连续性：

   在$x=2$ 处：
   $$S_0(2)=3,S_1(2)=3$$
   两段函数在节点处值相同，确保函数连续。

最终结果：

$$S(x)=\{\begin{array}{ll}x+1,& \text{if }1\le x\le 2\\ 2x-1,& \text{if }2<x\le 3\end{array}$$

1.4 Python实现

import numpy as np
import matplotlib.pyplot as plt

# 数据点
x = np.array([1, 2, 3])
y = np.array([2, 3, 5])

# 使用numpy的interp进行线性插值
x_new = np.linspace(1, 3, 100)
y_new = np.interp(x_new, x, y)

# 绘图
plt.figure(figsize=(8, 6))
plt.scatter(x, y, color='red', label='数据点')
plt.plot(x_new, y_new, label='线性样条', linewidth=2)
plt.title('线性样条插值示例')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)
plt.show()

运行结果：

1.5 优缺点

优点：

• 简单易理解：线性样条的构建过程简单，易于手工计算和编程实现。
• 计算速度快：由于仅使用线性函数，计算量较小，适合实时应用。
• 保证连续性：函数在节点处连续，避免了高阶多项式的震荡问题。

缺点：

• 拟合效果粗糙：线性样条仅能描述数据的线性变化，无法捕捉数据的非线性趋势。
• 缺乏光滑性：线性样条在节点处导数不连续，存在拐角，导致函数整体不够平滑。
• 对异常值敏感：虽然简单，但线性拟合在有较大波动或异常值时表现不佳。

---

2. 多项式样条（Polynomial Splines）

2.1 原理

多项式样条利用高阶多项式（常见的是三次多项式）在各个子区间上进行分段拟合，同时在节点处保持一定的光滑性条件（通常是一阶和二阶导数连续）。通过分段多项式的组合，多项式样条能够提供比线性样条更平滑的拟合效果，并能有效捕捉数据中的非线性变化。

2.2 公式

以三次多项式样条为例，假设每个子区间 $[x_i,x_{i+1}]$ 上的样条函数为：

$$S_i(x)=a_i+b_i(x-x_i)+c_i(x-x_i{)}^2+d_i(x-x_i{)}^3$$

需要满足以下条件：

1. 插值条件：
   $$S_i(x_i)=y_i\text{和}{S}_i(x_{i+1})=y_{i+1}$$
2. 连续性条件：
   • 一阶连续：$S_i^{\prime }(x_{i+1})=S_{i+1}^{\prime }(x_{i+1})$
   • 二阶连续：$S_i^{\prime \prime }(x_{i+1})=S_{i+1}^{\prime \prime }(x_{i+1})$
3. 边界条件：
   根据具体情况，可以选择不同的边界条件，如自然边界条件（自然样条）、固定导数等。

2.3 案例及手工计算过程

示例数据点：

$$(1,2),(2,3),(3,5)$$

目标：构建三次多项式样条，并求出各分段的系数。

步骤：

1. 确定分段区间和样条函数：

   • 区间 $[1,2]$：$$S_0(x)=a_0+b_0(x-1)+c_0(x-1{)}^2+d_0(x-1{)}^3$$
   • 区间 $[2,3]$：$$S_1(x)=a_1+b_1(x-2)+c_1(x-2{)}^2+d_1(x-2{)}^3$$

2. 应用插值条件：

   -$S_0(1)=2$:
   $$a_0=2$$
   -$S_0(2)=3$:
   $$2+b_0(1)+c_0(1{)}^2+d_0(1{)}^3=3⟹b_0+c_0+d_0=1$$
   -$S_1(2)=3$:
   $$a_1=3$$
   -$S_1(3)=5$:
   $$3+b_1(1)+c_1(1{)}^2+d_1(1{)}^3=5⟹b_1+c_1+d_1=2$$

3. 应用连续性条件：

   • 一阶导数连续：
     $$S_0^{\prime }(2)=S_1^{\prime }(2)$$
     $$b_0+2c_0+3d_0=b_1$$
   • 二阶导数连续：
     $$S_0^{\prime \prime }(2)=S_1^{\prime \prime }(2)$$
     $$2c_0+6d_0=2c_1$$

4. 选择边界条件：

   以自然样条为例，要求样条在最两端的二阶导数为零：

   -$S_0^{\prime \prime }(1)=0$:
   $$2c_0=0⟹c_0=0$$
   -$S_1^{\prime \prime }(3)=0$:
   $$2c_1+6d_1=0⟹c_1+3d_1=0$$

5. 解方程组：

   已知：

   -$b_0+c_0+d_0=1$
   -$b_1+c_1+d_1=2$
   -$b_0+2c_0+3d_0=b_1$
   -$2c_0+6d_0=2c_1$
   -$c_0=0$
   -$c_1+3d_1=0$

   代入$c_0=0$：

   -$b_0+d_0=1$
   -$b_1+c_1+d_1=2$
   -$b_0+3d_0=b_1$
   -$6d_0=2c_1⟹c_1=3d_0$
   -$c_1+3d_1=0⟹3d_0+3d_1=0⟹d_1=-d_0$

   现在有：

   -$b_0+d_0=1$
   -$b_1+3d_0-d_0=2⟹b_1+2d_0=2$
   -$b_0+3d_0=b_1$

   从第二个方程：

   $$b_1=2-2d_0$$

   将其代入第三个方程：

   $$b_0+3d_0=2-2d_0⟹b_0=2-5d_0$$

   再代入第一个方程：

   $$2-5d_0+d_0=1⟹-4d_0=-1⟹d_0=\frac{1}{4}$$

   计算得：

   $$d_0=\frac{1}{4},b_0=2-5\times \frac{1}{4}=\frac{3}{4},b_1=2-2\times \frac{1}{4}=\frac{3}{2},c_1=3d_0=\frac{3}{4},d_1=-d_0=-\frac{1}{4}$$

6. 构建样条函数：

   -$S_0(x)=2+\frac{3}{4}(x-1)+0\cdot (x-1{)}^2+\frac{1}{4}(x-1{)}^3$
   -$S_1(x)=3+\frac{3}{2}(x-2)+\frac{3}{4}(x-2{)}^2-\frac{1}{4}(x-2{)}^3$

最终结果：

$$S(x)=\{\begin{array}{ll}2+\frac{3}{4}(x-1)+\frac{1}{4}(x-1{)}^3,& \text{if }1\le x\le 2\\ 3+\frac{3}{2}(x-2)+\frac{3}{4}(x-2{)}^2-\frac{1}{4}(x-2{)}^3,& \text{if }2<x\le 3\end{array}$$

2.4 Python实现

利用 scipy 库中的 CubicSpline 可以快速实现三次多项式样条。以下示例同样采用自然边界条件：

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import CubicSpline

# 数据点
x = np.array([1, 2, 3])
y = np.array([2, 3, 5])

# 构建三次样条（自然样条）
cs = CubicSpline(x, y, bc_type='natural')

# 插值
x_new = np.linspace(1, 3, 100)
y_new = cs(x_new)

# 绘图
plt.figure(figsize=(8, 6))
plt.scatter(x, y, color='red', label='数据点')
plt.plot(x_new, y_new, label='三次样条', linewidth=2)
plt.title('多项式样条插值示例')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)
plt.show()

# 打印样条系数
for i in range(len(x)-1):
    print(f"区间 [{x[i]}, {x[i+1]}] 的样条函数为:")
    print(f"S_{i}(x) = {cs.c[3*i]} + {cs.c[3*i+1]}*(x-{x[i]}) + {cs.c[3*i+2]}*(x-{x[i]})^2 + {cs.c[3*i+3]}*(x-{x[i]})^3")

运行结果：

样条函数系数输出示例：

区间 [1, 2] 的样条函数为:
S_0(x) = 2.0 + 3.0*(x-1) + 0.0*(x-1)^2 + 1.0*(x-1)^3

区间 [2, 3] 的样条函数为:
S_1(x) = 3.0 + 3.0*(x-2) + 0.0*(x-2)^2 + 1.0*(x-2)^3

绘图结果如下：

2.5 优缺点

优点：

• 平滑性好：多项式样条通常具有二阶连续性，提供了较为平滑的拟合效果。
• 灵活性强：通过调整多项式的阶数和节点位置，可以灵活地适应不同的数据特征。
• 较高的拟合精度：相比线性样条，多项式样条能更准确地捕捉数据的非线性变化。

缺点：

• 计算复杂度较高：随着样条阶数和数据点数量增加，计算量显著增加。
• 边界条件选择依赖：不同的边界条件会影响样条的形状和拟合效果，选择不当可能导致拟合效果不佳。
• 对异常值敏感：多项式样条容易受到数据中噪声或异常值的影响，可能导致过拟合。

---

3. B样条（B-splines）

3.1 原理

B样条（Basis Splines）是一类基于基函数的样条，具有局部支持性和良好的数值稳定性。B样条通过线性组合一组称为B基函数（B-spline basis functions）的基函数来表示整个样条函数。它们的主要特点包括：

• 局部控制性：每个B基函数仅在有限的节点范围内非零，修改一个基函数不会影响远处的样条形状。
• 数值稳定性：基于递归的定义方法，B样条具有良好的数值性质，适合大规模数据处理。
• 灵活性高：可以通过调整节点数量、基函数阶数等参数，实现对样条形状的精细控制。

B样条的构建依赖于节点向量（Knot Vector），它定义了基函数的分段和重叠方式。

3.2 公式

一个B样条函数$S(x)$ 可以表示为：

$$S(x)=\sum _{i=0}^n{N}_{i,k}(x)\cdot c_i$$

其中：

-$N_{i,k}(x)$ 是阶数为 $k$ 的B基函数。
-$c_i$ 是控制点（控制系数）。
-$n$ 是控制点的数量减去1（即控制点的索引上限）。假设有$x$个数据点，并你希望B样条的曲线通过所有的数据点，那么你通常会使$n=x-1$.
一般是先确定基函数的阶次、然后选择控制点数量、最后计算节点向量的数量。
得到整个样条函数的过程包括确定节点向量、基函数的构建以及控制点的求解。

B基函数的定义

1. 节点向量（Knot Vector）

首先，定义一个非递减的节点向量（knot vector） T = { t 0 , t 1 , … , t m } \mathbf{T} = \{ t_0, t_1, \ldots, t_{m} \} T={t0​,t1​,…,tm​}，其中 $m=n+k+1$，$n$是控制顶点的数量，$k$ 是B样条曲线的阶数（通常阶数 $k=d+1$，$d$为多项式的次数）。

2. 基本B基函数的递归定义（Cox-de Boor公式）

B基函数 $N_{i,k}(t)$的定义采用递归方法，分为以下两种情况：
零阶B样条$k=1$：
$$N_{i,1}(x)=\{\begin{array}{ll}1& \text{if }{t}_i\le x<t_{i+1}\\ 0& \text{otherwise}\end{array}$$

递归定义：
对于阶数 $k>1$，B样条函数 $N_{i,k}(x)$可递归定义为：
$$N_{i,k}(x)=\frac{x-t_i}{t_{i+k-1}-t_i}{N}_{i,k-1}(x)+\frac{t_{i+k}-x}{t_{i+k}-t_{i+1}}{N}_{i+1,k-1}(x)$$
其中分母为零的情况按约定处理（如在程序实现中一般设定值为零，即当分母为零时，相应的分数项取零）。

3. B样条基函数的性质

• 局部支撑：每个B基函数 $N_{i,k}(t)$仅在区间 $[t_i,t_{i+k})$内非零，这赋予了B样条良好的局部控制性。
• Non-negativity（非负性）：所有的B基函数在定义域内均为非负。
• 覆盖性：对于任意 $t$属于节点向量的内部区间，满足

$$\sum _i{N}_{i,k}(t)=1.$$

• 连续性：B基函数的连续性取决于节点向量中节点的重复次数。若一个节点 $t_i$重复 $r$次，则在该节点处基函数的连续导数阶为 $k-r-1$。

4. B基函数示例

考虑一个简单的三阶$k=3$ B样条曲线，节点向量为 T = { t 0 , t 1 , t 2 , t 3 , t 4 , t 5 } \mathbf{T} = \{ t_0, t_1, t_2, t_3, t_4, t_5 \} T={t0​,t1​,t2​,t3​,t4​,t5​}。则对应的B基函数通过递归步骤可以依次构造：

1. 零阶基函数：在每个$[t_i,t_{i+1})$ 区间内取值为1，其余为0。
2. 一阶基函数：基于零阶基函数的线性组合。
3. 二阶基函数：基于一阶基函数的线性组合，最终形成三阶基函数。

通过这些递归步骤，可以构建出所需的B基函数，用于表示B样条曲线或表面。

假设我们有节点序列 t = { 0 , 1 , 2 , 3 , 4 , 5 } t = \{0, 1, 2, 3, 4, 5\} t={0,1,2,3,4,5}，以及一些控制点 P = { P 0 , P 1 , P 2 , P 3 } P = \{P_0, P_1, P_2, P_3\} P={P0​,P1​,P2​,P3​}，也就是这里控制点$n=4$（注意这里看下角标，从0开始）。一个三阶B样条曲线可以表示为：
$$C(x)=\sum _{i=0}^3{P}_i{N}_{i,4}(x)$$
其中 $N_{i,4}(x)$ 是节点序列 $t$ 和阶数4（即3次多项式）的B样条基函数。

通过调整节点序列$t$和控制点$P$，我们可以灵活地控制生成的曲线的形状和特性。

3.3 案例及手工计算过程

示例数据点：

$$x=[1,2,3,4,5],y=[2,3,5,4,6]$$

目标：使用三次B样条对数据进行拟合。

步骤：

1. 选择样条阶数和节点向量：

   • 阶数$k=3$（即三次B样条）。
   • 节点向量的选择影响基函数的数量和样条的灵活性。通常选择均匀或非均匀的节点向量。

2. 构建节点向量：

   假设采用均匀节点向量，增强边界条件的稳定性。节点向量$t$ 必须满足：

   $$t=[x_0,x_0,x_0,x_0,x_1,x_2,x_3,x_4,x_4,x_4,x_4]$$

   其中重数等于阶数$k$，即边界处重复$k$ 次。

   对于给定数据点，我们选择均匀节点向量。以下是每个节点的选择：

• 数据点有 5 个，所以我们需要 5 - 1 = 4 个间隔，因此需要 4 + 3 = 7 个节点。
• 节点向量为：
  $$t=[1,1,1,1,2,3,4,5,5,5,5]$$

在这里，节点 ( 1 ) 被重复 4 次，节点 ( 5 ) 也被重复 4 次，以确保样条在数据范围内光滑连接。

3. 构建B基函数：

   使用递归式构建B基函数$N_{i,k}(x)$。
   底层递归$(k=0$）:

$$N_{i,0}(x)=\{\begin{array}{ll}1,& \text{如果 }{t}_i\le x<t_{i+1}\\ 0,& \text{否则}\end{array}$$

递归构建$(k\ge 1$）:

$$N_{i,k}(x)=\frac{x-t_i}{t_{i+k}-t_i}{N}_{i,k-1}(x)+\frac{t_{i+k+1}-x}{t_{i+k+1}-t_{i+1}}{N}_{i+1,k-1}(x)$$

我们将依次构建$N_{i,k}(x)$，从$k=0$ 开始，直到$k=3$。

计算$N_{i,0}(x)$

1. 当$x<1$:
   -$N_{0,0}(x)=0$
   -$N_{1,0}(x)=0$
   -$N_{2,0}(x)=0$
   -$N_{3,0}(x)=0$

2. 当$1\le x<2$:
   -$N_{0,0}(x)=1$
   -$N_{1,0}(x)=0$
   -$N_{2,0}(x)=0$
   -$N_{3,0}(x)=0$

3. 当$2\le x<3$:
   -$N_{0,0}(x)=0$
   -$N_{1,0}(x)=1$
   -$N_{2,0}(x)=0$
   -$N_{3,0}(x)=0$

4. 当$3\le x<4$:
   -$N_{0,0}(x)=0$
   -$N_{1,0}(x)=0$
   -$N_{2,0}(x)=1$
   -$N_{3,0}(x)=0$

5. 当$4\le x<5$:
   -$N_{0,0}(x)=0$
   -$N_{1,0}(x)=0$
   -$N_{2,0}(x)=0$
   -$N_{3,0}(x)=1$

6. 当$x\ge 5$:
   -$N_{0,0}(x)=0$
   -$N_{1,0}(x)=0$
   -$N_{2,0}(x)=0$
   -$N_{3,0}(x)=0$

计算$N_{i,1}(x)$

通过定义计算$N_{i,1}(x)$：

$$N_{0,1}(x)=\frac{x-1}{1-1}{N}_{0,0}(x)+\frac{2-x}{2-1}{N}_{1,0}(x)(\text{当 }1\le x<2\text{ 时为 }0\text{，其余为 }0)$$
$$N_{1,1}(x)=\frac{x-2}{2-1}{N}_{1,0}(x)+\frac{3-x}{3-2}{N}_{2,0}(x)(\text{当 }2\le x<3)$$
$$N_{2,1}(x)=\frac{x-3}{3-2}{N}_{2,0}(x)+\frac{4-x}{4-3}{N}_{3,0}(x)(\text{当 }3\le x<4)$$
$$N_{3,1}(x)=\frac{x-4}{4-3}{N}_{3,0}(x)+\frac{5-x}{5-4}{N}_{4,0}(x)(\text{当 }4\le x<5)$$

注意: 在计算时，随后的$N_{i,k}$ 都基于前一个阶的$N_{i,k-1}$ 基函数构建，所以实际上可以通过这一过程继续计算$N_{i,2}$ 和$N_{i,3}$。

4. 求解控制点系数：

将$S(x)=\sum C_i{N}_{i,3}(x)$带入线性方程组中，来使得拟合通过给定的 ( (x,y) ) 数据点。

对于拟合样条，设定如下方程说明样条函数在数据点的值：

$$\begin{array}{rl}S(1)& =C_0{N}_{0,3}(1)+C_1{N}_{1,3}(1)+C_2{N}_{2,3}(1)+C_3{N}_{3,3}(1)+C_4{N}_{4,3}(1)=2\\ S(2)& =C_0{N}_{0,3}(2)+C_1{N}_{1,3}(2)+C_2{N}_{2,3}(2)+C_3{N}_{3,3}(2)+C_4{N}_{4,3}(2)=3\\ S(3)& =C_0{N}_{0,3}(3)+C_1{N}_{1,3}(3)+C_2{N}_{2,3}(3)+C_3{N}_{3,3}(3)+C_4{N}_{4,3}(3)=5\\ S(4)& =C_0{N}_{0,3}(4)+C_1{N}_{1,3}(4)+C_2{N}_{2,3}(4)+C_3{N}_{3,3}(4)+C_4{N}_{4,3}(4)=4\\ S(5)& =C_0{N}_{0,3}(5)+C_1{N}_{1,3}(5)+C_2{N}_{2,3}(5)+C_3{N}_{3,3}(5)+C_4{N}_{4,3}(5)=6\end{array}$$

将上述方程以 $C_i$形式组织，将能得到线性方程组$AX=B$，其中 $X$ 是控制点系数的向量，$B$ 是观测值 $y$。

示例简化：

实际手工计算B样条系数较为复杂，通常使用计算工具（如Python中的库）进行求解。

5. 构建样条函数：

   利用求得的控制点系数和基函数，组合成最终的样条函数。
   一旦我们得到了控制点 ( C_i )，我们可以将它们代入B样条的表达式中，得到最终的样条函数：
   $$S(x)=C_0{N}_{0,3}(x)+C_1{N}_{1,3}(x)+C_2{N}_{2,3}(x)+C_3{N}_{3,3}(x)+C_4{N}_{4,3}(x)$$

3.4 Python实现

利用 scipy 库中的 make_lsq_spline 和 BSpline 可以实现B样条的拟合。以下为具体实现示例：

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import make_lsq_spline, BSpline

# 数据点
x = np.array([1, 2, 3, 4, 5])
y = np.array([2, 3, 5, 4, 6])

# 设置样条阶数
k = 3  # 三次B样条

# 构建节点向量
# 为避免边界效应，通常在两端添加k个重复节点
t_internal = np.linspace(x.min(), x.max(), len(x) - k + 1)
t = np.concatenate((
    [x.min()] * k,  # 前k个节点
    t_internal,
    [x.max()] * k   # 后k个节点
))

# 构建最小二乘B样条
spl = make_lsq_spline(x, y, t, k)

# 插值用于绘图
x_new = np.linspace(x.min(), x.max(), 400)
y_new = spl(x_new)

# 获取基函数信息
knots = spl.t      # 节点向量
coeffs = spl.c     # 系数
degree = spl.k     # 样条阶数
n_basis = len(coeffs)  # 基函数数量

# 创建一个图形，显示拟合曲线和基函数
fig, ax = plt.subplots(figsize=(12, 8))

# 绘制原始数据点
ax.scatter(x, y, color='red', label='数据点')

# 绘制拟合的B样条曲线
ax.plot(x_new, y_new, label='B-spline 拟合曲线', linewidth=2, color='blue')

# 绘制每个基函数
for i in range(n_basis):
    # 创建一个系数向量，只有第i个系数为1，其余为0
    c = np.zeros(n_basis)
    c[i] = 1
    # 构建相应的基函数B样条
    basis_spline = BSpline(knots, c, degree)
    # 计算基函数在x_new处的值
    y_basis = basis_spline(x_new)
    # 绘制基函数，使用透明度以便更好地显示
    ax.plot(x_new, y_basis, linestyle='--', label=f'基函数 {i + 1}')

# 设置图例，避免重复标签
handles, labels = ax.get_legend_handles_labels()
by_label = dict(zip(labels, handles))
ax.legend(by_label.values(), by_label.keys(), loc='upper left', fontsize='small')

# 设置标题和标签
ax.set_title('B样条拟合及其基函数')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.grid(True)

plt.show()

# 可选：单独绘制每个基函数
fig, axes = plt.subplots(n_basis, 1, figsize=(12, 3 * n_basis), sharex=True)

for i in range(n_basis):
    c = np.zeros(n_basis)
    c[i] = 1
    basis_spline = BSpline(knots, c, degree)
    y_basis = basis_spline(x_new)
    axes[i].plot(x_new, y_basis, color='purple')
    axes[i].set_ylabel(f'B{i + 1}(x)')
    axes[i].set_title(f'基函数 {i + 1}')
    axes[i].grid(True)

axes[-1].set_xlabel('x')
plt.tight_layout()
plt.show()

运行结果：

3.5 优缺点

优点：

• 局部控制性强：调整一个控制点只影响样条曲线的局部区域，便于精细调整。
• 数值稳定性好：递归构建的B基函数具有良好的数值性质，适合大规模数据处理。
• 灵活性高：通过调整节点数量和基函数阶数，可以灵活地适应不同的数据特征。
• 高效性：适合计算机实现，能够高效处理复杂的拟合任务。

缺点：

• 参数选择复杂：需要选择合适的节点向量和阶数，参数选择不当可能影响拟合效果。
• 实现复杂：相对于简单的线性样条，多项式样条和B样条的手工计算和实现更为复杂。
• 对数据分布敏感：节点向量的选择对拟合结果影响较大，非均匀数据分布时需谨慎选择节点。

---

4. 自然样条（Natural Splines）

4.1 原理

自然样条（Natural Splines）是一种具有特定边界条件的三次多项式样条。在自然样条中，最两端的边界处二阶导数为零，即曲线在边界处呈线性行为。这种边界条件减少了样条的振荡，提高了拟合的稳定性，使得样条函数在两端更为“自然”。

4.2 公式

自然三次样条函数与一般的三次多项式样条类似，但在边界处满足以下条件：

$$S^{\prime \prime }(x_0)=0\text{和}{S}^{\prime \prime }(x_n)=0$$

其中，$x_0$ 和$x_n$ 分别为自变量的最小值和最大值。

4.3 案例及手工计算过程

示例数据点：

$$x=[0,1,2,3,4,5],y=[0,0.5,2,1.5,1,1.8]$$

目标：使用自然样条对数据进行拟合。

步骤：

1. 确定分段区间和样条函数：

   • 区间 $[0,1]$：$S_0(x)=a_0+b_{0}x+c_0{x}^2+d_0{x}^3$
   • 区间 $[1,2]$：$S_1(x)=a_1+b_1(x-1)+c_1(x-1{)}^2+d_1(x-1{)}^3$
   • …
   • 区间 $[4,5]$：$S_4(x)=a_4+b_4(x-4)+c_4(x-4{)}^2+d_4(x-4{)}^3$

2. 应用插值条件：

   对每个区间端点，确保样条函数通过数据点：

   $$S_i(x_i)=y_i\text{和}{S}_i(x_{i+1})=y_{i+1}$$

3. 应用连续性条件：

   保证样条函数在内部节点处的一阶和二阶导数连续：

   $$S_i^{\prime }(x_{i+1})=S_{i+1}^{\prime }(x_{i+1})$$
   $$S_i^{\prime \prime }(x_{i+1})=S_{i+1}^{\prime \prime }(x_{i+1})$$

4. 应用自然边界条件：

   在最左端$x_0$ 和最右端$x_n$，样条函数的二阶导数为零：

   $$S_0^{\prime \prime }(x_0)=0,S_{n-1}^{\prime \prime }(x_n)=0$$

5. 构建方程组并求解：

   对于多个区间和样条函数的系数，构建方程组并解出所有系数。这一过程手工计算繁琐，通常借助计算工具完成。

注：手工计算自然样条的全部过程涉及大量方程求解，建议使用计算工具（如Python的 scipy.interpolate.CubicSpline）进行实现。

4.4 Python实现

利用 scipy 库中的 CubicSpline 可以方便地构建自然样条。以下为具体实现示例：

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import CubicSpline

# 数据点
x = np.array([0, 1, 2, 3, 4, 5])
y = np.array([0, 0.5, 2, 1.5, 1, 1.8])

# 构建自然样条
cs = CubicSpline(x, y, bc_type='natural')

# 插值
x_new = np.linspace(0, 5, 500)
y_new = cs(x_new)

# 绘图
plt.figure(figsize=(8, 6))
plt.scatter(x, y, color='red', label='数据点')
plt.plot(x_new, y_new, label='自然样条', linewidth=2)
plt.title('自然样条拟合示例')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)
plt.show()

运行结果：

• 绘图如下：

  

4.5 优缺点

优点：

• 边界自然：在两端的二阶导数为零，使得样条函数在边界处呈线性行为，减少了端点处的振荡。
• 平滑性好：保持了一阶和二阶导数的连续性，提供了较为平滑的拟合曲线。
• 适用性广：适用于各种数据分布，尤其在要求边界条件自然的应用场景中表现良好。

缺点：

• 边界条件限制：仅适用于二阶导数为零的边界条件，限制了在边界处进一步的形状调整。
• 对异常值敏感：与多项式样条类似，受噪声或异常值影响较大，可能导致拟合效果不佳。
• 计算复杂度：对于大量数据点，求解三次样条系数需要较高的计算成本。

---

5. 平滑样条（Smoothing Splines）

5.1 原理

平滑样条（Smoothing Splines）不仅用于插值，还旨在通过平滑拟合来减小数据中的噪声影响。其主要思想是通过引入一个平滑参数（通常记为$\lambda$），在拟合过程中权衡数据拟合的精度与函数的平滑性。平滑样条通过最小化以下目标函数实现平衡：

$$\sum _{i=1}^n(y_i-S(x_i){)}^2+\lambda \int (S^{\prime \prime }(x){)}^{2}dx$$

其中，第一项保证了拟合精度，第二项是对函数二阶导数的惩罚项，用以控制函数的平滑性。参数$\lambda$控制两者之间的权衡关系：$\lambda$ 越大，函数越光滑；$\lambda$ 越小，函数对数据的拟合越精确。

5.2 公式

平滑样条的优化目标是：

$$\underset{S}{\min }\left\{\sum _{i=1}^n(y_i-S(x_i){)}^2+\lambda \int (S^{\prime \prime }(x){)}^{2}dx\right\}$$

其中，$S(x)$ 是样条函数，$\lambda \ge 0$ 是平滑参数。

5.3 案例及手工计算过程

示例：

假设有如下带噪声的数据点，用平滑样条进行拟合。

$$x=[0,1,2,3,4,5,6,7,8,9,10]$$
$$y=[0.1,0.9,2.1,2.9,4.2,5.1,6.0,7.2,8.1,9.0,10.1]$$

目标：使用平滑样条对数据进行拟合，减少噪声影响，捕捉整体趋势。

步骤：

1. 选择平滑参数$\lambda$：

   -$\lambda$ 的选择影响拟合的光滑程度。常通过交叉验证等方法选择最优参数。

2. 构建目标函数：

   $$\underset{S}{\min }\left\{\sum _{i=1}^{11}(y_i-S(x_i){)}^2+\lambda \int (S^{\prime \prime }(x){)}^{2}dx\right\}$$

3. 求解最优化问题：

   • 通过数学方法或数值算法（如基于自然样条的闭式解）求解最优的样条函数$S(x)$。

4. 分析拟合结果：

   • 观察拟合曲线的平滑性与数据拟合精度，调整$\lambda$ 以获得理想的平衡。

注：手工计算平滑样条涉及复杂的积分和最优化求解，通常采用计算工具实现。

5.4 Python实现

利用 scipy 库中的 UnivariateSpline 可以方便地构建平滑样条。以下为具体实现示例：

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import UnivariateSpline

# 生成带噪声的数据
np.random.seed(0)
x = np.linspace(0, 10, 50)
y_true = np.sin(x)
noise = np.random.normal(scale=0.5, size=x.shape)
y_noisy = y_true + noise

# 构建平滑样条
lambda_val = 1  # 平滑参数，可调整
spline = UnivariateSpline(x, y_noisy, s=lambda_val)

# 插值
x_new = np.linspace(0, 10, 1000)
y_new = spline(x_new)
y_smooth = np.sin(x_new)  # 真实曲线

# 绘图
plt.figure(figsize=(10, 6))
plt.scatter(x, y_noisy, color='gray', alpha=0.5, label='带噪声数据')
plt.plot(x_new, y_new, color='red', label='平滑样条', linewidth=2)
plt.plot(x_new, y_smooth, color='blue', linestyle='--', label='真实曲线')
plt.title('平滑样条拟合示例')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)
plt.show()

运行结果：

调整平滑参数$\lambda$：

• $ \lambda$ 较小：样条函数过于贴合数据，可能捕捉到噪声，拟合曲线较为复杂。
• $ \lambda$ 较大：样条函数更加平滑，可能忽略部分数据细节，减少拟合精度。

5.5 优缺点

优点：

• 有效平滑噪声：通过引入平滑参数，能够有效地减小数据中的噪声影响，捕捉数据的整体趋势。
• 参数控制灵活：平滑参数$\lambda$ 提供了对拟合平滑度的灵活控制，适应不同的应用需求。
• 适用范围广：适用于各种需要平滑拟合的应用场景，如信号处理、统计建模等。

缺点：

• 参数选择依赖：平滑参数$\lambda$ 的选择需要经验或交叉验证，增加了模型选择的复杂性。
• 可能过度平滑：若$\lambda$ 选择过大，可能导致拟合曲线过于平滑，忽略数据中的重要细节。
• 计算复杂度较高：与简单的插值样条相比，平滑样条需要额外的优化步骤，计算成本更高。

---

6. 惩罚样条（Penalized Splines, P-splines）

6.1 原理

惩罚样条（P-splines）是一种结合了B样条基函数和差分惩罚的样条方法。P-splines通过在拟合过程中对B样条系数施加差分惩罚，以控制样条的光滑性。它们利用低阶的B样条作为基函数，并通过调整惩罚项来平衡拟合精度与曲线的平滑度。P-splines在处理大规模数据和高维数据时表现良好，具有灵活性和稳定性。

6.2 公式

P-splines的目标函数为：

$$\sum _{i=1}^n(y_i-\sum _{j=1}^k{c}_j{B}_j(x_i){)}^2+\lambda \sum _{j=2}^k({\Delta }^m{c}_j{)}^2$$

其中：

-$B_j(x_i)$ 是B样条基函数。
-$c_j$ 是B样条的控制点系数。
-${\Delta }^m$ 是 $m$ 阶的差分操作符。
-$\lambda$ 是惩罚参数，控制惩罚项的强度。

通过引入差分惩罚项，P-splines能够在保持拟合精度的同时，确保样条的平滑性。

6.3 案例及手工计算过程

示例数据点：

$$\begin{array}{rl}x& =[0,1,2,3,4,5,6,7,8,9,10]\\ y& =[0.1,0.9,2.1,2.9,4.2,5.1,6.0,7.2,8.1,9.0,10.1]\end{array}$$

目标：使用P-splines对数据进行拟合，减小噪声影响，同时保持样条的平滑性。

步骤：

1. 选择B样条基函数：

   • 选择基函数的阶数和节点数量。通常，较多的基函数和较低的阶数可以提供更灵活的拟合。

   例如，选择三次B样条（阶数$k=3$），并设置适当数量的节点。

2. 构建差分惩罚矩阵：

   • 选择差分的阶数$m$。通常选择$m=2$，即二阶差分。

3. 设定惩罚参数$\lambda$：

   -$\lambda$ 控制惩罚项的强度。通过试验或交叉验证选择最优的$\lambda$。

4. 构建线性模型：

   表示为：

   $$y=Bc+ϵ$$

   其中，$B$ 是基函数矩阵，$c$ 是控制点系数，$ϵ$ 是误差项。

5. 最小化目标函数：

   采用岭回归（L2惩罚）求解控制点系数$c$。

6. 构建拟合样条函数：

   使用求得的控制点系数和B基函数，组合成最终的P-spline函数。

注：手工计算P-splines涉及复杂的矩阵运算和优化步骤，通常使用计算工具（如Python的 patsy 和 sklearn 库）进行实现。

6.4 Python实现

利用 patsy 和 sklearn 库可以方便地实现P-splines。以下为具体实现示例：

import numpy as np
import matplotlib.pyplot as plt
from patsy import bs
from sklearn.linear_model import Ridge

# 生成数据
np.random.seed(0)
x = np.linspace(0, 10, 100)
y = np.sin(x) + np.random.normal(scale=0.3, size=x.shape)

# 构建B样条基函数
degree = 3
n_splines = 20
B = bs(x, df=n_splines, degree=degree, include_intercept=False)

# 设置差分惩罚矩阵（这里使用二阶差分）
from scipy.sparse import diags
D = diags([1, -2, 1], [0, 1, 2], shape=(n_splines-2, n_splines)).toarray()

# 构建惩罚项
lambda_val = 10  # 惩罚参数，可调整

# 进行岭回归作为惩罚
model = Ridge(alpha=lambda_val, fit_intercept=False)
model.fit(B, y)
y_new = model.predict(B)

# 绘图
plt.figure(figsize=(10, 6))
plt.scatter(x, y, color='gray', alpha=0.5, label='带噪声数据')
plt.plot(x, y_new, color='blue', label='P-spline', linewidth=2)
plt.title('惩罚样条（P-splines）拟合示例')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)
plt.show()

运行结果：

调节惩罚参数$\lambda$：

• $ \lambda$ 较大：增加惩罚项，样条函数更加平滑，减少对噪声的敏感性。
• $ \lambda$ 较小：减小惩罚项，样条函数更贴合数据，可能捕捉到更多噪声细节。

6.5 优缺点

优点：

• 灵活的平滑控制：通过惩罚参数和差分阶数，能够灵活控制样条的平滑程度和曲线复杂度。
• 局部和全局控制：结合B样条的局部支持性和惩罚项的全局平滑性，P-splines兼具两者的优点。
• 适合大规模数据：数值稳定性好，适合处理大规模和高维数据。
• 自动化拟合：可通过交叉验证等方法自动选择最佳的惩罚参数，简化模型选择过程。

缺点：

• 参数选择复杂：需要同时选择基函数的数量、阶数及惩罚参数，调参过程较为复杂。
• 计算复杂度较高：尤其在基函数数量较多或数据量较大时，计算成本显著增加。
• 模型解释性较低：相比简单的多项式拟合，P-splines的控制点和惩罚参数使得模型解释性降低。

---

总结

样条函数是一类强大的数据拟合与插值工具，不同类型的样条函数各具特色，适用于不同的应用场景。本文系统介绍了六种典型的样条函数，包括线性样条、多项式样条、B样条、自然样条、平滑样条和惩罚样条，详细阐述了它们的原理、公式、具体案例及手工计算过程、Python实现代码以及各自的优缺点。

选择合适的样条函数类型及其参数需要根据具体的数据特性、应用需求及计算资源进行权衡。Python中强大的科学计算库（如 scipy、numpy、patsy、sklearn 等）为样条函数的实现提供了便捷的工具，使得在实际应用中能够高效地进行样条拟合与分析。

---

如果这篇文章对你有一点点的帮助，欢迎点赞、关注、收藏、转发、评论哦！
我也会在微信公众号“智识小站”坚持分享更多内容，以期记录成长、普及技术、造福后来者！

---