Устраните выбросы Используя идентификатор Хампеля

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

Сгенерируйте случайный сигнал, x, содержа 24 выборки. Сбросьте генератор случайных чисел для восстанавливаемых результатов.

rng default

lx = 24;
x = randn(1,lx);

Сгенерируйте окно наблюдения вокруг каждого элемента x. Возьмите k = 2 соседа в любой стороне выборки. Движущееся окно, которого результаты имеют длину 2×2+1=5 выборки.

k = 2;

iLo = (1:lx)-k;
iHi = (1:lx)+k;

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

iLo(iLo<1) = 1;
iHi(iHi>lx) = lx;

Запишите медиану каждого окружающего окна. Найдите медиану абсолютного отклонения каждого элемента относительно медианы окна.

for j = 1:lx
    w = x(iLo(j):iHi(j));
    medj = median(w);
    mmed(j) = medj;
    mmad(j) = median(abs(w-medj));
end

Масштабируйте среднее абсолютное отклонение с

12erf-1(1/2)1.4826

получить оценку стандартного отклонения нормального распределения.

sd = mmad/(erfinv(1/2)*sqrt(2));

Найдите выборки, которые отличаются от медианы больше, чем nd = 2 стандартных отклонения. Замените каждые из тех выбросов значением медианы его окружающего окна. Это - сущность алгоритма Hampel.

nd = 2;
ki = abs(x-mmed) > nd*sd;

yu = x;
yu(ki) = mmed(ki);

Используйте hampel функция, чтобы вычислить отфильтрованный сигнал и аннотировать выбросы. Наложите отфильтрованные значения, вычисленные в этом примере.

hampel(x,k,nd)

hold on
plot(yu,'o','HandleVisibility','off')
hold off

Figure contains an axes object. The axes object contains 3 objects of type line. These objects represent original signal, filtered signal, outliers.

Смотрите также