最小二乗近似

はじめに

Curve Fitting Toolbox™ ソフトウェアでは、データを近似するときに最小二乗法を使用します。近似では、1 つ以上の係数を使用して応答データを予測子データと関連付けるパラメトリックモデルが必要です。近似プロセスで得られる結果はモデル係数の推定値です。

最小二乗法では、係数推定値を得るために残差平方和を最小化します。i 番目のデータ点 r_iの残差は、観測応答値 y_iと近似応答値 ŷ_iの差として定義され、データに関連する誤差として認識されます。

$\begin{array}{l} r_{i} = y_{i} - {\hat{y}}_{i} \\ residual=data - fit \end{array}$

残差平方和は次で与えられます

$S = \sum_{i = 1}^{n} r_{i}^{2} = \sum_{i = 1}^{n} {(y_{i} - {\hat{y}}_{i})}^{2}$

ここで、n は近似に含まれるデータ点の数、S は誤差の二乗和の推定値です。以下のタイプの最小二乗近似がサポートされています。

線形最小二乗法
重み付き線形最小二乗法
ロバスト最小二乗法
非線形最小二乗法

誤差分布

ランダムな変動を含むデータを近似する場合、誤差に関して通常適用される重要な仮定として次の 2 つがあります。

誤差は応答データのみに存在し、予測子データには存在しない。
誤差はランダムであり、平均がゼロで分散 σ²が一定の正規 (ガウス) 分布に従う。

2 番目の仮定は多くの場合、次のように表現されます。

$e r r o r \sim N (0, σ^{2})$

正規分布は通常多くの測定量の分布に対して適切な近似となるため、誤差は正規分布すると仮定します。最小二乗近似法は、パラメーター推定値を計算するときに誤差が正規分布していることを仮定しませんが、極端な値の確率的誤差が大量に含まれていないデータに対しては最も効果的な方法です。正規分布は極端な確率的誤差が少ない確率分布の 1 つです。ただし、信頼限界や予測限界などの統計結果では、妥当性を確保するために誤差が正規分布している必要があります。

誤差の平均値がゼロの場合、誤差は純粋にランダムです。平均値がゼロでない場合、データに対するモデルの選択が適切でないか、誤差が純粋にランダムではなく系統的誤差を含んでいる可能性があります。

データの分散が一定であるのは、誤差の "広がり" が一定であることを意味します。分散が同じデータは、"等質"と呼ばれる場合があります。

重み付き最小二乗回帰では,確率的誤差の分散が一定であることを暗黙的に仮定しません。その代わりに,近似手順で指定される重みによってデータの品質レベルの違いが正しく示されることを仮定します。この重みは、近似係数の推定値に対する各データ点の影響の大きさを適切なレベルに調整するために使用されます。

線形最小二乗法

Curve Fitting Toolbox ソフトウェアでは、線形最小二乗法を使用して線形モデルをデータに当てはめます。"線形"モデルは、係数について線形である方程式として定義されます。たとえば、多項式は線形ですが、ガウス関数は線形ではありません。線形最小二乗近似のプロセスを説明するために、1 次多項式でモデル化できる n 個のデータ点があるとします。

$y = p_{1} x + p_{2}$

未知の係数 p₁と p₂についてこの方程式を解くために、S を未知数が 2 つの n 個の方程式から成る連立線形方程式系とします。n が未知数の数より大きい場合、この方程式系は "過決定"です。

$S = \sum_{i = 1}^{n} {(y_{i} - (p_{1} x_{i} + p_{2}))}^{2}$

最小二乗近似プロセスでは残差平方和を最小化するため、S を各パラメーターで微分し、その結果をゼロとすることにより係数を決定します。

$\begin{array}{l} \frac{\partial S}{\partial p_{1}} = - 2 \sum_{i = 1}^{n} x_{i} (y_{i} - (p_{1} x_{i} + p_{2})) = 0 \\ \frac{\partial S}{\partial p_{2}} = - 2 \sum_{i = 1}^{n} (y_{i} - (p_{1} x_{i} + p_{2})) = 0 \end{array}$

真のパラメーターの推定値は通常、b で表されます。b₁および b₂を p₁および p₂に代入すると、前の方程式は次のようになります。

$\begin{matrix} \sum x_{i} (y_{i} - (b_{1} x_{i} + b_{2})) = 0 \\ \sum (y_{i} - (b_{1} x_{i} + b_{2})) = 0 \end{matrix}$

ここで、和の範囲は i = 1 から n までです。"正規方程式"は次のように定義されます

$\begin{matrix} b_{1} \sum x_{i}^{2} + b_{2} \sum x_{i} = \sum x_{i} y_{i} \\ b_{1} \sum x_{i} + n b_{2} = \sum y_{i} \end{matrix}$

b₁について解きます

$b_{1} = \frac{n \sum x_{i} y_{i} - \sum x_{i} \sum y_{i}}{n \sum x_{i}^{2} - {(\sum x_{i})}^{2}}$

b₁の値を使用して b₂について解きます

$b_{2} = \frac{1}{n} (\sum y_{i} - b_{1} \sum x_{i})$

このように、係数 p₁および p₂の推定に必要なのは、数回の簡単な計算だけです。この例をより次数の高い多項式に拡張するのは簡単ですが、多少手間がかかります。必要なのは、モデルに追加した各線形項に対する正規方程式を追加することだけです。

行列形式では、線形モデルは次の式で与えられます

y = Xβ + ε

ここで、

y は n 行 1 列の応答ベクトルです。
β は m 行 1 列の係数ベクトルです。
X はこのモデルの n 行 m 列の計画行列です。
ε は n 行 1 列の誤差ベクトルです。

1 次多項式の場合、未知数が 2 つの n 個の方程式は y、X および β により次のように表されます。

$[\begin{matrix} y_{1} \\ y_{2} \\ y_{3} \\ . \\ . \\ . \\ y_{n} \end{matrix}] = [\begin{matrix} x_{1} 1 \\ x_{2} 1 \\ x_{3} 1 \\ . \\ . \\ . \\ x_{n} 1 \end{matrix}] \times [\begin{array}{l} p_{1} \\ p_{2} \end{array}]$

この問題の最小二乗解をベクトル b とします。これにより未知の係数ベクトル β が推定されます。正規方程式は次で与えられます

(X^TX)b = X^Ty

ここで、X^Tは計画行列 X の転置行列です。b について解きます

b = (X^TX)^–1X^Ty

MATLAB^®のバックスラッシュ演算子 (mldivide) を使用して、未知の係数について連立線形方程式系を解きます。X^TX の逆行列を計算すると許容できない丸め誤差が発生する場合があるため、バックスラッシュ演算子ではピボットによる QR 分解が使用されます。この分解は数値的に非常に安定なアルゴリズムです。バックスラッシュ演算子と QR 分解の詳細については、算術演算を参照してください。

b をモデルの式に戻すと予測された応答値 ŷ が得られます。

ŷ = Xb = Hy

H = X(X^TX)^–1X^T

文字の上にあるハット記号 (曲折アクセント符号) はパラメーターの推定値またはモデルの予測値を示します。射影行列 H は y にハット記号を付けるためハット行列と呼ばれます。

残差は次で与えられます

r = y – ŷ = (1–H)y

重み付き最小二乗法

通常、応答データは等質、つまり分散が一定であると仮定します。この仮定を満たしていない場合、近似は低品質のデータの影響を過度に受ける可能性があります。この近似を改良するために重み付き最小二乗回帰を使用できます。この場合、追加のスケール係数 (重み) が近似プロセスに含まれます。重み付き最小二乗回帰は次の誤差の推定値を最小化します

$s = \sum_{i = 1}^{n} w_{i} {(y_{i} - {\hat{y}}_{i})}^{2}$

ここで、w_iは重みです。この重みは、各応答値が最終的なパラメーター推定値にどの程度影響するかを決定します。品質の高いデータ点が、品質の低いデータ点よりも近似に大きく影響するようにします。重みが既知の場合またはそれらが特定の形式に従う正当な理由がある場合は、データに重みを付けることを推奨します。

重みを付けると、パラメーター推定値 b の式は次のように変わります。

$b = \hat{β} = {(X^{T} W X)}^{- 1} X^{T} W y$

ここで W は重み行列 w の対角要素で与えられます。

多くの場合、分散が一定かどうかはデータを近似して残差をプロットすることで判断できます。以下に示すプロットでは、データにさまざまな品質の複製データが含まれていて、近似は正しいと考えられます。残差のプロットからデータの品質が低いことがわかります。このプロットは "漏斗" 状であり、大きい予測子値よりも小さい予測子値の方で応答値が大きく散在しています。

重みは、応答の分散が定数値になるように指定する必要があります。データに含まれる測定誤差の分散が既知の場合、重みは次で与えられます

$w_{i} = 1 / σ_{i}^{2}$

また、各データ点の誤差変数の推定値しかわからない場合は、真の分散の代わりにこれらの推定値を使用することで通常は間に合います。分散が不明な場合、相対スケールの重みを指定すれば十分です。重みを指定しても全体の分散項を推定できることに注意してください。この場合、重みは近似の各点の相対的な重みを定義しますが、それぞれの点の正確な分散を指定するものとは見なされません。

たとえば、各データ点が複数の独立した測定値の平均値である場合、測定回数を重みとして使用するのが理にかなっていることがあります。

ロバスト最小二乗法

通常、応答誤差は正規分布に従い、極端な値はまれであると仮定します。ただし、"外れ値"と呼ばれる極端な値が発生します。

最小二乗近似の主な欠点は外れ値の影響を受けやすいことです。残差を 2 乗すると極端なデータ点の効果が拡大するため、外れ値は近似に大きく影響します。外れ値の影響を最小限に抑えるには、ロバスト最小二乗回帰を使用してデータに当てはめます。ツールボックスでは、次の 2 つのロバスト回帰法が用意されています。

最小絶対残差 (LAR) — LAR 法は、残差の二乗ではなく残差の絶対値を最小化する曲線を求めます。そのため、極端な値による近似への影響が小さくなります。
二重平方重み — この方法は、重み付き二乗和を最小化します。各データ点に与えられる重みはその点の近似直線からの距離によって決まります。直線に近い点は最大限の重み付けがなされます。直線から遠い点ほど重みが減少します。直線からの距離が偶然に起因すると考えられる距離を超えている点はゼロに重み付けされます。
二重平方重みは通常の最小二乗法を使用してデータの大部分に当てはめる曲線を求めると同時に外れ値の影響を最小化するため、多くの場合、LAR よりも二重平方重み法の方が好まれます。

二重平方重みによるロバスト近似は、反復重み付き最小二乗アルゴリズムを使用し、次の手順に従います。

重み付き最小二乗法によるモデルで近似します。
"調整済み残差"を計算し,それらを標準化します。調整済み残差は次で与えられます

$r_{a d j} = \frac{r_{i}}{\sqrt{1 - h_{i}}}$

r_iは通常の最小二乗残差、h_iは "てこ比"であり、最小二乗近似に大きな影響を与えるてこ比の大きいデータ点の重みを減らすことにより残差を調整します。標準化された調整済み残差は次で与えられます

$u = \frac{r_{a d j}}{K s}$

K は調整定数 4.685、s は MAD/0.6745 で与えられるロバスト標準偏差です。ここで MAD は残差の中央絶対偏差です。
ロバスト重みを u の関数として計算します。二重平方重みは次で与えられます

$w_{i} = {\begin{matrix} {(1 - {(u_{i})}^{2})}^{2} & | u_{i} | < 1 \\ 0 & | u_{i} | \geq 1 \end{matrix}$

独自の回帰重みベクトルを指定する場合、最終的な重みはロバスト重みと回帰重みの積であることに注意してください。
近似が収束した場合、これで完了です。そうでない場合は、最初の手順に戻り、近似手順の次の反復を実行します。

次に示すプロットは二重平方重みを使用するロバスト近似を通常の線形近似と比較しています。ロバスト近似は大部分のデータに従い、外れ値の影響をあまり受けていないことに注意してください。

ロバスト回帰を使用して外れ値の影響を最小化する代わりに、近似から排除するデータ点をマークすることができます。詳細については、外れ値の削除を参照してください。

非線形最小二乗法

Curve Fitting Toolbox ソフトウェアでは、非線形最小二乗法を使用して非線形モデルをデータに当てはめます。非線形モデルは、係数について非線形な方程式、または係数について線形と非線形の組み合わせである方程式として定義されます。たとえば、ガウス関数、多項式の比、べき関数はすべて非線形です。

行列形式では、非線形モデルは次の式で与えられます

y = f(X,β) + ε

ここで、

y は n 行 1 列の応答ベクトルです。
f は β と X の関数です。
β は m 行 1 列の係数ベクトルです。
X はこのモデルの n 行 m 列の計画行列です。
ε は n 行 1 列の誤差ベクトルです。

非線形モデルでは、単純な行列手法を使用した係数推定ができないため、線形モデルよりも近似が困難です。その代わりに、次の手順に従う反復法が必要になります。

各係数の初期推定から始めます。一部の非線形モデルでは、妥当な開始値を生成するヒューリスティックな方法が用意されています。その他のモデルでは、区間 [0,1] でランダムな値が提供されます。
現在の一連の係数について近似曲線を生成します。近似応答値 ŷ は次で与えられます

ŷ = f(X,b)

ここで f(X,b) の "ヤコビアン"を計算する必要があります。これは、係数についての偏導関数の行列として定義されます。
係数を調整し、近似を改良できるかどうかを判断します。調整の方向と大きさは近似アルゴリズムにより異なります。ツールボックスには次のアルゴリズムが用意されています。
- Trust-Region — これは既定のアルゴリズムであり、係数制約を指定する場合は使用しなければなりません。難しい非線形問題を他のアルゴリズムよりも効率的に解くことができ、一般に普及しているレーベンバーグ・マルカートアルゴリズムよりも改善が見られます。
- レーベンバーグ・マルカート — このアルゴリズムは長年にわたって使用されており、さまざまな非線形モデルと開始値に対して大体の場合は機能することが実証されています。信頼領域アルゴリズムで適切な近似が生成されず、係数制約がない場合は、レーベンバーグ・マルカートアルゴリズムを試してください。
指定された収束条件を近似が満たすまで、手順 2 に戻りこのプロセスを繰り返します。

非線形モデルに重みおよびロバスト近似を使用できますが、近似プロセスはそれに応じて変更されます。

この近似プロセスの性質上、すべての非線形モデル、データセット、開始点に対応する絶対確実なアルゴリズムはありません。そのため、既定の開始点、アルゴリズム、収束条件を使用して妥当な近似を実現できない場合は、別のオプションを試してください。既定のオプションを変更する方法については、近似オプションと最適化された開始点の指定を参照してください。非線形モデルは開始点の影響を特に受けやすい場合があるため、近似オプションの中でも開始点を最初に変更してください。

ロバスト近似

ライブスクリプトを開く

この例では、外れ値の排除による効果とロバスト近似の効果を比較する方法を示します。この例は、モデルからの距離が標準偏差の 1.5 倍を超える任意の外れ値を排除する方法を示します。さらに、その手順において、外れ値を削除する場合と外れ値に小さい重みを設定するロバスト近似を指定する場合を比較します。

ベースラインの正弦波信号を作成します。

xdata = (0:0.1:2*pi)'; y0 = sin(xdata);

分散が一定ではないノイズを信号に追加します。

% Response-dependent Gaussian noisegnoise = y0.*randn(size(y0));% Salt-and-pepper noisespnoise = zeros(size(y0)); p = randperm(length(y0)); sppoints = p(1:round(length(p)/5)); spnoise(sppoints) = 5*sign(y0(sppoints)); ydata = y0 + gnoise + spnoise;

ノイズを含むデータにベースラインの正弦波モデルを当てはめ、3 つの出力引数を指定して残差を含む近似情報を取得します。

f = fittype('a*sin(b*x)'); [fit1,gof,fitinfo] = fit(xdata,ydata,f,'StartPoint',[1 1]);

fitinfo 構造体の情報を確認します。

fitinfo

fitinfo =struct with fields:numobs: 63 numparam: 2 residuals: [63x1 double] Jacobian: [63x2 double] exitflag: 3 firstorderopt: 0.0883 iterations: 5 funcCount: 18 cgiterations: 0 algorithm: 'trust-region-reflective' stepsize: 4.1539e-04 message: 'Success, but fitting stopped because change in residuals less than tolerance (TolFun).'

fitinfo 構造体から残差を取得します。

residuals = fitinfo.residuals;

ベースラインモデルからの距離が標準偏差の 1.5 倍を超える任意の点を "外れ値" と特定し、外れ値を排除したデータを再近似します。

I = abs( residuals) > 1.5 * std( residuals ); outliers = excludedata(xdata,ydata,'indices',I); fit2 = fit(xdata,ydata,f,'StartPoint',[1 1],...'Exclude',outliers);

外れ値を排除する効果とロバスト近似で外れ値に小さい二重平方重みを設定する効果を比較します。

fit3 = fit(xdata,ydata,f,'StartPoint',[1 1],'Robust','on');

データ、外れ値、近似の結果をプロットします。内容を説明する凡例を指定します。

plot(fit1,'r-',xdata,ydata,'k.',outliers,'m*') holdonplot(fit2,'c--') plot(fit3,'b:') xlim([0 2*pi]) legend('Data','Data excluded from second fit','Original fit',...'Fit with points excluded','Robust fit') holdoff

Figure contains an axes object. The axes object contains 5 objects of type line. These objects represent Data, Data excluded from second fit, Original fit, Fit with points excluded, Robust fit.

外れ値を考慮した 2 つの近似の残差をプロットします。

figure plot(fit2,xdata,ydata,'co','residuals') holdonplot(fit3,xdata,ydata,“软”,'residuals') holdoff

Figure contains an axes object. The axes object contains 4 objects of type line. These objects represent data, zero line.