Аккуратная фильтрация

Пришлось проштудировать немало литературы, прежде чем написать статью о фильтрации выбросов. Вся беда в том, что подавляющее большинство методов фильтрации ориентировано на нормальный закон распределения случайной величины. Наша же случайная величина, которую мы рассматриваем, как ВРТ, распределена не по нормальному закону. Я показал это в нескольких предыдущих статьях.

И что же делать?

Я нашел один непараметрический метод, который не зависит от распределения случайной величины. Этот метод описал в 1959 году J.Walsh в статье «Large sample nonparametric rejection of outlying observations». Этот метод позволяет выявлять выбросы в выборках большого объёма. Саму статью я не нашел, однако разобраться с методом смог по тем материалам, что имеются в Интернет. Когда я написал скрипт по выявлению выбросов в ВРТ и их замене на среднее значение Q, то был разочарован. Нарушается симметрия эмпирической плотности распределения вероятностей. Похоже, что становится непостоянной дисперсия, данные перестают быть стационарными. По этой причине я отказался от этого метода.

Если кому-то интересно, то вот функция выявления выбросов на языке R по упомянутому методу:

#Функция фильтрации. Здесь V – это значения ВРТ, 
#a - уровень значимости гипотезы, k-#предполагаемое количество выбросов,
#d- значение на которое заменяем выброс
 fltV2<-function(V,a,k,d)
 {
 lenV<-length(V)
 n<-length(V)
 m<-trunc(0.5*(1+1/a)^2)
 #Требование к размеру выборки
 if(n<m) {return(V)}
 repeat
 {
 Vi<-getioV(V,a,k)
 if (is.null(Vi)) break;
 V[Vi]<-d
 }
 return(V)
 }
 #-----------------------------------------------------
 #Вспомогательная функция определения позиций выбросов
 getioV<-function(V,a,k)
 {
 n<-length(V)
 m<-trunc(0.5*(1+1/a)^2)
 #Требование к размеру выборки
 if(n<m) {return(NULL)}
 #Требование к количеству выбросов для тестирования
 c<-round(sqrt(2*n))
 if((k+c)>n) {return(NULL)}
 #Алгоритм
 V<-sort.int(V,index.return = TRUE)
 x<-V$x
 b<-sqrt(1/a)
 d<-(1+b*sqrt((c-b^2)/(c-1)))/(c-b^2-1)
 Vi<-NULL
 for(i in 1:k)
 {
 r<-i+c
 if ( (x[n+1-i]-(1+d)*x[n-i]+d*x[n+1-r])>0 ) {Vi<-c(Vi,(V$ix)[n+1-i])}
 if ( (x[i]-(1+d)*x[i+1]+d*x[r])<0 ) {Vi<-c(Vi,(V$ix)[i]) }
 }
 #Возвращает позиции выбросов
 return(Vi)
 }

Подробности этого метода ищите в Интернет на английском языке. Мне же метод не понравился по указанным выше причинам, и я сделал свой алгоритм. По сути этот алгоритм есть более «аккуратная» версия алгоритма, описанного мной здесь. Вот, что он из себя представляет.

  1. Проводится сортировка модулей значений ВРТ по возрастанию. Назову такой ВР Q={q0,…,qn-1}.
  2. Исключаются из Q все повторяющиеся значения.
  3. Строится эмпирическая функция распределения значений Q.
  4. Проводится кусочно-линейное сглаживание значений полученной функции распределения.
  5. Для заданного уровня значимости α (например, 0.05) по сглаженной функции распределения определяется граничное значение qe. Значениями верхней h и нижней l границ для выбросов будут значения h=qe и l=-qe.
  6. Выбросами считаются все значения ВРТ, которые удовлетворяют условию: |si|≥qe. В допустимых границах будут находиться все значения ВРТ с вероятностью 1- α.
  7. Все найденные выбросы внутри ВРТ заменить на среднее значение Q с соответствующим выбросу знаком (+/-).

Вот такой более аккуратный алгоритм. Он позволяет сохранить симметрию эмпирической функции плотности распределения вероятностей и стационарность самого ВРТ.

Вот необходимые функции на языке R:

#-----------------------------------------------
 #Построение эмпирической функции распределения
 getEDFx<-function(x)
 {
 lenx<-length(x)
 #Ограничение применимости
 if (lenx==0) return (NULL)
 #Строим функцию распределения
 x<-abs(x)
 x<-sort(unique(x))
 f<-(1:length(x))/length(x)
 EDFx<-data.frame(f,x)
 names(EDFx)<-c("EDF","x")
 return (EDFx)
 }
 #-----------------------------------------------------
 #Кусочно-линейное сглаживание функции распределения
 #и определение x в зависимости от доверительной вероятности a=1-α
 smV<-function(EDFx,a)
 {
 #Ограничения применимости
 if(!is.data.frame(EDFx)) {return (NA)}
 f<-as.vector(unlist(EDFx["EDF"])); x<-as.vector(unlist(EDFx["x"])) ; n<-length(x)
 if (a<=min(f)) {return (0)}
 if (a>=max(f)) {return (1)}
 #Определение интервала, куда попадает значение a
 i<-1
 while(i<n)
 {
 if((a>=f[i])&&(a<=f[i+1])) {break;}
 i<-i+1
 }
 #Строим линию регрессии
 x<-rbind(x[i],x[i+1])
 f<-rbind(f[i],f[i+1])
 M<-lm(x~f,data.frame(x,f))
 f<-a; f<-data.frame(f)
 return(as.vector(predict(M,f)))
 }

Вот график сглаженной функции распределения модулей значений ВРТ:

EDFsmooth

Получается, как видно по графику, вполне правдоподобно.

И для примера график ВРТ после фильтрации выбросов:

smoothVRT

А вот, что было до фильтрации:

VRT

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

Похожие статьи:

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *