Удалить значения за пределами кривой Лёсса

Я хочу удалить выбросы, прежде чем применять модель. Я использую кривую Лесса, чтобы разграничить линию тренда и установить пределы выбросов. Я хотел бы удалить строки, выходящие за установленные пределы. Помимо того, что это делается с помощью пользовательской функции, которая берет каждую точку по одной и проверяет локальный уклон лёсса и т. д., есть ли более простой способ?

Линия тренда лесса с ограничениями (1.2)

# Code generating image above
scatter.smooth( idam$T_d, idam$T_x10d)
loessline <- loess.smooth( idam$T_d, idam$T_x10d)
lines(loessline$x, loessline$y, lwd=3)
lines(loessline$x, loessline$y*1.2, lwd=3, col='red')
lines(loessline$x, loessline$y/1.2, lwd=3, col='red')

person Cyrille    schedule 07.11.2014    source источник
comment
Будь осторожен. Если физика не дает веских причин для выбрасывания данных, то, возможно, вы поступаете плохо, выбрасывая хорошие, но менее удобные данные.   -  person EngrStudent    schedule 05.10.2017


Ответы (3)


Обнаружение выбросов может быть выполнено с помощью пакета DBSCAN R, известный алгоритм, используемый для идентификации кластера (подробности см. на WIKIPEDIA).

Эта функция имеет три важных входа:

  • x: ваши данные (только числовое значение)
  • eps: целевое максимальное расстояние
  • minPts: количество минимальных точек, чтобы рассматривать их как кластер

Оценить eps можно с помощью функций knndist(...) и knndistplot(...):

  • knndistplot будет отображать значения eps в вашем наборе данных для заданного k (т.е. minPts) ==> вы можете визуально выбрать эффективное значение eps (обычно в части кривой колена)
  • knndist оценит значения eps и вернет их в матрицу из. ввод k будет генерировать оценки 1: 1: k, и вы можете использовать результат для программного определения точных значений eps и k

Далее вам нужно просто использовать dbscan(yourdata, eps, k) для получения объекта dbscan со следующими компонентами:

  • eps: количество eps, используемое для расчетов
  • minPts: минимальное количество точек для идентификации кластера
  • кластер: целочисленный вектор, определяющий точки, которые принадлежат кластеру (=1) или нет (=0). Последние соответствуют выбросам, которые вы хотите устранить.

Обратите внимание на следующие ограничения для dbscan:

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

Зная об этих ограничениях, пакет dbscan предлагает два альтернативных метода: LOF и OPTICS (расширение DBSCAN).

Редактировать 25 января 2016 г.

Следуя ответу @rawr, я привожу пример, основанный на наборе данных mtcars, чтобы показать, как использовать dbscan для выявления выбросов. Обратите внимание, что мой пример будет использовать отличный пакет data.table вместо классического data.frame.

Во-первых, я начинаю копировать подход rawr для иллюстрации использования data.table.

require(data.table)
require(ggplot2)
require(dbscan)
data(mtcars)
dt_mtcars <- as.data.table(mtcars)

# based on rawr's approach
plot(wt~mpg, data=dt_mtcars)
lo <- loess.smooth(dt_mtcars[,mpg], dt_mtcars[,wt])
lines(lo$x,lo$y, lwd=3)
lines(lo$x,lo$y * 1.2, lwd=3 , col=2 )
lines(lo$x,lo$y / 1.2, lwd=3 , col=2 )

введите здесь описание изображения

Таким образом, мы можем оценить, что получаем одинаковые результаты независимо от базовой поддержки.

Во-вторых, следующий код иллюстрирует подход DBSCAN, который начинается с определения eps и k, необходимого количества точек для идентификации кластера:

res_knn = kNNdist( dt_mtcars[, .(wt, mpg)] , k = 10)
dim_knn = dim(res_knn)
x_knn =  seq(1, dim_knn[1])
ggplot() + 
   geom_line( aes( x = x_knn , y = sort(res_knn[, 1])  , col = 1 ) ) + 
   geom_line( aes( x = x_knn , y = sort(res_knn[, 2])  , col = 2 ) ) + 
   geom_line( aes( x = x_knn , y = sort(res_knn[, 3])  , col = 3 ) ) + 
   geom_line( aes( x = x_knn , y = sort(res_knn[, 4])  , col = 4 ) ) + 
   geom_line( aes( x = x_knn , y = sort(res_knn[, 5])  , col = 5 ) ) + 
   geom_line( aes( x = x_knn , y = sort(res_knn[, 6])  , col = 6 ) ) + 
   geom_line( aes( x = x_knn , y = sort(res_knn[, 7])  , col = 7 ) ) + 
   geom_line( aes( x = x_knn , y = sort(res_knn[, 8])  , col = 8 ) ) + 
   geom_line( aes( x = x_knn , y = sort(res_knn[, 9])  , col = 9 ) ) + 
   geom_line( aes( x = x_knn , y = sort(res_knn[, 10]) , col = 10 ) )  +
   xlab('sorted results') + 
   ylab('kNN distance')

Результаты представлены на следующем графике:

введите здесь описание изображения

Он показывает, что рассчитанное расстояние kNN чувствительно к фактору k, однако точное значение eps для разделения выбросов находится в коленной части кривых ==> подходящее eps расположено между 2 и 4. Это визуальная оценка, которая может быть автоматизированы с помощью соответствующих алгоритмов поиска (например, см. эту ссылку). Что касается k, необходимо определить компромисс, зная, что чем ниже k, тем менее строгие результаты.

В следующей части мы параметризуем dbscan с помощью eps = 3 (на основе визуальной оценки) и k = 4 для получения слегка строгих результатов. Мы построим эти результаты с помощью кода rawr:

eps = 3
k = 4
res_dbscan = dbscan( dt_mtcars[, .(wt, mpg)] , eps , k )
plot(wt~mpg, data=dt_mtcars, col = res_dbscan$cluster)
lo <- loess.smooth(dt_mtcars[res_dbscan$cluster>0,mpg], dt_mtcars[res_dbscan$cluster>0,wt])
lines(lo$x,lo$y, lwd=3)
lines(lo$x,lo$y * 1.2, lwd=3 , col=2 )
lines(lo$x,lo$y / 1.2, lwd=3 , col=2 )

введите здесь описание изображения

мы получили эту цифру, где мы можем оценить, что мы получили разные результаты из подхода rawr, где точки, расположенные в mpg = [10,13], считаются выбросами.

Эти результаты можно рассматривать как странные по сравнению с решением rawr, которое работает в предположении наличия двумерных данных (Y ~ X). Однако mtcars представляет собой многомерный набор данных, в котором взаимосвязь между переменными может быть (или нет) линейной... Чтобы оценить эту точку, мы могли бы построить диаграмму рассеивания этого набора данных, отфильтрованную, например, по числовым значениям.

pairs(dt_mtcars[, .(mpg, disp, hp, drat, wt, qsec)])

введите здесь описание изображения

Если мы сосредоточимся только на результате wt ~ mpg, на первый взгляд мы можем подумать, что это антилинейная зависимость. Но с другими отношениями на графике это может быть не так, и поиск выбросов в среде N-Dim немного сложнее. Действительно, одна точка может рассматриваться как выброс при проецировании в конкретном 2D-сравнении... но наоборот, если мы добавим новое измерение сравнения. Действительно, у нас могла бы быть коллинеарность, которую можно было бы идентифицировать и, таким образом, усилить связь кластера или нет.

Друзья мои, я согласен, что это много if и чтобы проиллюстрировать эту ситуацию, мы перейдем к dbscan анализу числовых значений mtcars.

Поэтому я повторю процесс, представленный ранее, и начнем с анализа расстояния kNN:

res_knn = kNNdist( dt_mtcars[, .(mpg, disp, hp, drat, wt, qsec)] , k = 10)
dim_knn = dim(res_knn)
x_knn =  seq(1, dim_knn[1])
ggplot() + 
    geom_line( aes( x = x_knn , y = sort(res_knn[, 1])  , col = 1 ) ) + 
    geom_line( aes( x = x_knn , y = sort(res_knn[, 2])  , col = 2 ) ) + 
    geom_line( aes( x = x_knn , y = sort(res_knn[, 3])  , col = 3 ) ) + 
    geom_line( aes( x = x_knn , y = sort(res_knn[, 4])  , col = 4 ) ) + 
    geom_line( aes( x = x_knn , y = sort(res_knn[, 5])  , col = 5 ) ) + 
    geom_line( aes( x = x_knn , y = sort(res_knn[, 6])  , col = 6 ) ) + 
    geom_line( aes( x = x_knn , y = sort(res_knn[, 7])  , col = 7 ) ) + 
    geom_line( aes( x = x_knn , y = sort(res_knn[, 8])  , col = 8 ) ) + 
    geom_line( aes( x = x_knn , y = sort(res_knn[, 9])  , col = 9 ) ) + 
    geom_line( aes( x = x_knn , y = sort(res_knn[, 10]) , col = 10 ) )  +
    xlab('sorted results') + 
    ylab('kNN distance')

отсортированные расстояния kNN

По сравнению с анализом, проведенным на wt ~ mpg, мы видим, что kNNdist(...) дает гораздо более важное расстояние kNN (например, до 200 с k = 10). Однако у нас все еще есть часть колена, которая помогает нам оценить подходящее значение eps.

В следующей части мы будем использовать eps = 75 и k = 5 и

# optimal eps value is between 40 (k=1) and 130 (k=10)
eps = 75
k = 5
res_dbscan = dbscan( dt_mtcars[, .(mpg, disp, hp, drat, wt, qsec)] , eps , k )
pairs(dt_mtcars[, .(mpg, disp, hp, drat, wt, qsec)] , col = res_dbscan$cluster+2L)

введите здесь описание изображения

Таким образом, диаграмма рассеяния этого анализа показывает, что выявление выбросов может быть затруднено в среде N-Dim из-за сложных взаимосвязей между переменными. Но обратите внимание, что в большинстве случаев выбросы расположены в угловой части 2D-проекции, что усиливает результаты, полученные с помощью wt ~ mpg.

person Bruno Sarrant    schedule 19.01.2016
comment
Как можно применить алгоритм кластеризации к данным временных рядов? - person user2398029; 19.01.2016
comment
@ user2398029: данные временных рядов также можно рассматривать как числовые значения. Например, Excel или Matlab определяют точку отсчета времени и используют формат даты как псевдоним формата ГГГГ-ММ-ДД. Таким образом, вы можете переключаться со строкового значения на числовое. Еще одним способом интеграции данных временных рядов является наличие сезонности за период, поэтому временные данные можно рассматривать как фактор суммирования. Вы также можете преобразовать данные о времени как коэффициент, различая год, месяц, день и т. д. - person Bruno Sarrant; 20.01.2016
comment
Конечно... Я спрашиваю, как применить алгоритм кластеризации к данным, которые не образуют кластеры. - person user2398029; 20.01.2016
comment
Спасибо за разъяснение, и вот в чем вопрос: как вы можете оценить, что ваши данные не образуют кластер, даже если у вас есть двумерный случай (значение и время)? Например, у вас могут быть сезонные эффекты, которые объединяют результаты в определенные моменты времени... Конечно, если у вас строгий линейный эффект, использование DBSCAN может быть бесполезным в 2D-случае... но если у вас есть многомерные данные, некоторые корреляции может возникнуть и быть скрытым при отражении ваших данных в представлении с меньшими размерами. Я работаю над редактированием своего ответа на основе набора данных mtcars, чтобы выделить этот момент. - person Bruno Sarrant; 21.01.2016

Вы можете использовать approxfun

Вот пример с "выбросами"

plot(wt ~ mpg, data = mtcars)
lo <- loess.smooth(mtcars$mpg, mtcars$wt)
lines(lo$x, lo$y, lwd = 3)
lines(lo$x, lo$y * 1.2, lwd = 3, col = 2)
lines(lo$x, lo$y / 1.2, lwd = 3, col = 2)

введите здесь описание изображения

approxfun возвращает функцию, использующую наблюдаемые значения x, которые мы можем использовать для интерполяции набора новых значений y.

Затем вы можете установить порог для того, чтобы назвать точку выбросом; здесь я использую 1.2 * y, как и в исходном вопросе, для определения крайних наблюдений.

f1 <- approxfun(lo$x, lo$y * 1.2)
(wh1 <- which(mtcars$wt > f1(mtcars$mpg)))
# [1]  8 17 18

f2 <- approxfun(lo$x, lo$y / 1.2)
(wh2 <- which(mtcars$wt < f2(mtcars$mpg)))
# [1] 28

## identify points to exclude
mt <- mtcars[c(wh1, wh2), ]
points(mt$mpg, mt$wt, pch = 4, col = 2, cex = 2)

введите здесь описание изображения

## plot without points
plot(wt ~ mpg, data = mt2 <- mtcars[-c(wh1, wh2), ])
lo <- loess.smooth(mt2$mpg, mt2$wt)
lines(lo$x, lo$y, lwd = 3)
lines(lo$x, lo$y * 1.2, lwd = 3, col = 2)
lines(lo$x, lo$y / 1.2, lwd = 3, col = 2)

введите здесь описание изображения

Поскольку здесь есть пара шагов, вы можете упаковать это в функцию, чтобы немного упростить задачу:

par(mfrow = c(2,2))
with(mtcars, {
  plot_lo(mpg, wt)
  plot_lo(mpg, wt, limits = c(1 / 1.5, 1.5))
  dd <<- plot_lo(mpg, wt, limits = c(1 / 1.2, 1.2))
  plot_lo(mpg, wt, pch = 16, las = 1, tcl = .5, bty = 'l')
})

str(dd)
# List of 2
# $ x: num [1:28] 21 21 22.8 21.4 18.7 18.1 14.3 22.8 19.2 17.8 ...
# $ y: num [1:28] 2.62 2.88 2.32 3.21 3.44 ...

введите здесь описание изображения

plot_lo <- function(x, y, limits = c(-Inf, Inf), ...) {
  lo <- loess.smooth(x, y)
  fx <- approxfun(lo$x, lo$y * limits[1L])
  fy <- approxfun(lo$x, lo$y * limits[2L])

  idx <- which(y < fx(x) | y > fy(x))
  if (length(idx)) {
    x  <- x[-idx]
    y  <- y[-idx]
    lo <- loess.smooth(x, y)
  }

  op <- par(..., no.readonly = TRUE)
  on.exit(par(op))

  plot(x, y)
  lines(lo$x, lo$y, lwd = 3)
  lines(lo$x, lo$y * limits[1L], lwd = 3, col = 2L)
  lines(lo$x, lo$y * limits[2L], lwd = 3, col = 2L)

  invisible(list(x = x, y = y))
}
person rawr    schedule 19.01.2016
comment
Не является ли коэффициент x1.2 здесь несколько произвольным? - person user2398029; 19.01.2016
comment
@ user2398029 ты спрашиваешь меня или оператора? - person rawr; 19.01.2016
comment
Я только что заметил, что у ОП был этот фактор в его вопросе, поэтому я думаю, что ответ следует. Нет ли способа установить этот коэффициент непроизвольным образом? - person user2398029; 19.01.2016
comment
@user2398029 user2398029 Я думаю, что lo$y * limits[1] можно заменить чем угодно, используя quantile или sd для значений y вместо простого масштабирования. Прелесть статистики в том, что вы можете делать все что угодно! пока вы можете оправдать это - person rawr; 19.01.2016

Я предлагаю пойти посмотреть на outliers package. Этот пакет позволяет провести идентификацию до проведения анализа. Это очень ПРОСТОЙ пример:

library(outliers)
series<-c(runif(100,1,2),1000)
round(scores(series,prob=1,type="chisq"),3)

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

series<-series[which(series<0.95),]
person Hanjo Odendaal    schedule 22.01.2016