csaps

3 次平滑化スプライン

構文

pp = csaps(x,y)

pp = csaps(x,y,p)

pp = csaps(x,y,p,[],w)

values = csaps(x,y,p,xx)

values = csaps(x,y,p,xx,w)

[___] = csaps({x1,...,xm},y,___)

[___,P] = csaps(___)

説明

メモ

簡単で柔軟性が低い方法で平滑化スプラインを生成するには、曲線フィッターアプリまたは関数 fit を使用してみてください。

pp = csaps(x,y) は、与えられたデータ (x,y) に対して 3 次平滑化スプライン内挿の pp 型を返します。データ サイト x(j) におけるスプライン f の値は、j = 1:length(x) について、データ値 y(:,j) を近似します。

平滑化スプライン f は次を最小化します。

$$p\underset{\text{error measure}}{\underbrace{{\displaystyle \sum _{j=1}^n{w}_j}|y_j-f(x_j){|}^2}}+(1-p)\underset{\text{roughness measure}}{\underbrace{{\displaystyle \underset{}{\overset{}{\int }}\lambda (t)}|D^{2}f(t){|}^{2}dt}}$$

ここで、n は x のエントリ数であり、積分は x のすべてのエントリを含む最小の区間に対応しています。y_j および x_j はそれぞれ、y および x の j 番目のエントリを表します。D²f は、関数 f の 2 階微分を示します。

誤差測定の重み w_j の既定値 は 1 です。粗さ測定における区分定数の重み関数 λ の既定値は定数関数 1 です。既定では、csaps は、指定されたデータ サイト x に基づいて、平滑化パラメーター p の値を選択します。

基本区間外で平滑化スプラインを評価するには、まずその外挿を行わなければなりません。データ サイトによって決まる区間の外部では 2 階微分が必ず 0 となるよう、コマンド pp = fnxtr(pp) を使用します。

pp = csaps(x,y,p) は、平滑化パラメーター p を指定します。また、p を、最初のエントリが p で、i 番目のエントリが区間 (x(i-1),x(i)) における λ の値であるベクトルとして指定することで、粗さ測定の重み λ を指定できます。

pp = csaps(x,y,p,[],w) も、誤差測定の重み w を指定します。

values = csaps(x,y,p,xx) は、平滑化パラメーター p を使用し、点 xx において評価された平滑化スプラインの値を返します。この構文は、fnval(csaps(x,y,p),xx) と同じです。

values = 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(___) は、p を指定したかどうかに関係なく、最終的なスプラインの結果で使用された平滑化パラメーターも返します。この構文は、[pp,P] = csaps(x,y) で開始する可能性のある実験で、p から妥当な初期推定値を得るために役立ちます。

例

異なる平滑化パラメーターを使用したスプラインの近似

関数 csaps で、平滑化パラメーター p に異なる値を使用して、平滑化スプラインを近似します。0 と 1 の各極値の間の p 値を使用して、近似されたスプラインの形状および近似性への影響を確認します。

チタン データセットを読み込みます。

[x, y] = titanium();

p = 0 の場合、s0 はデータの最小二乗直線近似となります。p = 1 の場合、s1 は変分または自然な 3 次スプライン内挿となります。

0 < p < 1 の場合、sp は 2 つの極値間のトレードオフである平滑化スプラインとなります。すなわち、内挿 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 = ' num2str( p )], 'p = 1', 'Location', 'northwest')

平滑化パラメーターおよび重みの調整

平滑化パラメーター、誤差測定の重み、および粗さ測定の重みを調整します。

ノイズを伴う正弦曲線を作成します。

x = linspace(0,2*pi,21); y = sin(x)+(rand(1,21)-.5)*.3;

データに平滑化スプラインを近似します。平滑化パラメーター p = 0.4、およびデータ間で異なる誤差測定の重み w を指定します。

pp = csaps(x,y,0.4,[],[ones(1,10),repmat(5,1,10), 0]);

関数は、右半分の方が誤差測定の重みがはるかに大きいため、そこのデータにはるかに近い、ノイズの多いデータへの滑らかな近似を返します。最後のデータ点は誤差の重み付けが 0 であるため、この点は近似から除外されることに注意してください。

ここで、同じデータ、平滑化パラメーター、誤差測定の重み、および調整された粗さ測定の重みを使用して平滑化スプラインを当てはめます。

pp1 = csaps(x,y, [.4,ones(1,10),repmat(.2,1,10)], [], ...
                    [ones(1,10), repmat(5,1,10), 0]);

区間の右半分における粗さ測定の重みはわずか 0.2 です。そのため、近似は粗くなりますが、データの右側では近くなります (無視される最後のデータ点を除きます)。

比較のため、両方の近似をプロットします。

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 ',...
          'differently'])
legend('Larger error weight', 'Larger error and smaller roughness weight')

二変量データの平滑化

一様なノイズを追加して関数 peaks で生成した 二変量データに平滑化スプラインを近似させます。csaps を使用して、平滑化された新しいデータ点、および csaps によって特定される近似の平滑化パラメーターを取得します。

グリッドを作成します。この例では、グリッドは 51 行 61 列の等間隔グリッドです。

x = {linspace(-2,3,51),linspace(-3,3,61)};
[xx,yy] = ndgrid(x{1},x{2}); 

関数 peaks、および区間 $$[-\frac{1}{2},\frac{1}{2}]$$ の乱数を使用してノイズの多いデータを生成します。

y = peaks(xx, yy);
noisy = y + (rand(size(y)) - 0.5);
figure
surf(xx,yy,noisy)
axis off

データを近似します。csaps を使用して、グリッド x に対して評価される平滑化されたデータ値、および近似に使用される既定の平滑化パラメーターを取得します。

[sval,p] = csaps(x,noisy,[],x);

近似のプロットには、多少粗さが残っていることが示されています。配列 sval を転置しなければならない点に注意してください。

figure
surf(x{1},x{2},sval.')
axis off

やや滑らかな近似では、csaps の既定値より若干小さい p の値を指定します。

ssval = csaps(x,noisy,.996,x);
figure
surf(x{1},x{2},ssval.')
axis off

入力引数

x — データ サイト

ベクトル | cell 配列

近似対象のデータ値 y のデータ サイト。ベクトル、または多変量データの場合は cell 配列として指定します。スプライン f は、j のすべての値について f(x(j)) = y(:,j) を満たすよう、各データ サイト x の節点によって作成されます。

多変量グリッド データの場合、x を、各変数の次元のデータ サイトを指定する cell 配列として指定できます。f(x1(i),x2(j),...xn(k)) = y(:,i,j,...,k)

データ型: single | double

y — 近似対象のデータ値

ベクトル | 行列 | 配列

スプライン作成時の近似対象のデータ値。ベクトル、行列、または配列として指定します。データ値 y(:,j) は、スカラー、行列、または n 次元配列にすることができます。同じデータ サイト x のデータ値は平均化されます。

データ型: single | double

p — 平滑化パラメーター

範囲 [0,1] のスカラー | ベクトル | cell 配列 | 空の配列

平滑化パラメーター。0 ～ 1 のスカラー値、または多変量データ値の cell 配列として指定します。p をベクトルとして指定することで、粗さ測定の重み λ の値を指定することもできます。多変量データの粗さ測定の重みを指定するには、ベクトルの cell 配列を使用します。空の配列を指定した場合、関数は、p についてはデータ サイト x に基づく既定値、粗さ測定の重み λ については既定値 1 を選択します。

平滑化パラメーターは、f を滑らかにする、または f をデータに近づけるという相反する要求に適用する、相対的な重みを決定します。p = 0 の場合、f はデータの最小二乗直線近似となります。p = 1 の場合、f は変分または自然な 3 次スプライン内挿となります。p が 0 から 1 に移動すると、平滑化スプラインは一端からもう一方の端に変化します。

p の適切な範囲は多くの場合 1/(1 + h³/6) の近傍です。ここで、h はデータ サイトの平均間隔です。関数は、p の既定値をこの範囲内で選択します。等間隔のデータの場合、p = 1(1 + h³/60) では密接な近似を、p = 1/(1 + h³/0.6) ではある程度十分な平滑化を想定できます。p > 1 は入力できますが、これを選択した場合、平滑化スプラインが変分 3 次スプライン内挿よりもさらに粗くなります。

入力 p が負または空の場合、p の既定値が関数で使用されます。

p をベクトルとして指定することで、平滑化パラメーターと共に、粗さ測定の重み λ を指定できます。このベクトルは、x と同じサイズでなければなりません。ここで、i = 2:length(x) に対して、区間 (x(i-1)...x(i)) で i 番目のエントリを λ の値とします。入力ベクトル p の最初のエントリは平滑化パラメーター p の望ましい値となります。粗さ測定の重みを指定することで、区間のさまざまな場所において、結果の平滑化スプラインをより滑らかにしたり (重み値を大きくした場合)、データにより近づけたり (重み値を小さくした場合) することができます。粗さ測定の重みは非負でなければなりません。

p の選択が困難であっても、y のノイズのサイズをある程度は把握している場合は、代わりに spaps(x,y,tol) を使用することを検討します。この関数は、誤差測定を tol 以下にするという条件に従って粗さ測定ができるだけ小さくなるように p を選択します。この場合、誤差測定は通常 tol に指定した値と等しくなります。

データ型: single | double

w — 誤差測定の重み

ベクトル | cell 配列

誤差測定での誤差測定の重み w。x と同じサイズの非負のエントリのベクトルとして指定します。

誤差測定の重みベクトル w の既定値は ones(size(x)) です。

xx — 評価点

ベクトル | cell 配列

スプラインが評価される評価点。ベクトル、または多変量データの場合はベクトルの cell 配列として指定します。スプラインの評価は fnval を使用して実行されます。

データ型: single | double

出力引数

pp — pp 型のスプライン

スプラインの構造体

pp 型のスプライン。次のフィールドがある構造体として返されます。

Form — スプラインの形式

pp

スプラインの形式。pp として返されます。pp は、スプラインが区分的多項式型であることを示します。

Breaks — スプラインの節点の位置

ベクトル | cell 配列

スプラインの節点の位置。ベクトル、または多変量データの場合はベクトルの cell 配列として返されます。ベクトルには厳密に増加する要素が含まれます。これらの要素は、多項式区分が定義される各区間の開始と終了を表します。

Coefs — 多項式の係数

行列 | 配列

各区分の多項式の係数。行列、または多変量データの場合は配列として返されます。

Pieces — 多項式区分の数

スカラー | ベクトル

スプラインを記述する多項式区分の数。スカラー、または多変量データの場合は各変数の区分数のベクトルとして返されます。

Order — 多項式の次数

スカラー | ベクトル

スプラインの各多項式区分を記述する多項式関数の次数。スカラー、または多変量データの場合は各変数の次数を含むベクトルとして返されます。

Dim — 次元

スカラー

ターゲット関数の次元。スカラーとして返されます。

values — 評価されたスプライン

ベクトル | 行列 | 配列

評価されたスプライン。ベクトルまたは行列として、または多変量データの場合は配列として返されます。スプラインは特定の評価点 xx で評価されます。

P — 平滑化パラメーター

スカラー | cell 配列

スプラインの計算に使用された平滑化パラメーター。スカラー、または多変量データの場合はスカラー値の cell 配列として返されます。P は 0 と 1 の間になります。

アルゴリズム

csaps は PGS からの Fortran ルーチン SMOOTH の実装です。

平滑化スプラインの計算では、データ サイト x に応じて行列 A と B を使用した形式 p*A + (1-p)*B をとる係数行列をもつ線形システムを解くことが必要になります。p の既定値は、p*trace(A) を (1-p)*trace(B) と等しくします。

バージョン履歴

R2006a より前に導入

参考

spap2 | spaps | tpaps