Syslab使用手册
⋮
2026a

# spaps

平滑样条曲线

函数库: TyCurveFitting

# 语法

sp, values, rho = spaps(x,y,tol)
_ = spaps(x,y,tol,w,m)
_ = spaps((x1,...,xr),y,tol,___)

# 说明

sp, values, rho = spaps(x, y, tol) 以 B 型返回与给定点 (x[j],y[:,j]),j=1:length(x) 相差在给定容度 tol 以内最平滑的曲线 f，光滑值 values 以及实际使用的 。当仅需一个输出返回时，可以使用 sp = spaps(x,y,tol)[1]。

数据值 y[:,j] 为可以为标量、向量、矩阵或者多维数组。位于同一数据网格中的数据点将被其加权均值替代，其权为对应权的和，容度 tol 也会相应减少。

values 与 fnval(sp,x) 相同。此处曲线 f 与给定数据之间的距离为

其中权 w 默认选为使得 E(f) 为对 的复合梯形规则近似， 表示 z 中元素的平方和。

最光滑意味着下述粗糙测度的最小化

其中 表示 f 的 m 阶导数。m 默认为2，粗糙测度权函数 默认为 1， 这使得 f 为三次光滑样条。

当 tol 为非负值时，样条 f 为使得 最小的唯一值，其中光滑参数 （可选返回）选取使得 E(f) 与 tol 相同。因此，当 m 为 2 时，所得结果转为 pp 型后，应与 csaps(x, y, rho/(rho+1)) 仅有舍入误差的不同。当 tol 为 0 时，则会返回阶为 2m 的“普通”或者说是参变的样条插值。对于足够大的 tol，会返回对数据的度数小于 m 的最小二乘多项式估计。当 tol 为负值， 。

粗糙度量中的权函数 默认为取值为 1 的常值函数，但更一般地也可以选为分段常值函数，其断点尽在数据网格内。假设向量 x 严格增，可以通过输入 tol 为与 x 同样长度的向量来使 为分段常值函数，此时 tol[i] 为 在区间 (x[i-1], x[i]), i=2:length(x) 上取的常值，而 tol[1] 则为容度。

_ = spaps(x, y, tol, w, m) 让您决定权向量 w 以及/或者整数 m，此通过输入参数完成。对于这种用法，w 必须为与 x 相同长度的非负向量，m 必须为 1（表示分段线性光滑样条）、2（表示默认的三次光滑样条）或 3（表示五次光滑样条）。

如果结果的光滑样条 sp 需要在其基本区间以外求值，则需使用 fnxtr(sp, m) 来替代，以保证在基本区间以外，其 m 阶导数为 0。

_ = spaps((x1,...,xr), y, tol, ___)以 B 形式返回有 r 个变元的张量光滑样条，该样条与给定网格值大体在容度内。对于散点值，请使用 tpaps。现在 y 应给定对应的网格值，其中 size(y) 在 f 返回标量时与 (length(x1),...,length(xr)) 相同，在 f 返回 d 向量时与 (d,length(x1),...,length(xr)) 相同。tol 必须为有 r 个元素的元组，其中 tol[i] 表示在第 i 步，即对第 i 个变元构造单变元（返回向量）光滑样条，时使用的容度。可选输入 m 必须为 r 向量（其中每个元素必须为 1、2 或 3），可选输入 w 必须为长为 r 的元组，其中 w[i] 要么空（表示选择默认值）要么为与 xi 相同长度的正向量。示例

# 示例

来自噪声数据的两条三次光滑样条的比较

该代码返回对与非噪声数据足够接近的噪声数据的拟合，因为前者来源于一个变化很慢的函数，也因此使用的 TOL 的维数与噪声的维数相匹配。

using TyCurveFitting
using TyPlot
using TyMath
rng = MT19937ar(5489)
x = LinRange(0, 2 * pi, 21)'
y = sin.(x) .+ (rand(rng, 1, 21) .- 0.5) .* 0.2
sp, = spaps(x, y, (0.05)^2 * (x[end] - x[1]))
fnplt(sp)

该代码使用了与之前相同的数据与容度，但在区间右半使粗糙权仅为 0.1，这相对应地构造了更为粗糙但也更好的拟合。

sp1, = spaps(x, y, [(0.025)^2 * (x[end] - x[1]) ones(1, 10) repeat([0.1], 1, 10)])
fnplt(sp1)

最后，将之前的得到的两条三次光滑样条作比较。

fnplt(sp)
hold("on")
fnplt(sp1, "r")
plot(x, y, "ok")
hold("off")
title("Two cubic smoothing splines")
xlabel("The red one has reduced smoothness requirement in right half.")
对二元函数的噪声数据的光滑逼近

这个例子构造了对光滑二元函数的噪声的光滑逼近。注意 ndgrid 的使用，如果使用 meshgrid2 则会报错。

using TyCurveFitting
using TyBase
using TyMath
x = -2:0.2:2
y = -1:0.25:1
xx, yy = ndgrid(x, y)
rng = MersenneTwister(1234)
z = exp.(-(xx .^ 2 + yy .^ 2)).+ (rand(rng,size(xx))-.5)/30
sp, = spaps((x, y), z, 8 / (60^2))
fnplt(sp)

# 输入参数

x - 数据位点
向量 | 数组元组

需要拟合的数据值 y 的数据位点，对单变元指定为一维矩阵，对于多变元字段为元组。样条 f 使得对于每一个数据位点 x 其节点满足对于所有的j，f(x[j])=y[:,j]。

对多变元、网格数据，可以指定 x 为指定每个变元维度的数据位点的元组：f(x1[i],x2[j],...,xn[k])=y[:,i,j,...,k]。

数据类型： Int | Float

y - 拟合数据
向量 | 矩阵 | 多维数组

在样条构造中需要拟合的数据值，指定为一维矩阵或多维数组。数据值 y[:,j] 可以为标量、矩阵或者 n 维数组。在相同位点 x 中给出的数据值将会被平均。

数据类型： Int | Float

tol - 容度
标量

给定数据点的容度，由标量表示。

数据类型： Float

w - 数据权
向量 | 数组元组

数据点的权 w，指定为与 x 相同长度的、元素非负的向量。

数据类型： Int | Float

m - 微分阶数
标量

微分阶数，指定为标量。该值必须为 1，表示分段线性光滑样条，2，表示默认的三次光滑样条，或 3，表示五次光滑样条。

数据类型： Int

# 输出参数

sp - B 样条
样条结构体

样条，以有以下元素的结构体形式返回。

Form - 样条的类型
"B-"

样条的类型，返回为 "B-"，这意味着样条为 B 型。

数据类型： String

Knots - 样条节点的位置
向量 | 数组元组

样条的节点位置，单变元时以向量，多变元时以向量元组形式返回。其中向量含有严格增的元素，分别表示定义多项式片段的定义区间的左右端点。

数据类型： Int | Float

Coefs - 多项式系数
矩阵 | 数组

每个片段的多项式系数，单变元时以矩阵，多变元时以数组形式返回。

数据类型： Int | Float

Number - 多项式片段的数目
标量 | 向量

描述样条的多项式片段的数目，单变元时以标量，多变元时以由各变元的多项式片段的数目组成的向量形式返回。

数据类型： Int

Order - 多项式阶数
标量 | 向量

描述样条的分段多项式的阶数，单变元时以标量，多变元时以各变元阶数组成的向量形式返回。

数据类型： Int

Dim - 维数
标量

目标函数的维数，以标量形式返回。

数据类型： Int

values - 计算值
向量 | 矩阵 | 多维数组

样条在 x 中的点的值，单变元时以向量，多变元时以矩阵或数组形式返回。

数据类型： Int | Float

rho - 光滑参数
标量 |数组元组

计算样条时使用的光滑参数，单变元时以标量，多变元时以标量值元组形式返回。

数据类型： Float

# 算法

该函数使用 Reinsch 的方法[1]，包括他对最优光滑参数生成方程的选择，该选择会保证存在可用的好的初值且 Newton 法能够而且很快地收敛。

# 参考文献

[1] C. Reinsch. "Smoothing by spline functions." Numerical Math. 10, 1967(177-183).

# 另请参阅

csaps | spap2 | spapi | tpaps