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

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

Чтобы проиллюстрировать, как CUSUM работает, сначала исследуйте общие случаи, о которых сообщают, 2 014 западноафриканских вспышек Эболы.

Источник: центры по контролю и профилактике заболеваний

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 diem) in Liberia since March 25, 2014');
ylabel('Change in number of reported cases per diem');
xlabel('Number of days since outbreak began');

Если вы увеличиваете масштаб первой сотни дней данных, вы видите, что, в то время как был начальный приток случаев, многие из них были исключены после дня 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)

Нахождение существенного изменения в отклонении

Другой метод обнаружения резких изменений в статистике посредством обнаружения точек изменения, которое делит сигнал в смежные сегменты где статистическая величина (e.g. среднее значение, отклонение, наклон, и т.д.), является постоянным в каждом сегменте.

Следующий пример анализирует ежегодный минимальный уровень воды реки Нил в течение лет 622 - 1 281 AD, измеренный в приборе Рода под Каиром.

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

Конструкция начала на более новом более точном измерительном приборе приблизительно 715 AD. Не много известно перед этим временем, но на дальнейшем исследовании, вы видите, что существует значительно меньше изменчивости после приблизительно 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(datenum(xp),yp,[.5 .5 .5],'facealpha',0.1);

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

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

Наблюдение изменений многоступенчатого события, разделяемого между сигналами

Чтобы видеть улучшение, увеличьте число changepoints к 20 и наблюдайте изменения в близости переключения передач в демонстрационном номере 19000

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

Заметьте, что скорость вращения двигателя начала уменьшаться в демонстрационных 19035 и взяла 510 выборок, прежде чем она обосновалась в демонстрационных 19550. Поскольку интервал выборки составляет 1 мс, это - ~0.51 задержки с и является типичным количеством времени после изменений скорости.

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

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

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

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

xlim([19000 19050])

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