Обнаружение вспышек с помощью накопительных сумм

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

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

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

Следующий пример касается 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 образцов (миллисекунд), чтобы полностью разцепиться.