信号の急激な増加と大幅な変化の検出

累積和による急激な増加の検出

データを監視し、基本的な過程に変化があったときにできる限り迅速に警告を受けたい、というような実際的応用事例は数多くあります。これを行う非常に一般的な手法として、累積和 (CUSUM) 管理図の使用があります。

CUSUM の機能を説明するため、はじめに、Centers for Disease Control and Prevention によって記録された、2014 年に西アフリカで発生したエボラ大流行の全報告例を検証します。

load WestAfricanEbolaOutbreak2014
plot(WHOreportdate, [TotalCasesGuinea TotalCasesLiberia TotalCasesSierraLeone],'.-')
legend('Guinea','Liberia','Sierra Leone')
title('Total Suspected, Probable, and Confirmed Cases of Ebola Virus Disease')

ギニアにおける大流行の先端部を見てみると、2014 年 3 月 25 日頃に最初の 100 例が報告され、その日以後、大幅に件数が増加したことがわかります。興味深いのは、3 月にはリベリアにも少数の疑い例があるものの、発生件数は約 30 日後まで比較的抑制されていたという点です。

新しい患者の発生率を感覚的に把握するために、2014 年 3 月 25 日を開始時点として、合計症例数の日々の相対的変動をプロットします。

daysSinceOutbreak = datetime(2014, 3, 24+(0:400));
cases = interp1(WHOreportdate, TotalCasesLiberia, daysSinceOutbreak);
dayOverDayCases = diff(cases);

plot(dayOverDayCases)
title('Rate of New Cases (per Day) in Liberia since March 25, 2014');
ylabel('Change in Number of Reported Cases per Day');
xlabel('Number of Days from Start of Outbreak');

最初の 100 日間のデータを拡大すると、初期には症例数が急増しているものの、それらの多くは 30 日目の後に排除されており、変化率は一時的に 0 未満に低下しています。また、95 日目と 100 日目の間には著しい上昇トレンドが示されています。ここでは、1 日あたりの新規発症件数の割合が 7 件に達しています。

xlim([1 101])

入力データに対して CUSUM テストを実行すると、急激な増加発生のタイミングを簡単に特定できる場合があります。CUSUM は、局所的な平均が上向きに変化するタイミングを検出する上位の和と、平均が下向きに変化するタイミングを検出する下位の和の 2 種類の累積和を追跡します。CUSUM には、積分手法によって入力レートにおける大きな (過渡的な) スパイクを無視する機能がありますが、より安定した小さいレートの変化に対する感度は維持されます。

既定の引数で CUSUM を呼び出すと、最初の 25 サンプルのデータを検査し、初期データ内から 5 標準偏差を超える平均の変化を検知すると警告を発します。

cusum(dayOverDayCases(1:101))
legend('Upper sum','Lower sum')

CUSUM は、30 日目 (33 日目) の誤った報告件数を捕捉し、80 日目 (90 日目) に始まる急激な増加の最初の発生を検出していることに注意してください。これらの結果を前のプロットと注意深く比較すると、CUSUM は、29 日目の誤った上昇を無視する一方で、95 日目に始まる大きな上昇トレンドの 5 日前に警告を発することが可能であったことがわかります。

ターゲット平均値 0 件/日とターゲット ± 3 件/日を与えて CUSUM を調整すると、30 日目の誤検出警報を無視し、かつ 92 日目の急激な増加を検出することができます。

climit = 5;
mshift = 1;
tmean = 0;
tdev = 3;
cusum(dayOverDayCases(1:100),climit,mshift,tmean,tdev)

分散の著しい変化の検出

統計量の急激な変化を検出するもう 1 つの方法は、変化点の検出を使用するものです。この方法では、各セグメント内において統計量 (平均、分散、勾配など) が一定な隣接したセグメントに信号を分割します。

次の例では、カイロに近いローダ島の水位計で測定された西暦 622 年から 1281 年のナイル川の年間最低水位を解析します。

load nilometer

years = 622:1284;
plot(years,nileriverminima)
title('Yearly Minimum level of Nile River')
xlabel('Year')
ylabel('Level (m)')

西暦 715 年頃に、より新しく正確な測定機器の製作が始まりました。これ以前についてあまり多くのことは知られていませんが、さらに詳しく調べてみると、722 年頃以降の変動性が大幅に低下していることがわかります。新しい機器が使用されるようになった期間を特定するため、要素単位の微分を行って変化の緩やかなトレンドを除去した後で、平方根平均二乗水位における最良の変化を探索します。

i = findchangepts(diff(nileriverminima),'Statistic','rms');

ax = gca;
xp = [years(i) ax.XLim([2 2]) years(i)];
yp = ax.YLim([1 1 2 2]);
patch(xp,yp,[.5 .5 .5],'FaceAlpha',0.1)

サンプル単位の微分はトレンドを除去する簡単な方法ですが、この他、より大きなスケールで分散を調査する高度な方法もあります。このデータセットを使用した、ウェーブレットによる変化点検出の実行方法の例については、Wavelet Changepoint Detection (Wavelet Toolbox)を参照してください。

入力信号における複数の変化の検出

次の例は、1 ms 間隔でサンプリングされた CR-CR 4 速トランスミッション ブロックの 45 秒間のシミュレーションを扱っています。自動車のエンジン速度とトルクのシミュレーション データは次のとおりです。

load simcarsig

subplot(2,1,2)
plot(carTorqueNM)
xlabel('Samples')
ylabel('Torque (N m)')
title('Torque')

subplot(2,1,1);
plot(carEngineRPM)
xlabel('Samples')
ylabel('Speed (RPM)')
title('Engine Speed')

ここでは、車は、加速してギアを 3 度変え、ニュートラルに切り替えてからブレーキを使用しています。

エンジン速度は一連の線形セグメントとして自然にモデル化できるため、findchangepts を使用して車がギアを変えた時点のサンプルを探すことができます。

figure
findchangepts(carEngineRPM,'Statistic','linear','MaxNumChanges',4)
xlabel('Samples')
ylabel('Engine speed (RPM)')

ここでは、4 度の変化 (5 つの線形セグメントの間) が確認でき、それらが 10,000、20,000、30,000、40,000 のサンプル マーク付近で発生したことがわかります。波形のアイドル部分を拡大します。

xlim([18000 22000])

直線近似は入力波形を厳密に辿ります。ただし、この近似はさらなる改善が可能です。

信号間で共有される多段イベントの変化の観察

改善を確認するには、変化点の数を 20 に増やし、サンプル番号 19000 のギア チェンジの近傍内における変化を観察します。

findchangepts(carEngineRPM,'Statistic','linear','MaxNumChanges',20)
xlim([18000 22000])

エンジン速度がサンプル 19035 で低下し始め、510 サンプルを経過した後、サンプル 19550 で安定したことを確認します。サンプリング間隔は 1 ms のため、これは 0.51 秒程度の遅延であり、ギア チェンジ後の標準的な時間量です。

次に、同じ領域内のエンジン トルクの変化点を調べます。

findchangepts(carTorqueNM,'Statistic','Linear','MaxNumChanges',20)
xlim([19000 20000])

エンジン トルクは、エンジン速度が安定してから 55 ミリ秒後のサンプル 19605 で完全に車軸に伝達されたことを確認します。この時間は、エンジンの吸気行程とトルクの生成の間の遅延に関連しています。

クラッチを接続したタイミングを特定するために、さらに信号を拡大できます。

xlim([19000 19050])

クラッチは、サンプル 19011 で踏み込まれ、約 30 サンプル (ミリ秒) かけて完全に開放されました。