Запомните последнее правильное значение последовательности (для удаления выбросов)

У меня есть небольшая проблема в функции. Целью этого является удаление выбросов, которые я обнаружил в своем data.frame. Они обнаруживаются, когда есть слишком большая разница с предыдущим правильным значением (например, c(1,2,3,20,30,4,5,6): "20" и "30" являются выбросами). Но мои данные намного сложнее, чем это.

Моя идея состоит в том, чтобы считать первые два числовых значения моего столбца «правильными». Затем я хочу проверить каждое следующее значение:

  • если разница между тестируемым значением и предыдущим составляет ‹20, то это новое правильное значение, и тест должен начинаться снова с этого нового правильного значения (а не с предыдущего правильного)
  • если та же разница больше 20, то это неправильная. Индекс должен быть помещен рядом с неправильным значением, и тест должен продолжаться с того же правильного значения, пока не будет обнаружено новое правильное значение.

Вот пример с моей функцией и поддельным DF:

myts <- data.frame(x=c(12,12,35,39,46,45,33,5,26,28,29,34,15,15),z=NA) 

test <- function(x){
st1 = NULL
temp <- st1[1] <- x[1]
st1 <- numeric(length(x))
for (i in 2:(length(x))){ 
    if((!is.na(x[i])) & (!is.na(x[i-1]))& (abs((x[i])-(temp)) > 20)){
st1[i] <- 1
} } 
return(st1)
}

myts[,2] <- apply(as.data.frame(myts[,1]),2,test)  
myts[,2] <- as.numeric(myts[,2])

Он выполняет почти всю работу, но проблема в том, что последнее правильное значение не запоминается. Он по-прежнему выполняет тест с первого правильного значения. Из-за этого строки с 9 по 11 в моем примере не обнаруживаются. Я позволю вам представить проблему на 500 000 строк data.frame.

Как мне решить эту маленькую проблему? Остальные функции могут быть в порядке.


r function outliers
person jeff6868    schedule 03.03.2015    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Вам просто нужно обновить temp для любых индексов, которые не являются выбросами:

test <- function(x) {
  temp <- x[1]
  st1 <- numeric(length(x))
  for (i in 2:(length(x))){ 
    if(!is.na(x[i]) & !is.na(x[i-1]) & abs(x[i]-temp) > 20) {
      st1[i] <- 1
    } else {
      temp <- x[i]
    }
  } 
  return(st1)
}

myts[,2] <- apply(as.data.frame(myts[,1]),2,test)  
myts[,2] <- as.numeric(myts[,2])
myts
#     x z
# 1  12 0
# 2  12 0
# 3  35 1
# 4  39 1
# 5  46 1
# 6  45 1
# 7  33 1
# 8   5 0
# 9  26 1
# 10 28 1
# 11 29 1
# 12 34 1
# 13 15 0
# 14 15 0

Следует отметить, что циклы for в R будут довольно медленными по сравнению с векторизованными функциями. Однако, поскольку каждый элемент в вашем векторе сложным образом зависит от предыдущих, сложно использовать встроенные векторизованные функции R для эффективного вычисления вашего вектора. Вы можете преобразовать этот код почти дословно в C++ и использовать пакет Rcpp для восстановления эффективности:

library(Rcpp)
test2 <- cppFunction(
"IntegerVector test2(NumericVector x) {
  const int n = x.length();
  IntegerVector st1(n, 0);
  double temp = x[0];
  for (int i=1; i < n; ++i) {
    if (!R_IsNA(x[i]) && !R_IsNA(x[i]) && fabs(x[i] - temp) > 20.0) {
      st1[i] = 1;
    } else {
      temp = x[i];
    }
  }
  return st1;
}")
all.equal(test(myts[,1]), test2(myts[,1]))
# [1] TRUE

# Benchmark on large vector with some NA values:
set.seed(144)
large.vec <- c(0, sample(c(1:50, NA), 1000000, replace=T))
all.equal(test(large.vec), test2(large.vec))
# [1] TRUE
library(microbenchmark)
microbenchmark(test(large.vec), test2(large.vec))
# Unit: milliseconds
#              expr         min          lq       mean     median         uq        max neval
#   test(large.vec) 2343.684164 2468.873079 2667.67970 2604.22954 2747.23919 3753.54901   100
#  test2(large.vec)    9.596752    9.864069   10.97127   10.23011   11.68708   16.67855   100

Код Rcpp примерно в 250 раз быстрее на векторе длиной 1 миллион. В зависимости от вашего варианта использования это ускорение может быть или не быть важным.

person josliber♦    schedule 03.03.2015
comment
В моем случае моей функции будет достаточно. Для вычисления всего моего вектора требуется около 2 секунд. Но я вдохновлюсь этим для других циклов в моем коде, которые намного медленнее, чем эта функция. Еще раз спасибо Джозильбер! - person jeff6868; 04.03.2015