28.5 多項式補間
Octave はさまざまな種類の補間を適切にサポートしており、そのほとんどは補間で説明されています。前述の章で説明されている関数の簡単な代替手段の 1 つは、単一の多項式、または区分多項式 (スプライン) をいくつかの指定されたデータ ポイントに当てはめることです。変動の大きい多項式を避けるために、データに低次の多項式を当てはめることが最もよくあります。これは通常、最小二乗の意味で多項式を当てはめる必要があることを意味し、polyfit関数が行うことはまさにそれです。
: p = polyfit (x, y, n)
: [p, s] = polyfit (x, y, n)
: [p, s, mu] = polyfit (x, y, n)
点への近似の最小二乗誤差を最小化する n次多項式p ( x )の係数を返します。 [x(:), y(:)]
nは通常、近似多項式の次数を指定する 0 以上の整数です。n が論理ベクトルの場合、対応する多項式係数を選択的に使用または無視するように強制するマスクとして使用されます。
多項式係数は行ベクトルpに返されます。出力p はpolyval、近似多項式を使用して値を推定するため に直接使用できます。
オプションの出力sは、次のフィールドを含む構造体です。
'yf'
xの各値に対する多項式の値。
'バツ'
多項式係数を計算するために使用されるヴァンデルモンド行列。
'R'
QR 分解からの三角因子 R。
'Ｃ'
スケールされていない共分散行列。正式には x' * xの逆行列に等しいが、丸め誤差の伝播を最小限に抑える方法で計算されます。
'df'
自由度。
'ノルム'
残差のノルム。
2番目の出力はpolyval、予測値の統計的誤差限界を計算するために使用できます。特に、p係数の標準偏差は次のように与えられます。
sqrt (diag (s.C)/s.df) * s.normr。
3番目の出力muが存在する場合、元のデータは中心化され、スケーリングされ、フィッティングの数値安定性が向上します。係数pは、次の多項式に関連付けられています 。
xhat = (x - mu(1)) / mu(2) where mu(1) = mean (x), and mu(2) = std (x).
例 1: 論理値nと整数値n
f = @(x) x.^2 + 5; # data-generating function x = 0:5; y = f (x); ## Fit data to polynomial A*x^3 + B*x^1 p = polyfit (x, y, logical ([1, 0, 1, 0])) ⇒ p = [ 0.0680, 0, 4.2444, 0 ] ## Fit data to polynomial using all terms up to x^3 p = polyfit (x, y, 3) ⇒ p = [ -4.9608e-17, 1.0000e+00, -3.3813e-15, 5.0000e+00 ]
See also: polyval, polyaffine, roots, vander, zscore.
単一の多項式では不十分な状況では、複数の多項式を組み合わせて使用するのが解決策です。この関数は、 splinefit区分多項式 (スプライン) をデータのセットに当てはめます。
: pp = splinefit (x, y, breaks)
: pp = splinefit (x, y, p)
: pp = splinefit (…, "periodic", periodic)
: pp = splinefit (…, "robust", robust)
: pp = splinefit (…, "beta", beta)
: pp = splinefit (…, "order", order)
: pp = splinefit (…, "constraints", constraints)
ブレーク（ノット）ブレークを含む区分 3 次スプラインをノイズ データxとyに適合させます。
xはベクトルであり、y はベクトルまたは ND 配列です。y が ND 配列の場合、 x ( j) はy (:,…,:,j)と一致します。
pはx に沿った区間の数を定義する正の整数で 、p +1 はブレークの数です。各区間の点の数は 1 以内で異なります。
オプションのプロパティperiodic は、スプラインに周期境界条件を適用するかどうかを指定する論理値です。周期の長さは です。デフォルト値は です。 max (breaks) - min (breaks)false
オプションのプロパティrobustは、外れ値データ ポイントの影響を減らすためにロバスト フィッティングを適用するかどうかを指定する論理値です。重み付き最小二乗法が 3 回反復されます。重みは、以前の残差から計算されます。外れ値識別の感度は、プロパティbetaによって制御されます。betaの値は 、0 < beta < 1 の範囲に制限されます。既定値はbeta = 1/2 です。0 に近い値は、すべてのデータに等しい重みを与えます。beta の値を増やすと、外れ値データの影響が減ります。1 に近い値は、不安定性やランク不足を引き起こす可能性があります。
適合されたスプラインは区分多項式ppとして返され、 を使用して評価できますppval。
スプラインは、次数 の多項式で構成されています。デフォルトは、次数= 3 の 3 次です。P 個の部分を持つスプラインには、P+次数の自由度があります。周期境界条件では、自由度は P に減少します。
オプションのプロパティ、constraintsは、近似に対する線形制約を指定する構造体です。構造体には"xc"、、、 "yc"の3 つのフィールドがあります"cc"。
"xc"
制約の x 位置のベクトル。
"yc"
場所xcでの制約値。デフォルトはゼロの配列です。
"cc"
係数 (行列)。デフォルトは 1 の配列です。行数は区分多項式の次数 に制限されます。
制約は、次式に従って 0次から-1次の導関数の線形結合である。
cc(1,j) * y(xc(j)) + cc(2,j) * y'(xc(j)) + ... = yc(:,...,:,j).
See also: interp1, unmkpp, ppval, spline, pchip, ppder, ppint, ppjumps.
区分多項式の構築に使用されるブレーク（またはノット）の数は、入力データxとyに存在するノイズを抑制する上で重要な要素です。これは以下の例で示されています。
x = 2 * pi * rand (1, 200);
y = sin (x) + sin (2 * x) + 0.2 * randn (size (x));
## Uniform breaks
breaks = linspace (0, 2 * pi, 41); % 41 breaks, 40 pieces
pp1 = splinefit (x, y, breaks);
## Breaks interpolated from data
pp2 = splinefit (x, y, 10); % 11 breaks, 10 pieces
## Plot
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
plot (x, y, ".", xx, [y1; y2])
axis tight
ylim auto
legend ({"data", "41 breaks, 40 pieces", "11 breaks, 10 pieces"})
その結果は図28.1に示されています。
Figure 28.1:41 個のブレークを持つ区分多項式と 11 個のブレークを持つ区分多項式のフィッティングの比較。ブレークの数が多い場合のフィッティングでは、基礎となる関数には存在しない高速リップルが見られます。
によって提供される区分多項式近似は、次数splinefit-1までの連続導関数を持ちます。たとえば、3次近似は連続した1次導関数と2次導関数を持ちます。これはコードによって示されます。
## Data (200 points)
x = 2 * pi * rand (1, 200);
y = sin (x) + sin (2 * x) + 0.1 * randn (size (x));
## Piecewise constant
pp1 = splinefit (x, y, 8, "order", 0);
## Piecewise linear
pp2 = splinefit (x, y, 8, "order", 1);
## Piecewise quadratic
pp3 = splinefit (x, y, 8, "order", 2);
## Piecewise cubic
pp4 = splinefit (x, y, 8, "order", 3);
## Piecewise quartic
pp5 = splinefit (x, y, 8, "order", 4);
## Plot
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
y3 = ppval (pp3, xx);
y4 = ppval (pp4, xx);
y5 = ppval (pp5, xx);
plot (x, y, ".", xx, [y1; y2; y3; y4; y5])
axis tight
ylim auto
legend ({"data", "order 0", "order 1", "order 2", "order 3", "order 4"})
その結果は図28.2に示されています。
Figure 28.2:ノイズの多いデータに対する 8 回のブレークを持つ区分定数、線形、2 次、3 次、4 次多項式の比較。高次のソリューションは基礎となる関数をより正確に表しますが、計算の複雑さが増します。
適合を提供する基礎関数が周期的である場合、splinefit 周期的適合を明示するために必要な境界条件を適用できます。これは以下のコードで示されます。
## Data (100 points)
x = 2 * pi * [0, (rand (1, 98)), 1];
y = sin (x) - cos (2 * x) + 0.2 * randn (size (x));
## No constraints
pp1 = splinefit (x, y, 10, "order", 5);
## Periodic boundaries
pp2 = splinefit (x, y, 10, "order", 5, "periodic", true);
## Plot
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
plot (x, y, ".", xx, [y1; y2])
axis tight
ylim auto
legend ({"data", "no constraints", "periodic"})
その結果は図28.3に示されています。
Figure 28.3: 周期境界条件がある場合とない場合のノイズの多い周期関数に対する区分多項式近似の比較。
さらに複雑な制約を追加することもできます。たとえば、以下のコードは、エンドポイントでクランプされた値を使用した周期的フィットと、エンドポイントでヒンジされた 2 番目の周期的フィットを示しています。
## Data (200 points)
x = 2 * pi * rand (1, 200);
y = sin (2 * x) + 0.1 * randn (size (x));
## Breaks
breaks = linspace (0, 2 * pi, 10);
## Clamped endpoints, y = y' = 0
xc = [0, 0, 2*pi, 2*pi];
cc = [(eye (2)), (eye (2))];
con = struct ("xc", xc, "cc", cc);
pp1 = splinefit (x, y, breaks, "constraints", con);
## Hinged periodic endpoints, y = 0
con = struct ("xc", 0);
pp2 = splinefit (x, y, breaks, "constraints", con, "periodic", true);
## Plot
xx = linspace (0, 2 * pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
plot (x, y, ".", xx, [y1; y2])
axis tight
ylim auto
legend ({"data", "clamped", "hinged periodic"})
その結果は図28.4に示されています。
Figure 28.4: ノイズの多い周期信号に対する 2 つの周期区分 3 次近似の比較。1 つの近似ではエンドポイントがクランプされ、2 つ目の近似ではエンドポイントがヒンジされています。
この関数は、外れ値データの影響が軽減される堅牢なsplinefitフィッティングの利便性も提供します 。以下の例では、3 つの異なるフィッティングが提供されています。2 つは外れ値抑制のレベルが異なり、3 つ目は堅牢でないソリューションを示しています。
## Data
x = linspace (0, 2*pi, 200);
y = sin (x) + sin (2 * x) + 0.05 * randn (size (x));
## Add outliers
x = [x, linspace(0,2*pi,60)];
y = [y, -ones(1,60)];
## Fit splines with hinged conditions
con = struct ("xc", [0, 2*pi]);
## Robust fitting, beta = 0.25
pp1 = splinefit (x, y, 8, "constraints", con, "beta", 0.25);
## Robust fitting, beta = 0.75
pp2 = splinefit (x, y, 8, "constraints", con, "beta", 0.75);
## No robust fitting
pp3 = splinefit (x, y, 8, "constraints", con);
## Plot
xx = linspace (0, 2*pi, 400);
y1 = ppval (pp1, xx);
y2 = ppval (pp2, xx);
y3 = ppval (pp3, xx);
plot (x, y, ".", xx, [y1; y2; y3])
legend ({"data with outliers","robust, beta = 0.25", ...
"robust, beta = 0.75", "no robust fitting"})
axis tight
ylim auto
その結果は図28.5に示されています。
Figure 28.5: 外れ値データと組み合わされたノイズの多いデータに対する 2 つの異なるレベルの堅牢なフィッティング (ベータ= 0.25 と 0.75)の比較。堅牢なフィッティングなしの従来のフィッティング (ベータ= 0) も含まれています。
多項式解釈の非常に特殊な形式はパデ近似です。制御システムの場合、連続時間遅延は近似値を使用して非常に簡単にモデル化できます。
: [num, den] = padecoef (T) : [num, den] = padecoef (T, N)
伝達関数形式の 連続時間遅延TのN次パデ近似を計算します。
のパデ近似はexp (-sT)次の式で定義される。
Pn(s)
exp (-sT) ~ -------
Qn(s)
ここで、Pn(s)とQn(s)はともに次の式で定義される N次の有理関数である。
N (2N - k)!N! k
Pn(s) = SUM --------------- (-sT)
k=0 (2N)!k!(N - k)!
Qn(s) = Pn(-s)
入力TとN は負でない数値スカラーでなければなりません。N が指定されていない場合は、 デフォルトで1になります。
出力行ベクトルnumとden には、分子係数と分母係数が s の降順で格納されます。どちらも N次多項式です。
例えば：
t = 0.1;
n = 4;
[num, den] = padecoef (t, n)
⇒ num =
1.0000e-04 -2.0000e-02 1.8000e+00 -8.4000e+01 1.6800e+03
⇒ den =
1.0000e-04 2.0000e-02 1.8000e+00 8.4000e+01 1.6800e+03
関数 はppval、または他の手段によって作成された区分多項式を評価しmkpp、unmkpp区分多項式に関する詳細な情報を返します。
次の例は、2 つの線形関数と 2 次関数を 1 つの関数に結合する方法を示しています。これらの各関数は、隣接する区間で表現されます。
x = [-2, -1, 1, 2];
p = [ 0, 1, 0;
1, -2, 1;
0, -1, 1 ];
pp = mkpp (x, p);
xi = linspace (-2, 2, 50);
yi = ppval (pp, xi);
plot (xi, yi);
: pp = mkpp (breaks, coefs)
: pp = mkpp (breaks, coefs, d)
サンプルポイントbreaksと係数coefsから区分多項式(pp)構造を構築します 。
ブレークは厳密に増加する値のベクトルでなければなりません。間隔の数は で与えられます。 ni = length (breaks) - 1
mが多項式の次数の場合、係数のサイズはni行 ( m + 1)列でなければなりません 。
係数の i 行目 ( )には、 i番目の区間における多項式の係数が最高次数 ( m ) から最低次数 ( 0 ) の順に並べられて格納されます。 coefs(i,:)
coefs は、ベクトル値または配列値の多項式を指定する多次元配列にすることもできます。その場合、多項式の次数m は、 coefsの最後の次元の長さによって定義されます。最初の次元のサイズは、スカラーまたはベクトルdによって指定されます。dが指定されていない場合は、 に設定されます1。この場合、 には、 区間iで定義されたr番目の多項式の係数が含まれます 。いずれの場合も、coefs は サイズ の 2 次元行列に再形成されます。 p(r, i, 
[ni*prod(d) m]
プログラミングノート:ppvalは多項式を で評価します 。つまり、現在の区間の下端点をxiから減算します。 で区分多項式オブジェクトを作成するときは、これを考慮する必要があります。 xi - breaks(i)mkpp
See also: unmkpp, ppval, spline, pchip, ppder, ppint, ppjumps.
: [x, p, n, k, d] = unmkpp (pp)
区分多項式構造ppの要素を抽出します。
この関数は の逆関数です。区分多項式構造ppmkppを作成するために必要な入力を に抽出します 。以下のコードではこの関係を明示的に示しています。 mkpp
[breaks、coefs、numinter、order、dim] = unmkpp (pp);
pp2 = mkpp (breaks, coefs, dim);
このようにして得られた区分多項式構造はpp2、元の と同一ですpp。同じものは、構造 のフィールドに直接アクセスすることによっても得られますpp。
コンポーネントは次のとおりです。
x
サンプルポイントまたはブレーク。
p
サンプル区間内の点の多項式係数。 には、最高次数から最低次数の順に並べられた区間i上の多項式の係数が含まれます。の場合、pはサイズ の行列であり 、 行は 区間i内のすべてのd個の多項式の係数です。 p(i, d > 1[n*prod(d) m]i + (1:d)
ん
多項式の部分または区間の数 。 n = length (x) - 1
け
多項式の次数に 1 を加えたもの。
d
各間隔に定義された多項式の数。
See also: mkpp, ppval, spline, pchip.
: yi = ppval (pp, xi)
点xiにおける区分多項式構造pp を評価します。
pp がスカラー多項式関数を記述する場合、結果はxiと同じ形状の配列になります。それ以外の場合、結果のサイズは、 xiがベクトルの 場合、または多次元配列の場合になります。 [pp.dim, length(xi)][pp.dim, size(xi)]
See also: mkpp, unmkpp, spline, pchip.
: ppd = ppder (pp)
: ppd = ppder (pp, m)
区分多項式構造体ppの区分m次導関数を計算します。
mが省略された場合は、1 次導関数が計算されます。
See also: mkpp, ppval, ppint.
: ppi = ppint (pp)
: ppi = ppint (pp, c)
区分多項式 struct ppの積分を計算します。
cが指定されている場合は、積分定数です。
See also: mkpp, ppval, ppder.
: jumps = ppjumps (pp)
区分多項式の境界ジャンプを評価します。
区間がn個あり、 ppの次元がdの場合 、結果の配列の次元は になります[d, n-1]。
See also: mkpp.