# 三次平滑样条
此示例说明如何使用曲线拟合工具箱中 csaps 和 spaps 命令构建三次平滑样条曲线。
# CSAPS 指令
指令 csaps 提供平滑样条。 这是一个或多或少遵循噪声数据中假定的潜在趋势的三次样条曲线。 由您选择的平滑参数决定了平滑样条曲线与给定数据的紧密程度。以下是基本信息，这是帮助文档的简略版本：
CSPS 三次平滑样条。
values = CSPS(X, Y, P, XX)
返回给定数据 (X,Y)的三次平滑样条曲线在 XX 处的值并取决于介于 [0, 1] 之间的平滑参数 P。 该平滑样条 f 最小化
# 示例：来自三次多项式的噪声数据
这里有一些测试运行。 我们从简单三次的数据开始，q(x) := x^3，用一些噪声污染这些值，并选择平滑参数的值为 0.5。 然后绘制得到的平滑值以及基础的三次和受污染的数据。
using TyPlot
using TyCurveFitting
using TyMath
xi = 0:0.05:1
q(x) = x^3
yi = q.(xi)
randomStream = Random.seed!(23)
ybad = yi .+ 0.3 * (rand(randomStream, size(xi, 1), 1) .- 0.5)
p = 0.5
xxi = (0:100) ./ 100
ys, = csaps(xi, ybad, p, xxi)
plot(xi, yi, ":", xi, ybad, "x", xxi, ys, "r-")
title("Clean Data, Noisy Data, Smoothed Values")
legend(["Exact", "Noisy", "Smoothed", "Location", "NorthWest"])
这里的平滑过度了 通过选择更接近 1 的平滑参数 p，我们获得更接近给定数据的平滑样条。 我们尝试 p = .6, .7, .8, .9, 1，并绘制生成的平滑样条曲线。
yy = zeros(5, length(xxi))
p = [0.6 0.7 0.8 0.9 1]
for j = 1:5
yy[j, :], = csaps(xi, ybad, p[j], xxi)
end
hold("on")
plot(xxi,yy')
hold("off")
title("Smoothing Splines for Various Values of the Smoothing Parameter")
legend(["Exact", "Noisy", "p = 0.5", "p = 0.6", "p = 0.7", "p = 0.8",
"p = 0.9", "p = 1.0"])
我们发现平滑样条对于平滑参数的选择十分敏感。即使 p = 0.9，平滑样条也与隐藏趋势差得比较远，而 p = 1 时，我们就得到了（含噪）数据的插值。
事实上，csapi 使用的公式（样条实用指南的第 235 页）对自变量的缩放非常敏感。对所用方程的简单分析表明，p 的敏感范围约为 1/(1+epsilon)，其中 epsilon := h^3/16，h 是相邻位点之间的平均差异。 具体来说，当 p = 1/(1+epsilon/100) 时，拟合与数据很接近，而 p = 1/(1+epsilon*100) 时会得到满意的平滑效果。
下图显示了接近这个魔数 1/(1+epsilon) 的 p 值的平滑样条曲线。 对于这种情况，查看 1-p 会提供更多信息，因为魔数 1/(1+epsilon) 非常接近 1。
epsilon = ((xi[end] - xi[1]) / (length(xi) - 1))^3 / 16
1 - 1 / (1 + epsilon)
ans = 7.81243896530448e-6
plot(xi, yi, ":", xi, ybad, "x")
hold("on")
labels = ones(String, 5)
for j = 1:5
p = 1 / (1 + epsilon * (10.0^(j - 3)))
yy[j, :] = csaps(xi, ybad, p, xxi)[1]
pp = 1 - p
labels[j] = "1-p= $pp"
end
plot(xxi, yy')
title("Smoothing Splines for Smoothing Parameter Near Its 'Magic' Value")
legend([["Exact", "Noisy"]; labels],loc="northwest")
hold("off")
在这个示例中，平滑样条曲线对魔数附近的平滑参数变化非常敏感。 离 1 最远的那个似乎是这些中的最佳选择，但您可能更喜欢除此之外的下面这个。
p = 1 / (1 + epsilon * 10^3)
yy, = csaps(xi, ybad, p, xxi)
hold("on")
plot(xxi, yy, "y", linewidth=2)
pp = 1 - p
title("The Smoothing Spline For 1-p = $pp is Added, in Yellow")
hold("off")
您还可以为 csap 提供错误权重，比其他数据点更关注某些数据点。 此外，如果您不提供计算位点 xx，则 csaps 返回平滑样条的 pp 型。
最后，csaps 还可以处理向量值数据，甚至是多变量的网格数据。
# SPAPS 指令
由命令 spaps 提供的三次平滑样条曲线与在 csaps 中构造的样条曲线的不同之处仅在于它的选择方式。以下是 spaps 文档的缩写版本：
SPAPS 平滑样条曲线。
SP,VALUES = SPAPS(X,Y,TOL) 返回 B 型 sp 以及 给定数据 (X[i],Y[:,i]), i=1,2,...,n 的三次平滑样条 f 在 X 处的值 values。
平滑样条 f 最小化粗糙度测度
其中 f 满足误差测度
不大于给定的 TOL。其中
f 具有形式
其中平滑参数 RHO 满足 E(f) 等于 TOL。因此， fn2fm(SP,"pp") 应该（直到四舍五入）与 cpaps(X,Y,RHO/(1+RHO)) 相同。
# 容差与平滑参数
与 csaps 所需的平滑参数 p 相比，为 spaps 提供合适的容差可能更容易。在我们之前的示例中，我们从区间 0.3 .*[-0.5 .. 0.5] 中添加了均匀分布的随机噪声。 因此，我们可以将 tol 的合理值估计为在这种噪声下的误差测度的值。
tol = sum((0.3 * (rand(randomStream, size(yi)...) .- 0.5)) .^ 2)
该图显示了由 spaps 构建的平滑样条曲线。请注意，误差权重被指定为统一的，这也是它们在 csaps 中的默认值。
sp, ys, rho = spaps(xi, ybad, tol, ones(size(xi)))
plot(xi, yi, ":", xi, ybad, "x", xi, ys, "r-")
title("Clean Data, Noisy Data, Smoothed Values (1-p = $(1/(1+rho)) )")
legend(["Exact", "Noisy", "Smoothed"])
图标题显示了您将在 csaps 中使用的 p 值，以准确获得这些数据的平滑样条曲线。
此外，这是 csaps 在没有给定平滑参数时提供的平滑样条曲线。在这种情况下，csaps 通过某个特定过程选择参数，该过程试图定位平滑样条对平滑参数最敏感的区域（类似于之前的讨论）。
hold("on")
plot(xxi, fnval(csaps(xi, ybad)[1], xxi), "-")
title("Clean Data, Noisy Data, Smoothed Values")
legend(["Exact", "Noisy", "spaps, specified tolerance", "csaps, default smoothing parameter"])
hold("off")
# CSAPS vs. SPAPS
csaps 和 spaps 命令的不同之处在于您通过平滑参数还是容差指定特定平滑样条曲线的方式。 另一个区别是除了三次平滑样条之外，spaps 还可以提供线性或五次平滑样条。
在您希望二阶导数尽可能少移动的情况下，五次平滑样条曲线优于三次平滑样条曲线。