根据阈值在R中将向量分段,以确定值是否高于阈值。

8

我有一个长向量,需要根据阈值将其分成若干段。一段是连续的超过阈值的值。当值低于阈值时,该段结束并且下一段从值再次超过阈值的地方开始。我需要记录每个段的起始和结束索引。

以下是一个效率低下的实现。如何编写最快和最合适的方式?这很丑陋,我不得不假设有一个更简洁的实现方式。

set.seed(10)
test.vec <- rnorm(100, 8, 10)
threshold <- 0
segments <- list()
in.segment <- FALSE
for(i in 1:length(test.vec)){

    # If we're in a segment
    if(in.segment){
        if(test.vec[i] > threshold){
            next
        }else{
            end.ind <- i - 1
            in.segment <- FALSE
            segments[[length(segments) + 1]] <- c(start.ind, end.ind)
        }
    }

    # if not in segment
    else{
        if(test.vec[i] > threshold){        
            start.ind <- i
            in.segment <- TRUE
        }
    }
}

编辑:所有解决方案的运行时间

感谢所有回复,这对我很有帮助,也很有教育意义。下面是五种解决方案的小测试(包括提供的四种和原始示例)。正如您所看到的,所有四种解决方案都比原始解决方案要快得多,但Khashaa的解决方案是迄今为止最快的。

set.seed(1)
test.vec <- rnorm(1e6, 8, 10);threshold <- 0

originalFunction <- function(x, threshold){
    segments <- list()
    in.segment <- FALSE
    for(i in 1:length(test.vec)){

    # If we're in a segment
        if(in.segment){
            if(test.vec[i] > threshold){
                next
            }else{
                end.ind <- i - 1
                in.segment <- FALSE
                segments[[length(segments) + 1]] <- c(start.ind, end.ind)
            }
        }

    # if not in segment
        else{
            if(test.vec[i] > threshold){        
                start.ind <- i
                in.segment <- TRUE
            }
        }
    }
    segments
}

SimonG <- function(x, threshold){

  hit <- which(x > threshold)
  n <- length(hit)

  ind <- which(hit[-1] - hit[-n] > 1)

  starts <- c(hit[1], hit[ ind+1 ])
  ends <- c(hit[ ind ], hit[n])

  cbind(starts,ends)
}

Rcpp::cppFunction('DataFrame Khashaa(NumericVector x, double threshold) {
  x.push_back(-1);
  int n = x.size(), startind, endind; 
  std::vector<int> startinds, endinds;
  bool insegment = false;
  for(int i=0; i<n; i++){
    if(!insegment){
      if(x[i] > threshold){        
        startind = i + 1;
        insegment = true;          }
    }else{
      if(x[i] < threshold){
        endind = i;
        insegment = false;
        startinds.push_back(startind); 
        endinds.push_back(endind);
      }
    }
  }
  return DataFrame::create(_["start"]= startinds, _["end"]= endinds);
}')

bgoldst <- function(x, threshold){
    with(rle(x>threshold),
         t(matrix(c(0L,rep(cumsum(lengths),2L)[-length(lengths)]),2L,byrow=T)+1:0)[values,])   
}

ClausWilke <- function(x, threshold){
    suppressMessages(require(dplyr, quietly = TRUE))
    in.segment <- (x > threshold)
    start <- which(c(FALSE, in.segment) == TRUE & lag(c(FALSE, in.segment) == FALSE)) - 1
    end <- which(c(in.segment, FALSE) == TRUE & lead(c(in.segment, FALSE) == FALSE))
    data.frame(start, end)    
}

system.time({ originalFunction(test.vec, threshold); })
 ## user  system elapsed 
 ## 66.539   1.232  67.770 
system.time({ SimonG(test.vec, threshold); })
 ## user  system elapsed 
 ## 0.028   0.008   0.036 
system.time({ Khashaa(test.vec, threshold); })
 ## user  system elapsed 
 ## 0.008   0.000   0.008 
system.time({ bgoldst(test.vec, threshold); })
 ## user  system elapsed 
 ## 0.065   0.000   0.065 
system.time({ ClausWilke(test.vec, threshold); })
 ## user  system elapsed 
 ## 0.274   0.012   0.285 
4个回答

6

这里有另一个选项,主要使用 which。通过查找 hit 序列中的非连续元素来确定起始点和结束点。

test.vec <- rnorm(100, 8, 10)
threshold <- 0


findSegments <- function(x, threshold){

  hit <- which(x > threshold)
  n <- length(hit)

  ind <- which(hit[-1] - hit[-n] > 1)

  starts <- c(hit[1], hit[ ind+1 ])
  ends <- c(hit[ ind ], hit[n])

  cbind(starts,ends)

}

findSegments(test.vec, threshold=0)

这将得到类似以下的结果:
> findSegments(test.vec, threshold=0)
      starts ends
 [1,]      1    3
 [2,]      5    7
 [3,]      9   11
 [4,]     13   28
 [5,]     30   30
 [6,]     32   32
 [7,]     34   36
 [8,]     38   39
 [9,]     41   41
[10,]     43   43
[11,]     46   51
[12,]     54   54
[13,]     56   61
[14,]     63   67
[15,]     69   72
[16,]     76   77
[17,]     80   81
[18,]     83   84
[19,]     86   88
[20,]     90   92
[21,]     94   95
[22,]     97   97
[23,]    100  100

将其与原始序列进行比较:

> round(test.vec,1)
  [1]  20.7  15.7   4.3 -15.1  24.6   9.4  23.2  -4.5  16.9  20.9  13.2  -1.2
 [13]  22.6   7.7   6.0   6.6   4.1  21.3   5.3  16.7  11.4  16.7  19.6  16.7
 [25]  11.6   7.3   3.7   8.4  -4.5  11.7  -7.1   8.4 -18.5  12.8  22.5  11.0
 [37]  -3.3  11.1   6.9  -7.9  22.9  -3.7   3.5  -7.1  -5.9   3.5  13.2  20.0
 [49]  13.2  23.4  15.9  -5.0  -6.3  10.0  -6.2   4.7   2.1  26.4   5.9  27.3
 [61]  14.3 -12.4  28.4  30.9  18.2  11.4   5.7  -4.5   6.2  12.0  10.9  11.1
 [73]  -2.0  -9.0  -1.4  15.4  19.1  -1.6  -5.4   5.4   7.8  -5.6  15.2  13.8
 [85] -18.8   7.1  17.1   9.3  -3.9  22.6   1.7  28.9 -21.3  21.2   8.2 -15.4
 [97]   3.2 -10.2  -6.2  14.1

5

我喜欢使用 for 循环,因为将其翻译到 Rcpp 很简单。

Rcpp::cppFunction('DataFrame findSegment(NumericVector x, double threshold) {
  x.push_back(-1);
  int n = x.size(), startind, endind; 
  std::vector<int> startinds, endinds;
  bool insegment = false;
  for(int i=0; i<n; i++){
    if(!insegment){
      if(x[i] > threshold){        
        startind = i + 1;
        insegment = true;          }
    }else{
      if(x[i] < threshold){
        endind = i;
        insegment = false;
        startinds.push_back(startind); 
        endinds.push_back(endind);
      }
    }
  }
  return DataFrame::create(_["start"]= startinds, _["end"]= endinds);
}')
set.seed(1); test.vec <- rnorm(1e7,8,10); threshold <- 0;
system.time(findSegment(test.vec, threshold))

#   user  system elapsed 
#  0.045   0.000   0.045 

# @SimonG's solution
system.time(findSegments(test.vec, threshold))
#   user  system elapsed 
#  0.533   0.012   0.548 

4
with(rle(test.vec>threshold),t(matrix(c(0L,rep(cumsum(lengths),2L)[-length(lengths)]),2L,byrow=T)+1:0)[values,]);
##       [,1] [,2]
##  [1,]    1    8
##  [2,]   10   13
##  [3,]   16   17
##  [4,]   20   26
##  [5,]   28   28
##  [6,]   30   34
##  [7,]   36   38
##  [8,]   41   46
##  [9,]   48   49
## [10,]   51   53
## [11,]   55   81
## [12,]   84   90
## [13,]   92  100

翻译:解释
test.vec>threshold
##  [1]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE FALSE FALSE  TRUE  TRUE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE

使用向量化比较计算输入向量中哪些元素超过了阈值。
rle(...)
## Run Length Encoding
##   lengths: int [1:25] 8 1 4 2 2 2 7 1 1 1 ...
##   values : logi [1:25] TRUE FALSE TRUE FALSE TRUE FALSE ...

计算逻辑向量的运行长度编码。它返回一个列表,类别为'rle',其中包含两个命名组件:lengths,包含每个运行长度的长度,和values,包含运行该长度的值,这种情况下将是TRUEFALSE,前者表示感兴趣的段落,后者表示非段落运行长度。
with(...,...)

第一个参数是如上所述的运行长度编码。这将在虚拟环境中评估第二个参数,该环境由 'rle' 类列表组成,从而使得 lengthsvalues 组件可以作为词法变量访问。
下面我会深入介绍第二个参数的内容。
cumsum(lengths)
##  [1]   8   9  13  15  17  19  26  27  28  29  34  35  38  40  46  47  49  50  53  54  81  83  90  91 100

计算 lengths 的累积和。这将成为计算每个运行长度的起始索引和结束索引的基础。关键点:cumsum 的每个元素表示该运行长度的结束索引。
rep(...,2L)
##  [1]   8   9  13  15  17  19  26  27  28  29  34  35  38  40  46  47  49  50  53  54  81  83  90  91 100   8   9  13  15  17  19  26  27  28  29  34  35  38  40  46  47  49  50  53  54  81  83  90  91 100

复制累加和。第一次重复将作为起始索引的基础,第二次重复将作为结束索引的基础。从此我将把这些重复称为“起始索引重复”和“结束索引重复”。
c(0L,...[-length(lengths)])
##  [1]   0   8   9  13  15  17  19  26  27  28  29  34  35  38  40  46  47  49  50  53  54  81  83  90  91   8   9  13  15  17  19  26  27  28  29  34  35  38  40  46  47  49  50  53  54  81  83  90  91 100

这将删除开始索引重复结尾处的最后一个元素,并在其开头添加零。这实际上将开始索引重复向后延迟了一个元素。这是必要的,因为我们需要通过将前一个运行长度的结束索引加一来计算每个开始索引,将零作为第一个不存在的运行长度的结束索引。
matrix(...,2L,byrow=T)
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## [1,]    0    8    9   13   15   17   19   26   27    28    29    34    35    38    40    46    47    49    50    53    54    81    83    90    91
## [2,]    8    9   13   15   17   19   26   27   28    29    34    35    38    40    46    47    49    50    53    54    81    83    90    91   100

这将根据先前的结果构建一个两行矩阵。延迟的起始索引重复是顶行,结束索引重复是底行。
...+1:0
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25]
## [1,]    1    9   10   14   16   18   20   27   28    29    30    35    36    39    41    47    48    50    51    54    55    82    84    91    92
## [2,]    8    9   13   15   17   19   26   27   28    29    34    35    38    40    46    47    49    50    53    54    81    83    90    91   100

R循环将这两个元素加数首先沿行移动,然后沿列移动,这样可以将一添加到顶行。这完成了起始索引的计算。
t(...)
##       [,1] [,2]
##  [1,]    1    8
##  [2,]    9    9
##  [3,]   10   13
##  [4,]   14   15
##  [5,]   16   17
##  [6,]   18   19
##  [7,]   20   26
##  [8,]   27   27
##  [9,]   28   28
## [10,]   29   29
## [11,]   30   34
## [12,]   35   35
## [13,]   36   38
## [14,]   39   40
## [15,]   41   46
## [16,]   47   47
## [17,]   48   49
## [18,]   50   50
## [19,]   51   53
## [20,]   54   54
## [21,]   55   81
## [22,]   82   83
## [23,]   84   90
## [24,]   91   91
## [25,]   92  100

将其转换为二列矩阵。如果您可以接受将结果作为二行矩阵获得,这并非完全必要。
...[values,]
##       [,1] [,2]
##  [1,]    1    8
##  [2,]   10   13
##  [3,]   16   17
##  [4,]   20   26
##  [5,]   28   28
##  [6,]   30   34
##  [7,]   36   38
##  [8,]   41   46
##  [9,]   48   49
## [10,]   51   53
## [11,]   55   81
## [12,]   84   90
## [13,]   92  100

只选择感兴趣的片段。由于values是表示超过阈值的重复长度的逻辑向量,因此我们可以直接将其用作行索引向量。

性能

我猜我自己在这里搞砸了,但SimonG的解决方案的表现大约是我的两倍:

bgoldst <- function() with(rle(test.vec>threshold),t(matrix(c(0L,rep(cumsum(lengths),2L)[-length(lengths)]),2L,byrow=T)+1:0)[values,]);
simong <- function() findSegments(test.vec,threshold);
set.seed(1); test.vec <- rnorm(1e7,8,10); threshold <- 0;
identical(bgoldst(),unname(simong()));
## [1] TRUE
system.time({ bgoldst(); })
##    user  system elapsed
##   1.344   0.204   1.551
system.time({ simong(); })
##    user  system elapsed
##   0.656   0.109   0.762

我给你点赞(+1)...

3

我认为下面的解决方案更简单。请注意,设置随机数生成器的种子时,必须使用 set.seed(10),而不是 set.seed <- 10

require(dplyr) # for lead() and lag()

set.seed(10)
test.vec <- rnorm(100, 8, 10)
threshold <- 0

in.segment <- (test.vec > threshold)
start <- which(c(FALSE, in.segment) == TRUE & lag(c(FALSE, in.segment) == FALSE)) - 1
end <- which(c(in.segment, FALSE) == TRUE & lead(c(in.segment, FALSE) == FALSE))
segments <- data.frame(start, end)

head(segments)
##   start end
## 1     1   2
## 2     4   6
## 3     8   8
## 4    10  16
## 5    18  21
## 6    23  23

一般来说,在R语言中,如果你发现自己要写复杂的循环和条件语句,那么很可能是方法不对。大多数问题都可以用向量化的方式解决。


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接