robustcov

ロバスト多変量共分散および平均の推定

構文

sig = robustcov(x)

[sig,mu] = robustcov(x)

[sig,mu,mah] = robustcov(x)

[sig,mu,mah,outliers] = robustcov(x)

[sig,mu,mah,outliers,s] = robustcov(x)

[___] = robustcov(x,Name,Value)

説明

sig = robustcov(x) は、x に格納されている多変量データのロバスト共分散推定値 sig を返します。

[sig,mu] = robustcov(x) は、ロバスト最小共分散行列式 (MCD) の平均 mu の推定値も返します。

[sig,mu,mah] = robustcov(x) は、平均および共分散のロバスト推定値を使用して観測値のマハラノビス距離として計算したロバスト距離 mah も返します。

[sig,mu,mah,outliers] = robustcov(x) は、標本データ内で外れ値として記録されている観測値のインデックス outliers も返します。

[sig,mu,mah,outliers,s] = robustcov(x) は、推定に関する情報が格納されている構造体 s も返します。

[___] = robustcov(x,Name,Value) は、1 つ以上の Name,Value ペア引数で指定された追加オプションを使用して、前の構文に示されている引数のいずれかを返します。たとえば、使用するロバスト推定器や、アトラクターに対して使用する開始方法を指定できます。

例

距離-距離プロットの使用による外れ値の検出

ガウス型コピュラを使用して、二変量分布からランダムなデータ点を生成します。

rng default
rho = [1,0.05;0.05,1];
u = copularnd('Gaussian',rho,50);

無作為に選択した 5 つの観測値を外れ値に変更します。

noise = randperm(50,5);
u(noise,1) = u(noise,1)*5;

使用可能な 3 つの方式を使用して、ロバスト共分散行列を計算します。つまり、高速 MCD、直交 Gnanadesikan-Kettenring (OGK) および Olive-Hawkins です。

[Sfmcd, Mfmcd, dfmcd, Outfmcd] = robustcov(u);
[Sogk, Mogk, dogk, Outogk] = robustcov(u,'Method','ogk');
[Soh, Moh, doh, Outoh] = robustcov(u,'Method','olivehawkins');

マハラノビス尺度を使用して、標本データの従来型の距離の値を計算します。

d_classical = pdist2(u, mean(u),'mahal');
p = size(u,2);
chi2quantile = sqrt(chi2inv(0.975,p));

各ロバスト共分散計算法について DD プロットを作成します。

tiledlayout(2,2)
nexttile
plot(d_classical, dfmcd, 'o')
line([chi2quantile, chi2quantile], [0, 30], 'color', 'r')
line([0, 6], [chi2quantile, chi2quantile], 'color', 'r')
hold on
plot(d_classical(Outfmcd), dfmcd(Outfmcd), 'r+')
xlabel('Mahalanobis Distance')
ylabel('Robust Distance')
title('DD Plot, FMCD method')
hold off

nexttile
plot(d_classical, dogk, 'o')
line([chi2quantile, chi2quantile], [0, 30], 'color', 'r')
line([0, 6], [chi2quantile, chi2quantile], 'color', 'r')
hold on
plot(d_classical(Outogk), dogk(Outogk), 'r+')
xlabel('Mahalanobis Distance')
ylabel('Robust Distance')
title('DD Plot, OGK method')
hold off

nexttile
plot(d_classical, doh, 'o')
line([chi2quantile, chi2quantile], [0, 30], 'color', 'r')
line([0, 6], [chi2quantile, chi2quantile], 'color', 'r')
hold on
plot(d_classical(Outoh), doh(Outoh), 'r+')
xlabel('Mahalanobis Distance')
ylabel('Robust Distance')
title('DD Plot, Olive-Hawkins method')
hold off

DD プロットでは、原点を通過する直線にデータ点が集まる傾向があります。この直線から離れた位置にある点は、一般に外れ値であると考えられます。前の各プロットで、赤い '+' 記号は robustcov で外れ値と見なされたデータ点を示します。

多変量正規分布のデータの評価

この例では、robustcov を使用して標本データが多変量正規分布または他の楕円輪郭 (EC) 分布であるかを評価する方法を示します。

多変量正規分布から、ランダムな標本データを生成します。ロバスト共分散推定用 (Olive-Hawkins 法を使用) と従来型の共分散推定用のマハラノビス距離を計算します。

rng('default')
x1 = mvnrnd(zeros(1,3),eye(3),200);
[~, ~, d1] = robustcov(x1,'Method','olivehawkins');
d_classical1 = pdist2(x1,mean(x1),'mahalanobis');

楕円輪郭 (EC) 分布から、ランダムな標本データを生成します。ロバスト共分散推定用 (Olive-Hawkins 法を使用) と従来型の共分散推定用のマハラノビス距離を計算します。

mu1 = [0 0 0];
sig1 = eye(3);
mu2 = [0 0 0];
sig2 = 25*eye(3);
x2 = [mvnrnd(mu1,sig1,120);mvnrnd(mu2,sig2,80)];
[~, ~, d2] = robustcov(x2, 'Method','olivehawkins');
d_classical2 = pdist2(x2, mean(x2), 'mahalanobis');

多変量対数正規分布から、ランダムな標本データを生成します。これは、多変量正規分布でも楕円輪郭分布でもありません。ロバスト共分散推定用 (Olive-Hawkins 法を使用) と従来型の共分散推定用のマハラノビス距離を計算します。

x3 = exp(x1);
[~, ~, d3] = robustcov(x3, 'Method','olivehawkins');
d_classical3 = pdist2(x3, mean(x3), 'mahalanobis');

比較のため、3 組の標本データのそれぞれについて D-D プロットを作成します。

figure
subplot(2,2,1)
plot(d_classical1,d1, 'o')
line([0 4.5], [0, 4.5])
xlabel('Mahalanobis Distance')
ylabel('Robust Distance')
title('DD Plot, Multivariate Normal')

subplot(2,2,2)
plot(d_classical2, d2, 'o')
line([0 18], [0, 18])
xlabel('Mahalanobis Distance')
ylabel('Robust Distance')
title('DD Plot, Elliptically-Contoured')

subplot(2,2,3)
plot(d_classical3, d3, 'o')
line([0 18], [0, 18])
xlabel('Mahalanobis Distance')
ylabel('Robust Distance')
title('DD Plot, 200 Lognormal cases')

多変量正規分布のデータ (左上) では、プロットされた点は原点から伸びる 45°の直線に従っています。楕円輪郭分布のデータ (右上) では、プロットされた点は直線に従っていますが、角度は 45°ではありません。対数正規分布 (左下) では、プロットされた点は直線に従っていません。

ほとんどの点がプロットの左下にあるので、対数正規分布のパターンを特定することは困難です。重み付きの DD プロットを使用して、この隅を拡大し、大きいロバスト距離が存在すると明らかではなくなる特徴量を明らかにします。

d3_weighted = d3(d3 < sqrt(chi2inv(0.975,3)));
d_classical_weighted = d_classical3(d3 < sqrt(chi2inv(0.975,3)));

4 番目のサブプロットを図に追加して、対数正規分布データに重みを付けた結果を示します。

subplot(2,2,4)
plot(d_classical_weighted, d3_weighted, 'o')
line([0 3], [0, 3])
xlabel('Mahalanobis Distance')
ylabel('Robust Distance')
title('Weighted DD Plot, 200 Lognormal cases')

このプロットのスケールは、対数正規データが元の DD プロットの拡大表示になっていることを示しています。この表示は、プロットにパターンがないことをより明確に示しています。したがって、データが多変量正規分布にも楕円輪郭分布にもなっていないことがわかります。

ロバスト共分散の計算と外れ値のプロット

ガウス型コピュラを使用して、二変量分布からランダムなデータ点を生成します。

rng default
rho = [1,0.05;0.05,1];
u = copularnd('Gaussian',rho,50);

無作為に選択した 5 つの観測値を外れ値に変更します。

noise = randperm(50,5);
u(noise,1) = u(noise,1)*5;

散布図を使用して二変量データを可視化します。

figure
scatter(u(:,1),u(:,2))

ほとんどのデータ点がプロットの左側にあります。しかし、一部のデータ点は右の離れた位置にあります。これらの点は、共分散行列の計算に影響を与える外れ値である可能性があります。

従来型の共分散行列とロバスト共分散行列を比較します。

c = cov(u)  

c = 2×2

    0.5523    0.0000
    0.0000    0.0913

rc = robustcov(u)

rc = 2×2

    0.1117    0.0364
    0.0364    0.1695

標本データ内に存在する外れ値が結果に影響を与えるので、従来型の共分散行列とロバスト共分散行列は異なります。

robustcov が外れ値と見なすデータ点を特定してプロットします。

[sig,mu,mah,outliers] = robustcov(u);
figure
gscatter(u(:,1),u(:,2),outliers,'br','ox')
legend({'Not outliers','Outliers'})

robustcov は、プロットの右側にあるデータ点を潜在的な外れ値として識別し、ロバスト共分散行列を計算するときはこの結果に基づいてこれらのデータ点を扱います。

入力引数

x — 標本データ

数値の行列

ロバスト共分散行列の推定に使用する標本データ。数値の行列を指定します。x は、各行が観測値に、各列が変数に対応する n 行 p 列の行列です。

robustcov は、ロバスト共分散行列を計算するときに、予測子の値が欠損している行を削除します。

データ型: single | double

名前と値の引数

オプションの引数のペアを Name1=Value1,...,NameN=ValueN として指定します。ここで Name は引数名、Value は対応する値です。名前と値の引数は他の引数の後ろにする必要がありますが、ペアの順序は関係ありません。

R2021a より前では、名前と値をそれぞれコンマを使って区切り、Name を引用符で囲みます。

例: 'Method','ogk','NumOGKIterations',1 は、ロバスト推定器として直交 Gnanadesikan-Kettenring 法を指定し、直交化の反復回数を 1 に設定します。

すべての推定器

Method — ロバスト推定器

'fmcd' (既定値) | 'ogk' | 'olivehawkins'

ロバスト推定器。次のいずれかを指定します。

名前	値
'fmcd'	高速 MCD (最小共分散行列式) 法
'ogk'	直交 Gnanadesikan-Kettenring (OGK) 推定
'olivehawkins'	集中アルゴリズム手法 (高速で高い整合性をもつ、外れ値に対してロバストな一連の方式)

例: 'Method','ogk'

FMCD 法および OliveHawkins 法のみ

OutlierFraction — 外れ値の比率

0.5 (既定値) | 範囲 [0,0.5] の数値

外れ値の比率。'OutlierFraction' と 範囲 [0,0.5] の数値から構成されるコンマ区切りのペアとして指定します。共分散行列式を最小化する対象となる観測値の比率は、値 1 – OutlierFraction によって指定されます。

このアルゴリズムでは、サイズが h = ceiling(n + p + 1) / 2) である副標本が選択されます。ここで、n は観測値の個数、p は次元数です。OutlierFraction は、分解が最大になる値であり、共分散行列式を最小化する対象となるサブセット h のサイズを調節します。その後、サブセットごとにほぼ (1 – OutlierFraction) × n 個の観測値になるように h が選択されます。

例: 'OutlierFraction',0.25

データ型: single | double

NumTrials — 試行回数

正の整数値

試行回数。'NumTrials' と正の整数値から構成されるコンマ区切りのペアとして指定します。

'Method' が 'fmcd' である場合、NumTrials はアルゴリズムで開始点として標本データから無作為抽出される、サイズが p + 1 の副標本の個数です。p は標本データの次元数です。この場合、NumTrials の既定値は 500 です。

'Method' が 'olivehawkins' の場合、NumTrials は使用する試験近似、つまりアトラクターの数です。この場合、NumTrials の既定値は 2 です。このオプションは、非決定的な開始の場合のみ役立ちます。

例: 'NumTrials',300

データ型: single | double

FMCD 法のみ

BiasCorrection — 小規模標本補正係数を適用するためのフラグ

1 (既定値) | 0

小規模標本補正係数を適用するためのフラグ。'BiasCorrection' と 1 または 0 から構成されるコンマ区切りのペアとして指定します。値 1 は、小規模な標本の場合に robustcov が共分散推定のバイアスを補正することを示します。値 0 は、この補正を robustcov が適用しないことを示します。

例: 'BiasCorrection',0

データ型: logical

OGK 法のみ

NumOGKIterations — 直交化の反復回数

2 (既定値) | 正の整数値

直交化の反復回数。'NumOGKIterations' と正の整数値から構成されるコンマ区切りのペアとして指定します。通常、この値は 1 または 2 に設定します。ステップ数を増やしても、推定の精度が向上する可能性はほとんどありません。

例: 'NumIter',1

データ型: single | double

UnivariateEstimator — 一変量ロバスト推定を計算する関数

'tauscale' (既定値) | 'qn'

一変量ロバスト推定を計算する関数。'UnivariateEstimator' と次のいずれかから構成されるコンマ区切りのペアとして指定します。

名前	値
'tauscale'	Yohai および Zamar の "tau スケール" 推定を使用します。これは、切り捨てられた標準偏差と加重平均です。
'qn'	Croux および Rousseeuw の Qn スケール推定を使用します。

例: 'UnivariateEstimator','qn'

OliveHawkins 法のみ

ReweightingMethod — 再重み付けの方式

'rfch' (既定値) | 'rmvn'

効率性ステップにおける再重み付けの方式。'ReweightingMethod' と次のいずれかから構成されるコンマ区切りのペアとして指定します。

名前	値
'rfch'	2 つの再重み付けステップを使用します。これは、効率を向上させるための再重み付けの標準的な方式です。
'rmvn'	多変量正規を再重み付けします。クリーンなデータが多変量正規分布の場合、さまざまな外れ値構成の下で真の共分散行列を推定する場合に役立つ、2 つの再重み付けステップを使用します。

例: 'ReweightingMethod','rmvn'

NumConcentrationSteps — 集中ステップの数

10 (既定値) | 正の整数値

集中ステップの数。'NumConcentrationSteps' と正の整数値から構成されるコンマ区切りのペアとして指定します。

例: 'NumConcentrationSteps',8

データ型: single | double

StartMethod — 各アトラクターへの開始方法

'classical' (既定値) | 'medianball' | 'elemental' | 関数ハンドル | cell 配列

各アトラクターへの開始方法。'Start' と次のいずれかから構成されるコンマ区切りのペアとして指定します。

名前	値
'classical'	従来型の推定器を開始点として使用します。これは、単独で使用した場合は DGK 推定器として知られる DGK アトラクターです。
'medianball'	メディアン ボール (MB) を開始点として使用します。メディアン ボールは (med(x),eye(p)) です。したがって、MB の開始点を計算するため、標本の中央値からユークリッド距離で最も遠いデータが 50% のケースでトリミングされます。これは、単独で使用した場合は MB 推定器として知られる MB アトラクターです。
'elemental'	アトラクターは集中によって生成され、その開始点は無作為に選択された基本の開始点です。無作為に選択された p + 1 個のケースの "基本セット" に従来型の推定器が適用されます。この "基本" アトラクターは、計算効率は高くなりますが、整合性が低く分解がゼロであるという、理論的な欠点があります。

既定の設定では、アトラクターは次のように選択されます。アトラクターのいずれかが 'medianball' である場合、median(X) から推定位置までのユークリッド距離がデータの半数より大きい (つまり、メディアン ボールの外部にある) アトラクターは使用されません。そして、MCD の基準に基づいて最終的なアトラクターが選択されます。

初期位置と散乱の推定値を計算するための 2 つの出力引数を返す関数の関数ハンドルを指定することもできます。

前の表のオプションと関数ハンドルの任意の組み合わせが格納されている cell 配列を指定することもできます。cell 配列の長さと等しい数のアトラクターが使用されます。この方法を使用すると、アルゴリズムをより細かく調節でき、アトラクターおよび開始点の数を自由に指定できます。

例: 'StartMethod','medianball'

出力引数

sig — ロバスト共分散行列の推定値

数値行列

ロバスト共分散行列の推定値。p 行 p 列の数値行列として返されます。p は標本データに含まれている予測子の数です。

mu — ロバスト平均の推定値

数値の配列

ロバスト平均の推定値。1 行 p 列の数値配列として返されます。p は標本データに含まれている予測子の数です。

mah — ロバストなマハラノビス距離

数値の配列

ロバストなマハラノビス距離。1 行 n 列の数値配列として返されます。robustcov は欠損データが含まれている行を x から除外するので、mah の行数が x の行数より少なくなる場合があります。

outliers — 外れ値のインデックス

論理値の配列

標本データ x で外れ値として記録された観測値のインデックス。1 行 n 列の論理値配列として返されます。値 0 は、観測値が外れ値ではないことを示します。値 1 は、観測値が外れ値であることを示します。

robustcov は欠損データが含まれている行を x から除外するので、outliers の行数が x の行数より少なくなる場合があります。

s — 推定情報が格納される構造体

構造体

推定情報が格納される構造体。構造体として返されます。

詳細

マハラノビス距離

マハラノビス距離は、標本点と分布の間の尺度です。

ベクトル x から平均 μ および共分散 Σ をもつ分布までのマハラノビス距離は次のようになります。

$$d=\sqrt{(x-\mu ){\sum }^{-1}(x-\mu )\text{'}}.$$

この距離は、標準偏差単位で x が平均からどの程度離れているかを表します。

robustcov は、x の観測値から平均 mu および共分散 sig をもつ分布までのロバストなマハラノビス距離 (mah) を返します。

アルゴリズム

最小共分散行列式推定

"最小共分散行列式" (MCD) は、多変量の位置および散乱の最も高速な推定器であり、整合性が高くロバストです。ただし、可能な標本データのサブセットをすべて評価すると計算時間が非常に長くなるので、MCD を正確に評価することは困難です。robustcov は高速 MCD 法を使用して MCD を実装します [3]。

高速 MCD 法では、行列式が最小である従来型の共分散行列をもつ n 個の観測値から h 個の観測値を選択します (n/2 < h ≤ n)。MCD の平均は、選択した h 個の観測値の平均です。

MCD の共分散は、選択した h 個の点の共分散行列に、多変量正規分布で整合性を得るための一致係数と、小さい標本サイズのバイアスを補正するための補正係数を乗算した行列です。

直交 Gnanadesikan-Kettenring 推定

"直交 Gnanadesikan-Kettenring" (OGK) 推定は、Gnanadesikan-Kettenring (GK) 推定器 (非正定値である可能性があるペアワイズのロバストな散乱行列) から始まる、散乱の正定値推定です [1]。この推定では、固有値 (負の可能性があります) をロバストな分散に置き換えて、直交反復と呼ばれる主成分の形式をペアワイズの散乱行列に対して使用します。この手順は結果を改善するために繰り返すことができ、通常は 2 ～ 3 回の反復後に収束します。

Olive Hawkins 推定

Olive-Hawkins 推定では、Olive および Hawkins が提案した "集中アルゴリズム" 手法を使用します。これは、高速で高い整合性をもつ、外れ値に対して非常にロバストな一連の方式です。この推定は、4 次のモーメントをもつ楕円輪郭分布の共分散の、ロバストな root-n-consistent 推定量です。この推定は、まず試験推定 (開始点) を生成し、それぞれの試験近似からの集中手法を使用したアトラクターの取得によって得られます。

(T_0j,C_0j) が開始点であるとすると、次回の反復における従来型の平均と共分散の推定量は、前回の反復からの推定に基づくマハラノビス距離が最小である約 n / 2 個のケース (n は観測値の個数) から計算されます。この反復は固定回数 (k 回) のステップについて続けることができ、最終ステップ (k 回目) の推定がアトラクターになります。最終的な推定は、与えられた基準に基づいて選択されます。

既定の設定では、2 つのアトラクターが使用されます。1 番目のアトラクターは Devlin-Gnanadesikan-Kettering (DGK) アトラクターです。使用される開始点は従来型の推定器です。2 番目のアトラクターはメディアン ボール (MB) アトラクターです。使用される開始点は (median(x),eye(p))、つまりユークリッド距離で median(x) に最も近い半分のデータです。DGK アトラクターの位置推定器がメディアン ボールの外部にある場合は MB アトラクターが使用され、それ以外の場合は行列式が最も小さいアトラクターが使用されます。最終的な推定平均は、選択したアトラクターの推定平均です。最終的な推定共分散は、選択したアトラクターの推定共分散に対し、正規分布において推定を一致させるスケーリング係数を乗算することにより求められます。

参照

[1] Maronna, R. and Zamar, R.H.. “Robust estimates of location and dispersion for high dimensional datasets.” Technometrics, Vol. 50, 2002.

[2] Pison, S. Van Aelst and G. Willems. “Small Sample Corrections for LTS and MCD.” Metrika, Vol. 55, 2002.

[3] Rousseeuw, P.J. and Van Driessen, K. “A fast algorithm for the minimum covariance determinant estimator.” Technometrics, Vol. 41, 1999.

[4] Olive, D.J. “A resistant estimator of multivariate location and dispersion.” Computational Statistics and Data Analysis, Vol. 46, pp. 99–102, 2004.

拡張機能

スレッドベースの環境

MATLAB® の backgroundPool を使用してバックグラウンドでコードを実行するか、Parallel Computing Toolbox™ の ThreadPool を使用してコードを高速化します。

この関数は、スレッドベースの環境を完全にサポートします。詳細については、スレッドベースの環境での MATLAB 関数の実行を参照してください。

バージョン履歴

R2016a で導入

参考

cov | mahal | iforest | rrcforest | lof | ocsvm

トピック

• 教師なし異常検出