신호의 돌연한 출현 및 급격한 변화 감지하기

누적합을 통해 발병 감지하기

다수의 실제 응용 사례에서 데이터를 모니터링할 때 기본 과정이 변경된 경우 가능한 한 빠르게 알림을 받아야 하는 경우가 있습니다. 이를 위해 널리 사용되는 기법 중 하나는 누적합(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')

기니에서 첫 번째 발병이 시작하는 구간을 살펴보면 최초 100건의 발병 사례가 2014년 3월 25일경에 보고되었고 이날 이후 발병 건이 급격히 증가한 것을 볼 수 있습니다. 흥미로운 점은, 라이베리아에서도 3월에 몇 건의 의심되는 사례가 보고되긴 했으나 그로부터 약 30일 후까지는 발병 건수가 상대적으로 통제된 상태에 머물렀다는 것입니다.

신규 환자의 발생 속도를 가늠해 보려면 발병이 개시된 2015년 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 아래로 떨어진 것을 볼 수 있습니다. 일일 신규 건수가 7건에 도달한 95일 차와 100일 차 사이에 급격한 상승 추세가 있는 것도 볼 수 있습니다.

xlim([1 101])

입력 데이터에 대해 CUSUM 테스트를 수행하면 발병이 발생한 시점을 빠르게 확인할 수 있습니다. CUSUM은 국소 평균이 상승하는 시점을 감지하는 상합(upper sum)과 국소 평균이 하강하는 시점을 감지하는 하합(lower sum)의 두 가지 누적합을 추적합니다. 이러한 적분 기법 덕분에 CUSUM은 신규 환자 발생 속도의 큰(과도) 스파이크는 무시하면서도 비교적 평탄한 작은 속도 변화에 대한 민감도를 가질 수 있습니다.

디폴트 인수와 함께 CUSUM을 호출하면 처음 25개 샘플 데이터에 대한 조사가 이루어지고, 초기 데이터 내에서 평균이 5 표준편차를 초과하여 변경되면 경고가 발생됩니다.

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

여기서 CUSUM이 30일 차에 보고된 거짓 사례를 33일 차에 발견했고 80일차에 시작된 발병의 초기 개시를 90일 차에 알아차렸음을 알 수 있습니다. 이 결과를 직전 플롯과 신중하게 비교해 보면 CUSUM이 29일 차의 스퓨리어스 증가는 무시하면서도 큰 상향 추세가 시작되는 95일 차로부터 5일 전에 경고를 작동시켰음을 볼 수 있습니다.

일일 ±3건의 사례라는 목표치와 일일 0건의 사례라는 목표 평균을 갖도록 CUSUM을 조정하면 30일 차의 거짓 경고를 무시하고 92일 차의 발병을 발견할 수 있습니다.

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

분산의 급격한 변화 찾기

통계량의 급격한 변화를 감지하는 또 다른 방법은 변화 지점 감지를 통한 것입니다. 변화 지점 감지는 신호를 각 세그먼트 내의 통계량(예: 평균, 분산, 기울기 등)이 일정한 인접한 세그먼트들로 분할합니다.

다음 예제에서는 카이로 인근의 로다 측정소에서 측정된 서기 622년~1281년의 연간 나일강 최저 수위를 분석합니다.

load nilometer

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

서기 715년 즈음부터 보다 정확한 측정 장치를 사용한 건설이 시작되었습니다. 이 시점 전까지는 알려진 데이터가 많지 않지만, 자세히 살펴보면 722년 즈음부터 변동성이 상당히 줄어드는 것을 볼 수 있습니다. 새로운 장치가 사용되기 시작한 시기를 찾기 위해, 요소별 미분을 수행하여 느리게 변화하는 추세를 제거한 후 RMS(제곱평균제곱근) 수위의 최적의 변화를 탐색할 수 있습니다.

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) 항목을 참조하십시오.

입력 신호에서 여러 개의 변경 감지하기

다음 예제에서는 1ms 간격으로 샘플링된 CR-CR 방식 4단 변속기 블록의 45초 시뮬레이션을 사용합니다. 아래에 자동차 엔진 RPM과 토크의 시뮬레이션 데이터가 나와 있습니다.

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

여기서 자동차는 가속을 하고 기어를 세 번 바꾼 후 중립으로 전환한 다음 브레이크를 적용합니다.

엔진 속도는 자연스레 선형 세그먼트로 구성된 계열로 모델링될 수 있으므로, findchangepts를 사용하여 자동차의 변속이 이루어진 샘플을 찾을 수 있습니다.

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

여기서 (다섯 개의 선형 세그먼트 간에) 네 번의 변경을 볼 수 있으며, 각각 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번째 샘플 지점에서 안정화되기 시작했음을 관찰할 수 있습니다. 샘플링 간격이 1ms이므로 이는 약 0.51초의 지연에 해당하는데, 이는 변속 후에 일반적으로 소요되는 시간입니다.

이번에는 같은 영역 내에서 엔진 토크의 변화 지점을 살펴보겠습니다.

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

엔진 토크가 19605번째 샘플에서, 즉 엔진 속도의 안정화가 이루어지고 나서 55밀리초 후에 차축으로 완전히 전달되었음을 관찰할 수 있습니다. 이 시간은 엔진의 흡입 행정과 토크 생성 사이의 지연과 관련이 있습니다.

클러치가 결합되는 시점을 확인하려면 신호를 더 확대해 볼 수 있습니다.

xlim([19000 19050])

클러치는 19011번째 샘플에서 눌려서 약 30개 샘플(밀리초) 후에 완전히 분리되었습니다.