robustfit

ロバスト線形回帰の当てはめ

構文

b = robustfit(X,y)

b = robustfit(X,y,wfun,tune,const)

[b,stats] = robustfit(___)

説明

b = robustfit(X,y) は、行列 X 内の予測子に対するベクトル y 内の応答についてのロバスト多重線形回帰の係数推定値が格納されているベクトル b を返します。

b = robustfit(X,y,wfun,tune,const) は、近似重み関数オプション wfun および tune と、モデルが定数項を含むかどうかを決定するインジケーター const を指定します。wfun、tune および const について [] を渡すと、既定値を使用できます。

[b,stats] = robustfit(___) は、前の構文におけるいずれかの入力引数の組み合わせを使用して、推定統計量を格納する構造体 stats も返します。

例

ロバスト回帰係数の推定

複数の線形モデルについて、ロバスト回帰係数を推定します。

carsmall データ セットを読み込みます。予測子として自動車の重量と馬力を、応答としてガロンあたりのマイル数による燃費を指定します。

load carsmall
x1 = Weight;
x2 = Horsepower;
X = [x1 x2];
y = MPG;

ロバスト回帰係数を計算します。

b = robustfit(X,y)

b = 3×1

   47.1975
   -0.0068
   -0.0333

近似モデルをプロットします。

x1fit = linspace(min(x1),max(x1),20);
x2fit = linspace(min(x2),max(x2),20);
[X1FIT,X2FIT] = meshgrid(x1fit,x2fit);
YFIT = b(1) + b(2)*X1FIT + b(3)*X2FIT;
mesh(X1FIT,X2FIT,YFIT)

データをプロットする。

hold on
scatter3(x1,x2,y,'filled')
hold off
xlabel('Weight')
ylabel('Horsepower')
zlabel('MPG')
legend('Model','Data')
view(50,10)
axis tight

ロバスト重み関数の調整

異なる調整定数を使用して、ロバスト回帰の重み関数を調整します。

トレンド $$y=10-2x$$ でデータを生成し、1 つの値を変更して外れ値をシミュレートします。

x = (1:10)';
rng ('default') % For reproducibility
y = 10 - 2*x + randn(10,1);
y(10) = 0;

3 つの異なる調整定数について、二重平方重み関数を使用してロバスト回帰残差を計算します。既定の調整定数は 4.685 です。

tune_const = [3 4.685 6];

for i = 1:length(tune_const)
    [~,stats] = robustfit(x,y,'bisquare',tune_const(i));
    resids(:,i) = stats.resid;
end

残差のプロットを作成します。

scatter(x,resids(:,1),'b','filled')
hold on
plot(resids(:,2),'rx','MarkerSize',10,'LineWidth',2)
scatter(x,resids(:,3),'g','filled')
plot([min(x) max(x)],[0 0],'--k')
hold off 
grid on
xlabel('x')
ylabel('Residuals')
legend('tune = 3','tune = 4.685','tune = 6','Location','best')

3 つの異なる調整定数について、残差の平方根平均二乗誤差 (RMSE) を計算します。

rmse = sqrt(mean(resids.^2))

rmse = 1×3

    3.2577    2.7576    2.7099

調整定数を大きくすると、外れ値に割り当てられる重みの削減量が減るため、調整定数が大きくなるほど RMSE は小さくなります。

ロバスト回帰と最小二乗回帰の比較

トレンド $$y=10-2x$$ でデータを生成し、1 つの値を変更して外れ値をシミュレートします。

x = (1:10)';
rng('default') % For reproducibility
y = 10 - 2*x + randn(10,1);
y(10) = 0;

通常の最小二乗回帰を使用して直線を近似します。定数項があるモデルの係数推定値を計算するには、1 の列を x に含めます。

bls = regress(y,[ones(10,1) x])

bls = 2×1

    7.8518
   -1.3644

ロバスト回帰を使用して直線の近似を推定します。robustfit は、既定で定数項をモデルに追加します。

[brob,stats] = robustfit(x,y);
brob

brob = 2×1

    8.4504
   -1.5278

残差と残差の中央絶対偏差を比較することにより、潜在的な外れ値を特定します。

outliers_ind = find(abs(stats.resid)>stats.mad_s);

ロバスト回帰の残差の棒グラフをプロットします。

bar(abs(stats.resid))
hold on
yline(stats.mad_s,'k--')
hold off
xlabel('x')
ylabel('Residuals')

データの散布図を作成します。

scatter(x,y,'filled')

外れ値をプロットします。

hold on 
plot(x(outliers_ind),y(outliers_ind),'mo','LineWidth',2)

最小二乗とロバスト近似をプロットします。

plot(x,bls(1)+bls(2)*x,'r')
plot(x,brob(1)+brob(2)*x,'g')
hold off
xlabel('x')
ylabel('y')
legend('Data','Outlier','Ordinary Least Squares','Robust Regression')
grid on

ロバスト近似に対する外れ値の影響は、最小二乗近似よりも小さくなります。

入力引数

X — 予測子データ

数値行列

予測子データ。n 行 p 列の数値行列を指定します。X の行は観測値に、列は予測子変数に対応します。X の行数は y と同じでなければなりません。

既定では、const を 'off' として指定することによって明示的に削除しない限り、robustfit は定数項をモデルに追加します。そのため、1 の列を X に含めないでください。

データ型: single | double

y — 応答データ

数値ベクトル

応答データ。n 行 1 列の数値ベクトルを指定します。y の各行は異なる観測値に対応します。y の行数は X と同じでなければなりません。

データ型: single | double

wfun — ロバスト近似重み関数

'bisquare' (既定値) | 'andrews' | 'cauchy' | 'fair' | 関数ハンドル | ...

ロバスト近似重み関数。次の表に記載されている重み関数、または関数ハンドルの名前として指定します。robustfit は、tune によって特に指定されない限り、対応する既定の調整定数を使用します。

重み関数	説明	既定の設定の調整定数
'andrews'	w = (abs(r)<pi) .* sin(r) ./ r	1.339
'bisquare'	w = (abs(r)<1) .* (1 - r.^2).^2 (biweight とも呼ばれます)	4.685
'cauchy'	w = 1 ./ (1 + r.^2)	2.385
'fair'	w = 1 ./ (1 + abs(r))	1.400
'huber'	w = 1 ./ max(1, abs(r))	1.345
'logistic'	w = tanh(r) ./ r	1.205
'ols'	通常最小二乗 (重み関数なし)	なし
'talwar'	w = 1 * (abs(r)<1)	2.795
'welsch'	w = exp(-(r.^2))	2.985
関数ハンドル	スケーリングされた残差のベクトル r を受け入れ、r と同じサイズの重みのベクトルを返す、カスタム重み関数	1

重み関数の値 r は次のようになります。

r = resid/(tune*s*sqrt(1–h)),

ここで

• resid は、前回の反復の残差のベクトルです。

• tune は、調整定数です。

• h は、最小二乗近似のてこ比値のベクトルです。

• s は、s = MAD/0.6745 によって与えられる誤差項の標準偏差の推定値です。

MAD は、残差の中央値に対する残差の中央絶対偏差です。定数 0.6745 は、推定を正規分布に対して不偏にします。X に p 個の列がある場合、中央値を計算するときに、最小のものから p 個の絶対偏差が除外されます。

データ型: char | string | function handle

tune — 調整定数

正のスカラー

調整定数。正のスカラーとして指定します。tune を設定しない場合、robustfit は、各重み関数に対応する既定の調整定数を使用します (wfun の表を参照)。

応答が正規分布に従っており、外れ値がない場合、組み込み重み関数の既定の調整定数を使用すると、統計的な効率が通常の最小二乗推定の約 95% である係数推定値が得られます。調整定数を小さくすると、大きな残差に割り当てられる重みの削減量が増え、調整定数を大きくすると、大きな残差に割り当てられる重みの削減量が減ります。

データ型: single | double

const — 定数項のインジケーター

'on' (既定値) | 'off'

当てはめにおける定数項のインジケーター。'on' または 'off' として指定します。const が 'on' である場合、robustfit は 1 列の最初の列を予測子行列 X に追加し、出力 b は (p + 1) 行 1 列のベクトルになります。const が 'off' である場合、X は変更されず、b は p 行 1 列のベクトルになります。

データ型: char | string

出力引数

b — ロバスト多重線形回帰の係数推定値

数値ベクトル

ロバスト多重線形回帰の係数推定値。数値ベクトルとして返されます。b は p 行 1 列のベクトルです。p は X 内の予測子の数です。

既定では、const を 'off' として指定することによって明示的に削除しない限り、robustfit は定数項をモデルに追加します。

stats — モデルの統計量

構造体

モデルの統計量。構造体として返されます。次の表は、ロバスト回帰の診断統計量の構造体に関するフィールドについて説明しています。

フィールド	説明
ols_s	通常の最小二乗からのシグマ推定 (平方根平均二乗誤差)
robust_s	シグマのロバスト推定
mad_s	中央値からの残差の中央絶対偏差を使用して計算されたシグマの推定。反復近似中の残差のスケーリングに使用される
s	シグマの最終推定。robust_s と、ols_s および robust_s の加重平均の中で大きい値
resid	残差。観測した値から近似した値を減算した値 (生の残差を参照)
rstud	スチューデント化残差。残差を残差標準偏差の独立した推定値で除算した値 (スチューデント化残差を参照)
se	推定される係数値 b の標準誤差
covb	係数推定値の推定共分散行列
coeffcorr	係数推定値の推定相関
t	各係数の t 統計量。モデル内の他の予測子が与えられた場合に、対応する係数が 0 ではないという対立仮説に対して係数が 0 であるという帰無仮説を検定するために使用します。t = b/se であることに注意してください。
p	対応する係数が 0 であるかどうかについての仮説検定の t 統計量に対する p 値
w	ロバスト近似の重みベクトル
R	X の QR 分解の R 因子
dfe	誤差 (残差) の自由度。推定された係数の数を観測値の数から減算した値に等しくなる
h	最小二乗近似のてこ比値のベクトル

詳細

てこ比

てこ比は、入力空間で特定の観測の位置が原因で発生した、回帰予測におけるその観測値の影響を測定します。

観測値 i のてこ比はハット行列 H の i 番目の対角項 h_ii の値です。ハット行列 H は、データ行列 X に関して次のように定義されます。

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

ハット行列は "射影行列" とも呼ばれます。これは、観測値のベクトル y を予測値のベクトル $$\widehat{y}$$ に射影するので "ハット" が y の上に置かれるためです。

てこ比値の合計は p (回帰モデルの係数の個数) なので、てこ比が p/n (n は観測値の個数) を大幅に超える場合、観測値 i は外れ値であると考えることができます。

詳細は、ハット行列とてこ比を参照してください。

ヒント

• robustfit は、X または y 内の NaN 値を欠損値として扱います。robustfit は、欠損値がある観測値をロバスト近似から省略します。

アルゴリズム

• robustfit は、反復的に再重み付けした最小二乗を使用して、係数 b を計算します。入力 wfun は重みを指定します。

• robustfit は、式 inv(X'*X)*stats.s^2 を使用して、係数推定値 stats.covb の分散共分散行列を推定します。この推定により、標準誤差 stats.se と相関 stats.coeffcorr が生成されます。

• 線形モデルにおいて、y の観測値およびその残差は確率変数です。残差はゼロ平均を伴う正規分布ですが、予測子値によって分散が異なります。残差を比較可能なスケールにするため、robustfit は残差を "スチューデント化" します。つまり robustfit は、残差の値に依存しない標準偏差の推定値で残差を除算します。スチューデント化残差は、既知の自由度をもつ t 分布に従います。robustfit は、stats.rstud のスチューデント化残差を返します。

代替機能

robustfit は、関数の出力引数のみが必要である場合、またはループ内でモデルの当てはめを複数回繰り返す場合に便利です。ロバストを当てはめた回帰モデルをさらに調べる必要がある場合は、fitlm を使用して線形回帰モデル オブジェクト LinearModel を作成します。名前と値のペアの引数 'RobustOpts' の値を 'on' に設定します。

参照

[1] DuMouchel, W. H., and F. L. O'Brien. “Integrating a Robust Option into a Multiple Regression Computing Environment.” Computer Science and Statistics: Proceedings of the 21st Symposium on the Interface. Alexandria, VA: American Statistical Association, 1989.

[2] Holland, P. W., and R. E. Welsch. “Robust Regression Using Iteratively Reweighted Least-Squares.” Communications in Statistics: Theory and Methods, A6, 1977, pp. 813–827.

[3] Huber, P. J. Robust Statistics. Hoboken, NJ: John Wiley & Sons, Inc., 1981.

[4] Street, J. O., R. J. Carroll, and D. Ruppert. “A Note on Computing Robust Regression Estimates via Iteratively Reweighted Least Squares.” The American Statistician. Vol. 42, 1988, pp. 152–154.

バージョン履歴

R2006a より前に導入

参考

LinearModel | fitlm | regress | robustdemo

トピック

• 線形回帰モデルとは
• ロバスト回帰を使用した外れ値の影響の低減