Обнаружение вспышек через совокупные суммы

Существует много практических приложений, в которых вы отслеживаете данные и хотите быть предупреждены как можно быстрее, когда базовый процесс изменился. Очень популярный метод достижения этого - с помощью совокупной суммы (CUSUM) диаграммы управления.

Чтобы проиллюстрировать, как работает CUSUM, сначала исследуйте общее количество зарегистрированных случаев вспышки Эболы в Западной Африке в 2014 году, как это было зарегистрировано Центрами по контролю и профилактике заболеваний.

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

Если вы посмотрите на переднее ребро первой вспышки в Гвинее, то увидите, что первые сто случаев были зарегистрированы примерно 25 марта 2014 года и значительно увеличились после этой даты. Интересно отметить, что, хотя в марте в Либерии было также несколько предполагаемых случаев, число случаев оставалось относительно контролируемым примерно тридцать дней спустя.

Чтобы получить представление о приеме новых пациентов, постройте график относительных повседневных изменений общего количества случаев, начиная с начала 25 марта 2015 года.

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

Если вы увеличите масштаб первых ста дней данных, то увидите, что, хотя был первоначальный приток случаев, многие из них были исключены после 30-го дня, где частота изменений временно опустилась ниже нуля. Вы также видите значительный тренд к росту между 95 и 100 днями, где был достигнут показатель в семь новых случаев в день.

xlim([1 101])

Выполнение теста CUSUM на входных данных может быть быстрым способом определить, когда происходит вспышка. CUSUM отслеживает две совокупные суммы: верхнюю сумму, которая обнаруживает, когда локальное среднее смещается вверх, и более низкую сумму, которая обнаруживает, когда среднее смещается вниз. Метод интегрирования предоставляет CUSUM возможность игнорировать большой (переходный) всплеск входящей скорости, но все еще имеет чувствительность к более устойчивым небольшим изменениям скорости.

Вызов CUSUM с аргументами по умолчанию будет проверять данные первых двадцати пяти выборок и предупреждение, когда он сталкивается со сдвигом в среднем более чем на пять стандартных отклонений от исходных данных.

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

Обратите внимание, что CUSUM поймал ложные зарегистрированные случаи на 30 день (на 33 день) и подобрал начальное начало вспышки, начиная с 80 дня (на 90 день). Если вы тщательно сравните эти результаты с предыдущим графиком, можно увидеть, что CUSUM смог проигнорировать ложное повышение в день 29, но все еще запускает предупреждение за пять дней до большого восходящего тренда, начиная с дня 95.

Если вы корректируете 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 значительно меньше. Чтобы найти период времени, когда новое устройство заработало, можно искать лучшее изменение среднеквадратичного уровня воды после выполнения поэлементной дифференциации, чтобы удалить любые медленно изменяющиеся тренды.

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

Обнаружение нескольких изменений в входном сигнале

Следующий пример касается 45-секундной симуляции блока 4-ступенчатой передачи CR-CR, дискретизированной с интервалами 1 мс. Данные моделирования оборотов двигателя автомобиля и крутящего момента показаны ниже.

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 выборки. Поскольку интервал дискретизации составляет 1 мс, это ~ задержка 0,51 с и является типичным количеством времени после переключения передач.

Теперь посмотрите на точки изменения крутящего момента двигателя в той же области:

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

Заметьте, что крутящий момент двигателя был полностью подан на ось через выборку 19605, 55 миллисекунд после завершения осаждения скорости вращения двигателя. Это время связано с задержкой между штрихом впуска двигателя и производством крутящего момента.

Чтобы найти, когда муфта зацепилась, вы можете масштабировать дальше в сигнал.

xlim([19000 19050])

Муфта была нажата на выборку 19011 и взяла около 30 выборок (миллисекунд), чтобы полностью отсоединиться.