ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

hampel

Hampel 識別子を使用した外れ値の削除

構文

y = hampel(x)
y = hampel(x,k)
y = hampel(x,k,nsigma)
[y,j] = hampel(___)
[y,j,xmedian,xsigma] = hampel(___)
hampel(___)

説明

y = hampel(x) は、Hampel フィルターを入力ベクトル x に適用し、外れ値の検出および削除を行います。この関数は、x のサンプルごとに、そのサンプルと周囲の 6 個のサンプル (左右の各 3 個) で構成されるウィンドウの中央値を計算します。また、中央絶対偏差を使用して、ウィンドウ中央値に対する各サンプルの標準偏差も推定します。サンプルが中央値から標準偏差の 3 倍より大きく離れている場合には、中央値と置き換えられます。x が行列の場合、hampelx の各列を独立チャネルとして扱います。

y = hampel(x,k) は、x の各サンプルに対する、測定ウィンドウにおける右側と左側の各近傍の数 k を指定します。

y = hampel(x,k,nsigma) は、x のサンプルが、局所的中央値から標準偏差の何倍分離れていたら中央値と置き換えられるかを指示する標準偏差の個数 nsigma を指定します。

[y,j] = hampel(___) も、外れ値として識別されたすべての点の位置で true である logical 行列を返します。この構文は、前の構文の任意の入力引数を受け入れます。

[y,j,xmedian,xsigma] = hampel(___) も、x の各要素の局所的中央値と推定標準偏差を返します。

出力引数なしの hampel(___) は、フィルター処理した信号をプロットして、削除された外れ値に注釈を付けます。

すべて折りたたむ

正弦波信号のサンプルを 100 個生成します。6 番目と 20 番目のサンプルをスパイクに置き換えます。

x = sin(2*pi*(0:99)/100);
x(6) = 2;
x(20) = -2;

hampel を使用して、局所的中央値から 3 標準偏差を超えて離れた位置にあるすべてのサンプルを特定します。測定ウィンドウは、そのサンプルと周囲の 6 個のサンプル (片側ごとに 3 個) で構成されます。

[y,i,xmedian,xsigma] = hampel(x);

フィルター処理された信号をプロットし、外れ値に注釈を付けます。

n = 1:length(x);
plot(n,x)
hold on
plot(n,xmedian-3*xsigma,n,xmedian+3*xsigma)
plot(find(i),x(i),'sk')
hold off
legend('Original signal','Lower limit','Upper limit','Outliers')

計算を繰り返しますが、今度は中央値の計算に、両側の隣接サンプルをそれぞれ 1 つのみ使用します。関数は、極値を外れ値と見なします。

hampel(x,1)

異なる周波数の正弦波で構成される 2 チャネル信号を生成します。ランダムな位置にスパイクを配置します。NaN を使用し、欠損サンプルをランダムに追加します。再現可能な結果が必要な場合は、乱数発生器をリセットします。信号をプロットします。

rng('default')

n = 59;
x = sin(pi./[15 10]'*(1:n)+pi/3)';

spk = randi(2*n,9,1);
x(spk) = x(spk)*2;
x(randi(2*n,6,1)) = NaN;

plot(x)

hampel を既定の設定で使用して、信号をフィルター処理します。

y = hampel(x);
plot(y)

移動ウィンドウの長さを増やし、サンプルを外れ値として扱うしきい値を下げます。

y = hampel(x,4,2);
plot(y)

各チャネルの移動中央値を出力します。信号のプロットに中央値を重ね合わせます。

[y,j,xmd,xsd] = hampel(x,4,2);
plot(x)
hold on
plot(xmd,'--')

単位分散のホワイト ガウス ノイズ内の、周波数が異なる 2 つの正弦波から構成されるマルチチャネル信号を生成します。

rng('default')

t = 0:60;
x = sin(pi./[10;2]*t)'+randn(numel(t),2);

信号に Hampel フィルターを適用します。周囲 9 サンプルのウィンドウの中央値から 2 標準偏差を超えて離れている点を外れ値とします。外れ値の位置で true となる logical 行列を出力します。

k = 4;
nsig = 2;

[y,h] = hampel(x,k,nsig);

信号の各チャネルを、それ自体の座標軸セットにプロットします。元の信号、フィルター処理した信号および外れ値を描画します。外れ値の位置に注釈を付けます。

for k = 1:2
    hk = h(:,k);
    ax = subplot(2,1,k);
    plot(t,x(:,k))
    hold on
    plot(t,y(:,k))
    plot(t(hk),x(hk,k),'*')
    hold off
    ax.XTick = t(hk);
end

正弦波信号のサンプルを 100 個生成します。6 番目と 20 番目のサンプルをスパイクに置き換えます。

n = 1:100;
x = sin(2*pi*n/100);
x(6) = 2;
x(20) = -2;

hampel を使用して、すべてのサンプルに対する局所的中央値と推定標準偏差を計算します。入力パラメーターの既定値を使用します。

  • ウィンドウ サイズは です。

  • ウィンドウの中央値から 3 標準偏差を超えて離れている点は外れ値と見なされます。

結果をプロットします。

[y,i,xmedian,xsigma] = hampel(x);

plot(n,x)
hold on
plot(n,[1;1]*xmedian+3*[-1;1]*xsigma)
plot(find(i),x(i),'sk')
hold off
legend('Signal','Lower','Upper','Outliers')

のウィンドウ サイズと、外れ値の識別基準として 2 標準偏差を使用して、計算を繰り返します。

sds = 2;
adj = 10;
[y,i,xmedian,xsigma] = hampel(x,adj,sds);

plot(n,x)
hold on
plot(n,[1;1]*xmedian+sds*[-1;1]*xsigma)
plot(find(i),x(i),'sk')
hold off
legend('Signal','Lower','Upper','Outliers')

入力引数

すべて折りたたむ

ベクトルまたは行列として指定される入力信号。x が行列の場合、hampelx の各列を独立チャネルとして扱います。

例: cos(pi/4*(0:159))+randn(1,160) は単一チャネルの行ベクトル信号です。

例: cos(pi./[4;2]*(0:159))'+randn(160,2) は 2 チャネル信号です。

データ型: single | double

サンプル xs のそれぞれの側の近傍の数。整数値スカラーで指定します。片側のサンプル数が k より少ない信号エッジの近傍のサンプルは、より小さいウィンドウの中央値と比較されます。

データ型: single | double

x のサンプルが、サンプルの局所的中央値から標準偏差の何倍分離れていたら外れ値と見なされるかを指示する標準偏差の個数。nsigma を実数スカラーで指定します。この関数は、局所的な中央絶対偏差 (MAD) を κ=12erfc11/21.4826 の係数でスケーリングすることによって標準偏差を推定します。

データ型: single | double

出力引数

すべて折りたたむ

x と同じサイズのベクトルまたは行列として返されるフィルター処理された信号。

データ型: single | double

外れ値のインデックス。x と同じサイズのベクトルまたは行列として返されます。

データ型: logical

局所的な中央値。x と同じサイズのベクトルまたは行列として返されます。

データ型: single | double

推定標準偏差。x と同じサイズのベクトルまたは行列として返されます。

データ型: single | double

詳細

すべて折りたたむ

Hampel 識別子

Hampel 識別子は、外れ値に対してロバストな、統計の 3 シグマ ルールの一種です。

たとえば、シーケンス x1, x2, x3, …, xn と長さ k のスライディング ウィンドウの場合、ポイントツーポイントの中央値と標準偏差推定は以下を使用して定義します。

  • 局所的中央値 — mi=median(xik,xik+1,xik+2,,xi,,xi+k2,xi+k1,xi+k)

  • 標準偏差 — σi=κmedian(|xikmi|,,|xi+kmi|), ここで κ=12erfc11/21.4826

σi /κ  の量は中央絶対偏差 (MAD) として知られています。

サンプル xi が、与えられたしきい値 nσ に対して以下の条件を満たすとします。

|ximi|>nσσi

このとき、Hampel 識別子は xi を外れ値として宣言し、mi で置き換えます。

シーケンス端点付近で、mi および σi の計算に使用するウィンドウが関数によって打ち切られます。

  • i < k + 1

    mi=median(x1,x2,x3,,xi,,xi+k2,xi+k1,xi+k)

    σi=κmedian(|x1m1|,,|xi+kmi|)

  • i > n – k

    mi=median(xik,xik+1,xik+2,,xi,,xn2,xn1,xn)

    σi=κmedian(|xikmi|,,|xnmn|)

参照

[1] Liu, Hancong, Sirish Shah, and Wei Jiang. “On-line outlier detection and data cleaning.” Computers and Chemical Engineering. Vol. 28, March 2004, pp. 1635–1647.

R2015b で導入