Main Content

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

이 예제에서는 누적합과 변화 지점 감지를 통해 신호의 변화 또는 돌연한 출현을 확인하는 방법을 보여줍니다.

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

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

Figure contains an axes object. The axes object with title Total Suspected, Probable, and Confirmed Cases of Ebola Virus Disease contains 3 objects of type line. These objects represent Guinea, Liberia, Sierra Leone.

기니에서 첫 번째 발병이 시작하는 구간을 살펴보면 최초 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');

Figure contains an axes object. The axes object with title Rate of New Cases (per Day) in Liberia since March 25, 2014, xlabel Number of Days from Start of Outbreak, ylabel Change in Number of Reported Cases per Day contains an object of type line.

데이터의 처음 100일을 확대해 보면, 초기에 발병 건수가 확 늘어나긴 했으나, 30일 차 후에는 그중 다수가 배제되어 변동률이 일시적으로 0 아래로 떨어진 것을 볼 수 있습니다. 일일 신규 건수가 7건에 도달한 95일 차와 100일 차 사이에 급격한 상승 추세가 있는 것도 볼 수 있습니다.

xlim([1 101])

Figure contains an axes object. The axes object with title Rate of New Cases (per Day) in Liberia since March 25, 2014, xlabel Number of Days from Start of Outbreak, ylabel Change in Number of Reported Cases per Day contains an object of type line.

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

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

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

Figure contains an axes object. The axes object with title CUSUM Control Chart mu indexOf target baseline blank = blank 1 . 080000 blank sigma indexOf target baseline blank = blank 1 . 629734, xlabel Samples, ylabel Standard Errors contains 6 objects of type line. One or more of the lines displays its values using only markers These objects represent 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)

Figure contains an axes object. The axes object with title CUSUM Control Chart mu indexOf target baseline blank = blank 0 . 000000 blank sigma indexOf target baseline blank = blank 3 . 000000, xlabel Samples, ylabel Standard Errors contains 5 objects of type line. One or more of the lines displays its values using only markers

분산의 급격한 변화 찾기

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

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

load nilometer

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

Figure contains an axes object. The axes object with title Yearly Minimum level of Nile River, xlabel Year, ylabel Level (m) contains an object of type line.

서기 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)

Figure contains an axes object. The axes object with title Yearly Minimum level of Nile River, xlabel Year, ylabel Level (m) contains 2 objects of type line, patch.

샘플별 미분은 추세를 제거하는 간단한 방법이긴 하나 대규모 척도에 대해 분산을 조사하는 더욱 정밀한 방법도 존재합니다. 이 데이터셋을 사용하여 웨이블릿을 통해 변화 지점 감지를 수행하는 예제는 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')

Figure contains 2 axes objects. Axes object 1 with title Torque, xlabel Samples, ylabel Torque (N m) contains an object of type line. Axes object 2 with title Engine Speed, xlabel Samples, ylabel Speed (RPM) contains an object of type line.

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

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

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

Figure contains an axes object. The axes object with title Number of changepoints = 4 Total residual error = 342535085.7709, xlabel Samples, ylabel Engine speed (RPM) contains 3 objects of type line.

여기서 (다섯 개의 선형 세그먼트 간에) 네 번의 변경을 볼 수 있으며, 각각 10,000, 20,000, 30,000, 40,000번째 샘플 지점에서 발생했음을 알 수 있습니다. 파형의 유휴 부분을 확대합니다.

xlim([18000 22000])

Figure contains an axes object. The axes object with title Number of changepoints = 4 Total residual error = 342535085.7709, xlabel Samples, ylabel Engine speed (RPM) contains 3 objects of type line.

직선 피팅이 입력 파형을 근접하게 추적합니다. 하지만 피팅을 더 개선할 수 있습니다.

신호 간 공유되는 다단계 이벤트의 변화 관찰하기

이를 개선하려면 변화 지점의 개수를 20으로 늘리고 19000번째 샘플에서 발생한 변속 인근에서 변경 사항을 관찰하십시오.

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

Figure contains an axes object. The axes object with title Number of changepoints = 20 Total residual error = 983541.4074 contains 3 objects of type line.

19035번째 샘플에서 엔진 속도가 감소하기 시작했고 510개의 샘플 후인 19550번째 샘플 지점에서 안정화되기 시작했음을 관찰할 수 있습니다. 샘플링 간격이 1ms이므로 이는 약 0.51초의 지연에 해당하는데, 이는 변속 후에 일반적으로 소요되는 시간입니다.

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

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

Figure contains an axes object. The axes object with title Number of changepoints = 20 Total residual error = 435882.0922 contains 3 objects of type line.

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

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

xlim([19000 19050])

Figure contains an axes object. The axes object with title Number of changepoints = 20 Total residual error = 435882.0922 contains 3 objects of type line.

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

참고 항목

| |