データのフィルター処理と平滑化

データのフィルター処理と平滑化について

このトピックでは、この関数を使用して応答データを平滑化する方法を説明します。関数 smooth では、オプションの方法として、移動平均、Savitzky-Golay フィルター、重みやロバスト性を使用するまたは使用しない局所回帰 (lowess、loess、rlowess および rloess) を使用できます。行列、table、timetable や移動中央値法およびガウス法のサポートを含む、その他の機能については smoothdata を参照してください。

移動平均フィルター

移動平均フィルターは、所定の範囲内で定義された隣接データ点の平均で各データ点を置き換えることによりデータを平滑化します。このプロセスは、次の差分方程式で与えられる平滑化の応答を使用するローパスフィルター処理と等価です。

$y_{s} (i) = \frac{1}{2 N + 1} (y (i + N) + y (i + N - 1) + ... + y (i - N))$

ここで、y_s(i) は i 番目のデータ点の平滑化値、N は y_s(i) の片側の隣接データ点の数、2N+1 は範囲です。

Curve Fitting Toolbox™ で使用する移動平均平滑化法は以下の規則に従います。

範囲は奇数でなければなりません。
平滑化するデータ点は範囲の中央になければなりません。
片側に指定された数の隣接点が存在しないデータ点では、範囲が調整されます。
端点では範囲を定義できないため、端点は平滑化されません。

関数 filter を使用すると、上に示したような差分方程式を実装できることに注意してください。ただし、端点の扱い方が原因で、ツールボックスの移動平均の結果は filter が返す結果と異なります。詳細については、差分方程式とフィルター処理を参照してください。

たとえば、範囲が 5 の移動平均フィルターを使用してデータを平滑化することを考えます。上で説明した規則を使用すると、y_s の最初の 4 つの要素は次で与えられます。

y_s(1) = y(1)
y_s(2) = (y(1)+y(2)+y(3))/3
y_s(3) = (y(1)+y(2)+y(3)+y(4)+y(5))/5
y_s(4) = (y(2)+y(3)+y(4)+y(5)+y(6))/5

y_s(1)、y_s(2)、...、y_s(end) は並べ替えた後のデータの順序を示し、必ずしも元の順序ではないことに注意してください。

生成されたデータセットの最初の 4 つのデータ点に関する平滑化値と範囲を次に示します。

プロット (a) は、範囲を構成できないため最初のデータ点が平滑化されていないことを示します。プロット (b) は、2 番目のデータ点が範囲 3 を使用して平滑化されていることを示します。プロット (c) と (d) は、範囲 5 を使用して平滑化値を計算していることを示します。

Savitzky-Golay フィルター

Savitzky-Golay フィルターは一般化された移動平均と考えられます。フィルター係数を導くために、指定した次数の多項式を使用する重みなし線形最小二乗近似を実行します。このため、Savitzky-Golay フィルターは、デジタル平滑化多項式フィルターまたは最小二乗平滑化フィルターとも呼ばれます。高次多項式を使用するとデータの特徴を損なわずに高いレベルの平滑化を実現できることに注意してください。

Savitzky-Golay フィルター法は多くの場合、周波数データまたは分光 (ピーク) データに使用されます。周波数データの場合、この手法は信号の高周波数成分の維持に効果があります。分光データの場合、この手法は線幅などピークの高次モーメントの維持に効果があります。これに対して、移動平均フィルターは、信号の高周波数成分からかなりの部分を除去する傾向があり、重心などピークの低次モーメントのみを維持できます。ただし、Savitzky-Golay フィルターは移動平均に比べるとノイズをあまり適切に抑制できません。

Curve Fitting Toolbox ソフトウェアで使用する Savitzky-Golay 平滑化法は以下の規則に従います。

範囲は奇数でなければなりません。
多項式の次数は範囲より小さくなければなりません。
データ点が等間隔である必要はありません。
通常、Savitzky-Golay フィルターでは予測子データが等間隔である必要があります。ただし、Curve Fitting Toolbox のアルゴリズムでは、等間隔でないデータをサポートしています。そのため、追加のフィルター処理手順を実行して等間隔のデータを作成する必要はありません。

以下に示すプロットには、生成されたガウスデータと Savitzky-Golay 法を使用したいくつかの平滑化の結果を示しています。このデータにはかなりのノイズが含まれており、ピーク幅は広いものから狭いものまでさまざまです。範囲はデータ点数の 5% に相当します。

プロット (a) にはノイズの多いデータを示しています。平滑化した結果を比較しやすくするために、プロット (b) と (c) には余分なノイズのないデータを示しています。

プロット (b) には 2 次多項式による平滑化の結果を示しています。この手法は幅の狭いピークでパフォーマンスが悪化することに注意してください。プロット (c) には 4 次多項式による平滑化の結果を示しています。一般的に、高次多項式は、幅の狭いピークの高さと幅を正確に捉えますが、幅の広いピークを適切に平滑化できません。

Lowess と Loess

"lowess" および "loess" という名前は "locally weighted scatter plot smooth (局所的に重み付けされた散布図平滑化)" から来ており、どちらの手法も局所的に重み付けされた線形回帰を使用してデータを平滑化します。

平滑化プロセスは、移動平均法のように、所定の範囲内で定義された隣接データ点によって平滑化した値が決まるため、局所的であると見なされます。所定の範囲に含まれるデータ点に対して回帰重み関数が定義されるため、このプロセスは重み付けを伴います。回帰重み関数に加え、ロバスト重み関数を使用して外れ値に対する耐性をこのプロセスにもたせることができます。さらに、この手法は回帰に使用するモデルによって区別されます。lowess は線形多項式を使用し、loess は 2 次多項式を使用します。

Curve Fitting Toolbox ソフトウェアで使用する局所回帰平滑化法は以下の規則に従います。

範囲は偶数でも奇数でもかまいません。
データセットの合計データ点数の割合として範囲を指定できます。たとえば、範囲を 0.1 とすると 10% のデータ点が使用されます。

局所回帰法

局所回帰平滑化プロセスは、各データ点について以下の手順に従います。

範囲内の各データ点で "回帰重み" を計算します。この重みは、次に示す Tricube 関数で与えられます。
$w_{i} = {(1 - {| \frac{x - x_{i}}{d (x)} |}^{3})}^{3}$
x は平滑化される応答値に関連する予測子値、x_i は範囲によって定義される x の最近傍、d(x) は横座標に沿った x から範囲内で最も遠くにある予測子値までの距離です。重みには次の特性があります。
- 平滑されるデータ点は重みが最大で、近似に最も影響します。
- 範囲外にあるデータ点は重みがゼロで、近似には影響しません。
重み付き線形最小二乗回帰を実行します。lowess の場合は、回帰に 1 次多項式を使用します。loess の場合は、回帰に 2 次多項式を使用します。
平滑化値は、関心のある予測子値での加重回帰によって与えられます。

平滑化されるデータ点の両側で同じ数の隣接データ点について平滑化計算が行われる場合、重み関数は対称です。ただし、隣接点の数が、平滑化されるデータ点に関して対称でない場合、重み関数は対称ではありません。移動平均平滑化プロセスとは異なり、範囲が変化しないことに注意してください。たとえば、予測子値が最小であるデータ点を平滑化する場合、重み関数の形状は半分に切り取られた形状になり、範囲内の左端のデータ点の重みが最大で、すべての隣接点が平滑化値の右になります。

範囲が 31 個のデータ点である場合について端点と内側の点の重み関数を次に示します。

範囲が 5 の lowess を使用した場合について、生成されたデータセットの最初の 4 つのデータ点に関する平滑化値と関連する回帰を次に示します。

平滑化プロセスが各データ点を順に進んでも範囲は変化しないことに注意してください。ただし、最近傍点の数に応じて、回帰重み関数が平滑化されるデータ点に関して対称でない場合があります。特に、プロット (a) と (b) では非対称な重み関数が使用され、プロット (c) と (d) では対称な重み関数が使用されています。

loess 法の場合、2 次多項式によって平滑化値が生成されることを除き、グラフは同じようになります。

ロバスト局所回帰

データに外れ値が含まれている場合、平滑化値がゆがみ、大部分の隣接データ点の振る舞いが反映されない可能性があります。この問題を解決するために、ごく一部の外れ値からの影響を受けないロバスト手順を使用してデータを平滑化できます。外れ値の説明については、残差分析を参照してください。

Curve Fitting Toolbox ソフトウェアでは、lowess 平滑化法と loess 平滑化法の両方のロバストバージョンが用意されています。これらのロバスト法には、外れ値に耐性のあるロバスト重みの計算が追加で含まれています。ロバスト平滑化手順は以下の手順に従います。

前セクションで説明した平滑化手順で残差を計算します。
範囲内の各データ点で "ロバスト重み" を計算します。この重みは、二重平方関数で与えられます。
$w_{i} = {\begin{matrix} {(1 - {(r_{i} / 6 M A D)}^{2})}^{2}, & | r i | < 6 M A D, \\ 0, & | r i | \geq 6 M A D, \end{matrix}$
ここで、r_i は回帰平滑化手順で生成された i 番目のデータ点の残差、MAD は残差の中央絶対偏差です。
$M A D = median (| r |) .$
中央絶対偏差は残差の広がりの尺度です。r_i が 6MAD より小さい場合、ロバスト重みは 1 に近づきます。r_i が 6MAD より大きい場合、ロバスト重みは 0 であり、関連するデータ点は平滑化計算から排除されます。
ロバスト重みを使用してデータを再度平滑化します。最終的な平滑化値は、局所回帰重みとロバスト重みの両方を使用して計算されます。
前の 2 つの手順を合計 5 回繰り返します。

生成されたデータセットが 1 つの外れ値を含む場合について、 lowess 手順の平滑化結果とロバスト lowess 手順の結果を以下で比較します。どちらの手順の範囲もデータ点数は 11 です。

プロット (a) は、いくつかの最近傍の平滑化値に外れ値が影響することを示しています。プロット (b) は、外れ値の残差が中央絶対偏差の 6 倍よりも大きいことを示しています。そのため、このデータ点でロバスト重みはゼロになります。プロット (c) は、外れ値の近くの平滑化値が大部分のデータを反映していることを示しています。

例: データの平滑化

count.dat のデータを読み込みます。

load count.dat

24 行 3 列の配列 count には、3 つの交差点におけるある日の 1 時間ごとの交通量が含まれています。

最初に、範囲が 5 時間の移動平均フィルターを使用して、すべてのデータを (線形インデックスによって) 一度に平滑化します。

c = smooth(count(:));
C1 = reshape(c,24,3);

元のデータと平滑化したデータをプロットします。

subplot(3,1,1)
plot(count,':');
hold on
plot(C1,'-');
title('Smooth C1 (All Data)')

次に、同じフィルターを使用してデータの列を別々に平滑化します。

C2 = zeros(24,3);
for I = 1:3,
    C2(:,I) = smooth(count(:,I));
end

元のデータと平滑化したデータを再度プロットします。

subplot(3,1,2)
plot(count,':');
hold on
plot(C2,'-');
title('Smooth C2 (Each Column)')

2 つの平滑化データセットの差をプロットします。

subplot(3,1,3)
plot(C2 - C1,'o-')
title('Difference C2 - C1')

3 列の平滑化の場合に余分な末端効果があることに注意してください。

例: Loess および Robust Loess によるデータの平滑化

外れ値を含むノイズの多いデータを作成します。

x = 15*rand(150,1); 
y = sin(x) + 0.5*(rand(size(x))-0.5);
y(ceil(length(x)*rand(2,1))) = 3;

loess および rloess の手法を使用し、幅 10% でデータを平滑化します。

yy1 = smooth(x,y,0.1,'loess');
yy2 = smooth(x,y,0.1,'rloess');

元のデータと平滑化したデータをプロットします。

[xx,ind] = sort(x);
subplot(2,1,1)
plot(xx,y(ind),'b.',xx,yy1(ind),'r-')
set(gca,'YLim',[-1.5 3.5])
legend('Original Data','Smoothed Data Using ''loess''',...
       'Location','NW')
subplot(2,1,2)
plot(xx,y(ind),'b.',xx,yy2(ind),'r-')
set(gca,'YLim',[-1.5 3.5])
legend('Original Data','Smoothed Data Using ''rloess''',...
       'Location','NW')