三次平滑样条
此示例说明如何使用 Curve Fitting Toolbox™ 中的 csaps
和 spaps
命令来构造三次平滑样条。
CSAPS 命令
使用 csaps
命令可提供平滑样条。这是一个三次样条,它在一定程度上遵循了含噪数据中假定的基本趋势。要由您选择的平滑参数确定平滑样条与给定数据的接近程度。以下是基本信息和简要的文档描述:
CSAPS Cubic smoothing spline.
VALUES = CSAPS(X, Y, P, XX)
Returns the values at XX of the cubic smoothing spline for the
given data (X,Y) and depending on the smoothing parameter P, chosen
from the interval [0 .. 1]. This smoothing spline f minimizes
P * sum_i W(i)(Y(i) - f(X(i)))^2 + (1-P) * integral (D^2 f)^2
示例:三次多项式中的含噪数据
下面是一些试运行。我们从一个简单的三次多项式 q(x) := x^3
中的数据开始,用一些噪声污染这些值,并选择 0.5 作为平滑参数值。然后绘制生成的经过平滑处理的值,以及基础三次多项式数据和受污染数据。
xi = (0:.05:1); q = @(x) x.^3; yi = q(xi); randomStream = RandStream.create( 'mcg16807', 'Seed', 23 ); ybad = yi+.3*(rand(randomStream, size(xi))-.5); p = .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 = [.6 .7 .8 .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'}, 'Location', 'NorthWest' )
我们看到平滑样条对平滑参数的选择非常敏感。即使 p
= 0.9,平滑样条仍远离基本趋势,而对于 p
= 1,我们得到了(含噪)数据的插值。
事实上,csapi
(文献 A Practical Guide to Splines 的第 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))/(numel(xi)-1))^3/16; 1 - 1/(1+epsilon)
ans = 7.8124e-06
plot(xi,yi,':',xi,ybad,'x') hold on labels = cell(1,5); for j=1:5 p = 1/(1+epsilon*10^(j-3)); yy(j,:) = csaps(xi,ybad,p,xxi); labels{j} = ['1-p= ',num2str(1-p)]; end plot(xxi,yy) title('Smoothing Splines for Smoothing Parameter Near Its ''Magic'' Value') legend( [{'Exact', 'Noisy'}, labels], 'Location', 'NorthWest' ) hold off
在此示例中,平滑样条对幻数附近的平滑参数的变化非常敏感。离 1 最远的平滑参数似乎是最好的选项,但您可能更喜欢比它更远的那个平滑参数。
p = 1/(1+epsilon*10^3); yy = csaps(xi,ybad,p,xxi); hold on plot( xxi, yy, 'y', 'LineWidth', 2 ) title( sprintf( 'The Smoothing Spline For 1-p = %s is Added, in Yellow', num2str(1-p) ) ) hold off
您还可以为 csaps
提供误差权重,以便对某些数据点给予更多的关注。此外,如果不提供计算位点 xx
,则 csaps
会返回 pp 型平滑样条。
最后,csaps
也可以处理向量值数据,甚至多元网格数据。
SPAPS 命令
spaps
命令提供的三次平滑样条与在 csaps
中构造的样条仅在选择方式上有所不同。以下是 spaps
的简要文档描述:
SPAPS Smoothing spline.
[SP,VALUES] = SPAPS(X,Y,TOL) returns the B-form and, if asked,
the values at X, of a cubic smoothing spline f for the given
data (X(i),Y(:,i)), i=1,2, ..., n.
The smoothing spline f minimizes the roughness measure
F(D^2 f) := integral ( D^2 f(t) )^2 dt on X(1) < t < X(n)
over all functions f for which the error measure
E(f) := sum_j { W(j)*( Y(:,j) - f(X(j)) )^2 : j=1,...,n }
is no bigger than the given TOL. Here, D^M f denotes the M-th
derivative of f. The weights W are chosen so that E(f) is the
Composite Trapezoid Rule approximation for F(y-f).
f is constructed as the unique minimizer of
rho*E(f) + F(D^2 f),
with the smoothing parameter RHO so chosen so that E(f) equals
TOL. Hence, FN2FM(SP,'pp') should be (up to roundoff) the same
as the output from CPAPS(X,Y,RHO/(1+RHO)).
容差与平滑参数
与 csaps
要求的平滑参数 p
相比,为 spaps
提供合适的容差可能更容易一些。在之前的示例中,我们添加了来自区间 0.3*[-0.5 .. 0.5]
的均匀分布的随机噪声。因此,我们可以将该噪声下的误差测度 e
的值作为噪声 tol
的合理值。
tol = sum((.3*(rand(randomStream, size(yi))-.5)).^2);
此图显示由 spaps
构造而得到的平滑样条。请注意,误差权重指定为统一的,即它们在 csaps
中的默认值。
[sp,ys,rho] = spaps(xi,ybad,tol,ones(size(xi))); plot(xi,yi,':',xi,ybad,'x',xi,ys,'r-') title( sprintf( 'Clean Data, Noisy Data, Smoothed Values (1-p = %s )', num2str(1/(1+rho)) ) ); legend( {'Exact','Noisy','Smoothed'}, 'location', 'NorthWest' )
图窗标题显示 p
的值,您将在 csaps
中使用该值来获得这些数据的此平滑样条。
此外,此处是在没有给定平滑参数时由 csaps
提供的平滑样条。在这种情况下,csaps
通过特定的专门过程选择参数,该过程尝试定位平滑样条对平滑参数最敏感的区域(类似于前面所讨论的)。
hold on plot(xxi,fnval(csaps(xi,ybad),xxi),'-') title('Clean Data, Noisy Data, Smoothed Values') legend({'Exact' 'Noisy' 'spaps, specified tolerance' ... 'csaps, default smoothing parameter'}, 'Location', 'NorthWest' ) hold off
CSAPS 与 SPAPS
csaps
和 spaps
命令的不同之处在于:一个是通过平滑参数指定特定平滑样条,一个是通过容差指定。另一个区别是 spaps
除了提供三次平滑样条以外,还可以提供线性或五次平滑样条。
在希望二阶导数尽可能少移动的情况下,五次平滑样条优于三次平滑样条。