robustcov
ロバスト多変量共分散および平均の推定
構文
説明
[___] = robustcov(
は、1 つ以上の x
,Name,Value
)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'
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
BiasCorrection
— 小規模標本補正係数を適用するためのフラグ
1
(既定値) | 0
小規模標本補正係数を適用するためのフラグ。'BiasCorrection'
と 1
または 0
から構成されるコンマ区切りのペアとして指定します。値 1
は、小規模な標本の場合に robustcov
が共分散推定のバイアスを補正することを示します。値 0
は、この補正を robustcov
が適用しないことを示します。
例: 'BiasCorrection',0
データ型: logical
NumOGKIterations
— 直交化の反復回数
2 (既定値) | 正の整数値
直交化の反復回数。'NumOGKIterations'
と正の整数値から構成されるコンマ区切りのペアとして指定します。通常、この値は 1 または 2 に設定します。ステップ数を増やしても、推定の精度が向上する可能性はほとんどありません。
例: 'NumIter',1
データ型: single
| double
UnivariateEstimator
— 一変量ロバスト推定を計算する関数
'tauscale'
(既定値) | 'qn'
一変量ロバスト推定を計算する関数。'UnivariateEstimator'
と次のいずれかから構成されるコンマ区切りのペアとして指定します。
名前 | 値 |
---|---|
'tauscale' | Yohai および Zamar の "tau スケール" 推定を使用します。これは、切り捨てられた標準偏差と加重平均です。 |
'qn' | Croux および Rousseeuw の Qn スケール推定を使用します。 |
例: 'UnivariateEstimator','qn'
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 は標本データに含まれている予測子の数です。
outliers
— 外れ値のインデックス
論理値の配列
標本データ x
で外れ値として記録された観測値のインデックス。1 行 n 列の論理値配列として返されます。値 0
は、観測値が外れ値ではないことを示します。値 1
は、観測値が外れ値であることを示します。
robustcov
は欠損データが含まれている行を x
から除外するので、outliers
の行数が x
の行数より少なくなる場合があります。
s
— 推定情報が格納される構造体
構造体
推定情報が格納される構造体。構造体として返されます。
詳細
マハラノビス距離
マハラノビス距離は、標本点と分布の間の尺度です。
ベクトル x から平均 μ および共分散 Σ をもつ分布までのマハラノビス距離は次のようになります。
この距離は、標準偏差単位で 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 推定量です。この推定は、まず試験推定 (開始点) を生成し、それぞれの試験近似からの集中手法を使用したアトラクターの取得によって得られます。
(T0j,C0j) が開始点であるとすると、次回の反復における従来型の平均と共分散の推定量は、前回の反復からの推定に基づくマハラノビス距離が最小である約 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 で導入
MATLAB コマンド
次の MATLAB コマンドに対応するリンクがクリックされました。
コマンドを MATLAB コマンド ウィンドウに入力して実行してください。Web ブラウザーは MATLAB コマンドをサポートしていません。
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list:
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)