Main Content

三次平滑样条

此示例说明如何使用 Curve Fitting Toolbox™ 中的 csapsspaps 命令来构造三次平滑样条。

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/16h 为相邻位点之间的平均差。具体来说,当 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

csapsspaps 命令的不同之处在于:一个是通过平滑参数指定特定平滑样条,一个是通过容差指定。另一个区别是 spaps 除了提供三次平滑样条以外,还可以提供线性或五次平滑样条。

在希望二阶导数尽可能少移动的情况下,五次平滑样条优于三次平滑样条。