3 次平滑化スプライン

CSAPS コマンド

コマンド csaps は、"平滑化" スプラインを提供します。これは、ノイズを含むデータの基となるトレンドの推定をある程度近似する 3 次スプラインです。選択した平滑化パラメーターによって、与えられたデータを平滑化スプラインがどの程度近似するかが決まります。基本情報として、ドキュメンテーションの要約版を次に示します。

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

例: 3 次多項式のノイズを含むデータ

ここでは、試しにいくつかを実行してみます。単純な 3 次式 q(x) := x^3 のデータから開始し、値にノイズを追加して、平滑化パラメーターが .5 になるように値を選択します。次に、基にした 3 次式およびノイズ入りデータと共に、結果の平滑化された値をプロットします。

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' )

この場合は平滑化されすぎています。平滑化パラメーター p を 1 に近づけることで、与えられたデータに近い平滑化スプラインが得られます。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』の 235ff ページ) は、独立変数のスケーリングに大きな影響を受けます。使用される式を簡単に解析すると、epsilon := h^3/16 および隣接するサイト間の平均差異を h として、p の影響を受ける範囲はおよそ 1/(1+epsilon) となります。特に、p = 1/(1+epsilon/100) の場合にデータの忠実な近似、p = 1/(1+epsilon*100) の場合に良好な平滑化が予想されます。

次のプロットに、このマジック ナンバー 1/(1+epsilon) 付近の p の値の平滑化スプラインを示します。この場合、マジック ナンバー 1/(1+epsilon) が非常に 1 に近いため、1-p を確認するとより多くの情報が得られます。

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 により得られる 3 次平滑化スプラインは、その選択方法のみが 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' )

この Figure のタイトルは、、これらのデータの正確な平滑化スプラインを取得するために csaps で使用する p の値を示しています。

さらに、平滑化パラメーターが指定されなかった場合に 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 コマンドには、個々の平滑化スプラインを指定する方法に平滑化パラメーターを使用するか、許容誤差を使用するかの違いがあります。もう 1 つの違いは、spaps は 3 次平滑化スプライン以外にも、線形平滑化スプラインまたは 5 次平滑化スプラインを提供できることです。

2 階微分をできるだけ動かしたくない状況では、5 次平滑化スプラインが 3 次平滑化スプラインよりも優れています。