# csaps
三次平滑样条
函数库: TyCurveFitting
# 语法
pp,P = csaps(x,y)
pp,P = csaps(x,y,p)
pp,P = csaps(x,y,p,[],w)
values,P = csaps(x,y,p,xx)
values,P = csaps(x,y,p,xx,w)
_ = csaps((x1,...,xm),y,___)
# 说明
提示
使用 fit 函数可以使用更简单但更不灵活的方式来生成光滑样条。
pp,P = csaps(x,y) 将三次平滑样条插值返回到 pp 型中的给定数据 (x,y) 以及在最终样条结果中使用的平滑参数的值，无论是否指定 p。返回的 P 值对于实验很有用，可以从 pp,P=csaps(x,y) 开始并获得对 p 的合理的第一次猜测。数据点 x[j] 处的样条 f 的值近似于 j = 1:length(x) 时的数据值 y[:,j]。
平滑样条 f 最小化下式：
这里 n 是 x 的元素数，积分区间为包含 x 的所有元素的最小区间。
误差测度权重
为了在基本区间外求值平滑样条，必须首先将其外推。使用指令pp=fnxtr(pp)来确保数据位点扩张的区间外的二阶导数为零。
pp, P = csaps(x,y,p) 指定平滑参数 p 。还可以提供向量形式的 p 作为粗糙测度中的权重 λ，该向量的第一个元素是 p，第 i 个元素是区间 (x[i-1],x[i]) 上的 λ 值。示例
pp, P = csaps(x,y,p,[],w)同样指定了误差测度中的权重w。示例
values, P = csaps(x,y,p,xx) 使用平滑参数 p 并返回在点 xx 处计算的平滑样条的值以及最终样条结果中使用的平滑参数的值 P。values 的值与 fnval(csaps(x,y,p),xx) 相同。 示例
values, P = csaps(x,y,p,xx,w) 使用平滑参数 p 和误差测度权重 w，并返回在点 xx 处计算的平滑样条的值。此语法与 fnval(csaps(x,y,p,[],w),xx) 相同。
_ = csaps((x1,...,xm),y,___) 为 (x1,...,xm) 描述的符合矩形网格上的数据提供 m 变元张量积平滑的pp型样条。可以将此语法与先前语法中的任何参数一起使用。
# 示例
使用不同的平滑参数拟合样条
使用具有不同平滑参数值 p 的 csaps 函数拟合平滑样条。使用介于 0 和 1 之间的值 p 来查看它们如何影响拟合样条的形状和紧密度。
加载 titanium 数据集。
using TyCurveFitting
using TyPlot
x, y = titanium()
当 p = 0 ，s0 是拟合数据的最小二乘直线。当 p = 1, s1 是变分或自然三次样条插值。
对于 0 < p < 1 , sp 是一个平滑样条，它是两个极端之间的权衡：比插值 s1 更平滑，比直线 s0 更接近数据。
p = 0.00009
s0, = csaps(x,y,0)
sp, = csaps(x,y,p)
s1, = csaps(x,y,1)
figure()
fnplt(s0)
hold("on")
fnplt(sp)
fnplt(s1)
plot(x,y,"ko")
hold("off")
title("Smoothing splines with different values for p")
legend(["p = 0","p = $p", "p = 1"])
调整平滑参数和权重
调整平滑参数、误差测度权重和粗糙测度权重。
创建带有噪声的正弦曲线。
using TyCurveFitting
using TyPlot
using TyMath
rng = MT19937ar(5489)
x = LinRange(0,2*pi,21)
y = sin.(x) + (vec(rand(rng,1,21)).-0.5).*0.3
创建数据的拟合平滑样条。指定平滑参数 p = 0.4 和与数据相关的误差测度权重。
pp, = csaps(x,y,0.4,[],[ones(1,10) repeat([5],1,10) 0])
该函数返回对噪声数据的平滑拟合，该数据更接近右半部分的数据，因为那里的误差测度权重要大得多。 请注意，最后一个数据点的误差测度为零，这从拟合中排除了该点。
现在使用相同的数据、平滑参数和误差测度权重拟合平滑样条，但调整粗糙度测度权重。
p = [.4 ones(1,10) repeat([.2],1,10)]
xx = []
w = [ones(1,10) repeat([5],1,10) 0]
在区间的右半部分，粗糙度测度权重仅为 0.2。相应地，拟合更粗糙，但在数据的右侧更接近（除了最后一个数据点，它被忽略）。
绘制两个拟合进行比较。
pp1, =csaps(x,y, [.4 ones(1,10) repeat([.2],1,10)], [], [ones(1,10) repeat([5],1,10) 0])
figure()
hold("on")
fnplt(pp,"b")
fnplt(pp1,"r--")
plot(x,y,"ok")
hold("off")
ylim([-1.5 1.5])
title("Cubic smoothing spline, with right half treated")
legend(["Larger error weight", "Larger error and smaller roughness weight"])
平滑二元数据
将平滑样条拟合到由增加了均匀噪声的峰值函数生成的双变元数据。使用 csaps 获取新的平滑数据点，并使用 csaps 为拟合确定的平滑参数。
创建网格。在本例中，网格是一个 51×61 的均匀网格。
using TyCurveFitting
using TyPlot
using TyBase
using TyMath
rng = MT19937ar(5489)
x = (LinRange(-2,3,61),LinRange(-3,3,61))
xx,yy = ndgrid(x[1],x[2])
使用 peaks 函数和区间
y = peaks(xx, yy)
noisy = y .+ (rand(rng, size(y,1),size(y,2)) .- 0.5)
figure()
surf(xx,yy,noisy)
axis("off")
拟合数据。使用 csaps 获取在网格 x 上的平滑数据值，拟合中使用默认平滑参数。拟合图显示仍然存在一些粗糙度。请注意，必须转置数组 sval。
figure()
X1,X2 = meshgrid2(x[1],x[2])
sval,p = csaps(x,noisy,[],x)
surf(X1,X2,Matrix(sval'))
axis("off")
要获得更平滑的近似值，请为 p 指定一个略小于 csaps 默认值的值。
ssval = csaps(x,noisy,0.996,x)[1]
figure()
X3,X4 = meshgrid2(x[1],x[2])
surf(X3,X4,Matrix(ssval'))
axis("off")
# 输入参数
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
p - 平滑参数[0,1] 区间内的标量 | 向量 | 空数组
平滑参数，指定为 0 到 1 之间的标量值。还可以通过提供 p 作为向量来指定粗糙度测度权重 λ 的值。要为多元数据提供粗糙度测度权重，请使用向量。如果提供一个空数组，该函数会根据数据位点 x 选择 p 的默认值，以及粗糙度测度权重 λ 的默认值 1。
平滑参数决定了对 f 平滑或 f 接近数据的矛盾要求的相对权重。对于 p = 0，f 是数据的最小二乘直线拟合。对于 p = 1，f 是变分或自然三次样条插值。随着 p 从 0 移动到 1，平滑样条从一个极端变为另一个极端。
p 的有利范围通常接近 1/(1 +
如果输入 p 为负或为空，则函数使用 p 的默认值。
通过提供 p 作为向量，可以在平滑参数旁边指定粗糙度测度权重 λ。该向量必须与 x 的大小相同，第 i 个元素是区间 (x(i-1)...x(i)) 上的 λ 值， i = 2:length(x)。输入向量 p 的第一个元素是平滑参数 p 的值。通过提供粗糙度测度权重，可以使生成的平滑样条在区间的不同部分更平滑（具有较大的权重值）或更接近于数据（具有较小的权重值）。粗糙度测度权重必须是非负的。
数据类型： Int | Float
w - 误差测度权重向量 | 多维数组
误差测度中的权重 w，指定为与 x 大小相同的非负元素向量。
误差测度中的权重向量 w 的默认值是 ones(size(x))。
数据类型： Float
xx - 计算点一维数组 | 元组
计算样条的计算点，单变元时为一维数组，多变元时为向量元组。样条计算是使用fnval。
数据类型： Int | Float
# 输出参数
pp - 样条结构体样条结构体
pp型的样条，为带有以下元素的结构体。关于 pp 样条的更多信息，参照 pp 型。
Form - 样条格式"pp"
样条的格式，以 "pp" 返回。"pp" 样条是以分段多项式的形式给出。
数据类型： String
Breaks-样条节点位置向量 | 数组元组
样条节点位置，单变元时以向量，多变元时以向量元组形式返回。向量包含严格增元素，逐点分别表示各段多项式定义区间的左右端点。
数据类型： Int | Float
Coefs - 多项式的系数矩阵 | 数组
各段多项式的系数，单变元时以矩阵，多变元时以数组形式返回。
数据类型： Int | Float
Pieces - 多项式分段段数标量 | 向量
描述样条的多项式分段的段数，单变元时以标量，多变元时以各变元段数组成的向量形式返回。
数据类型： Int
Order - 多项式阶数标量 | 向量
描述样条的分段多项式的阶数，单变元时以标量，多变元时以各变元阶数组成的向量形式返回。
数据类型： Int
Dim - 维数标量
目标函数的维数，以标量形式返回。
数据类型： Int
values - 计算样条向量 | 矩阵 | 多维数组
计算样条，单变元时以向量，多变元时以矩阵或多维数组形式返回。在给定的计算点 xx 处计算样条。
数据类型： Int | Float
P - 平滑参数标量 | 向量
用于计算样条曲线的平滑参数，单变元时以标量，多变元以向量形式返回。 P 介于 0 和 1 之间。
数据类型： Float
# 算法
csaps 为 PGS 的 Fortran 进程 SMOOTH 的一个实现。
平滑样条的计算需要求解一个线性系统，其系数矩阵的形式为 p*A + (1-p)*B，矩阵 A 和 B 取决于数据位点 x。 p 的默认值使 p*trace(A) = (1-p)*trace(B)。