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

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

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

rng default

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

Сгенерируйте окно наблюдения вокруг каждого элемента x. Примите k = 2 соседа по обе стороны выборки. Скользящее окно, которое результаты, имеют длину $$2\times 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

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

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

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

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

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

См. также

hampel