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 が行列の場合、この関数は x の各列を独立チャネルとして扱います。

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

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

[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)

マルチチャネル信号の Hampel フィルター処理

異なる周波数の正弦波で構成される 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

Hampel フィルターによって返される信号統計

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

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

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

• ウィンドウ サイズは $$2\times 3+1=7$$ です。

• ウィンドウの中央値から 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\times 10+1=21$$ のウィンドウ サイズと、外れ値の識別基準として 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 — 入力信号

ベクトル | 行列

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

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

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

データ型: single | double

k — それぞれの側の近傍の数

3 (既定値) | 整数スカラー

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

データ型: single | double

nsigma — 標準偏差の数

3 (既定値) | 実数スカラー

x のサンプルが、サンプルの局所的中央値から標準偏差の何倍分離れていたら外れ値と見なされるかを指示する標準偏差の個数。nsigma を実数スカラーで指定します。この関数は、局所的な中央絶対偏差 (MAD) を $$\kappa =\frac{1}{\sqrt{2}{\mathrm{erf}}^{-1}\left(1/2\right)}\approx 1.4826$$ の係数でスケーリングすることによって標準偏差を推定します。

データ型: single | double

出力引数

y — フィルター処理された信号

ベクトル | 行列

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

データ型: single | double

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

ベクトル | 行列

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

データ型: logical

xmedian — 局所的中央値

ベクトル | 行列

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

データ型: single | double

xsigma — 推定標準偏差

ベクトル | 行列

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

データ型: single | double

詳細

Hampel 識別子

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

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

• 局所的中央値 — $$m_i=\mathrm{median}\left(x_{i-k},x_{i-k+1},x_{i-k+2},\dots ,x_i,\dots ,x_{i+k-2},x_{i+k-1},x_{i+k}\right)$$

• 標準偏差 — $${\sigma }_i=\kappa \mathrm{median}\left(\left|x_{i-k}-m_i\right|,\dots ,\left|x_{i+k}-m_i\right|\right),$$ ここで、$$\kappa =\frac{1}{\sqrt{2}{\mathrm{erf}}^{-1}\left(1/2\right)}\approx 1.4826$$ です。

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

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

$$\left|x_i-m_i\right|>n_{\sigma }{\sigma }_i$$

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

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

• i < k + 1

  $$m_i=\mathrm{median}\left(x_1,x_2,x_3,\dots ,x_i,\dots ,x_{i+k-2},x_{i+k-1},x_{i+k}\right)$$

  $${\sigma }_i=\kappa \mathrm{median}\left(\left|x_1-m_1\right|,\dots ,\left|x_{i+k}-m_i\right|\right)$$

• i > n – k

  $$m_i=\mathrm{median}\left(x_{i-k},x_{i-k+1},x_{i-k+2},\dots ,x_i,\dots ,x_{n-2},x_{n-1},x_n\right)$$

  $${\sigma }_i=\kappa \mathrm{median}\left(\left|x_{i-k}-m_i\right|,\dots ,\left|x_n-m_n\right|\right)$$

参照

[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.

[2] Suomela, Jukka. “Median Filtering Is Equivalent to Sorting.” 2014.

拡張機能

C/C++ コード生成

MATLAB® Coder™ を使用して C および C++ コードを生成します。

使用上の注意および制限:

• 可変サイズの入力行列または N 次元配列は実行時に行ベクトルに変換されてはいけません。

• プロットは MATLAB^® の実行でのみサポートされています。

バージョン履歴

R2015b で導入

参考

medfilt1 | median | filloutliers | filter | isoutlier | mad (Statistics and Machine Learning Toolbox) | movmad | movmedian | sgolayfilt

トピック

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