查找出现在另一个向量值范围内的向量值答案

【问题标题】：Finding values of vector that occur within range of another vector's values查找出现在另一个向量值范围内的向量值
【发布时间】：2015-12-19 02:44:05
【问题描述】：

我有两个序列。它们是以秒为单位的时间。我想知道序列 b 中的哪些值出现在序列 a 中任何值的 10 秒内。

seqa = c(4.53333333333333, 7.43333333333334, 9.03333333333333, 20.6166666666667, 
20.6333333333333, 42.5666666666667, 48.3166666666667, 48.8, 49.75, 
55.1, 56.7833333333333, 59.3833333333333, 110.15, 113.95, 114.6)

seqb = c(18.3833333333333, 18.3833333333333, 63.8833333333333, 72.3166666666667, 
76.7166666666667, 85.2166666666667, 91.25, 91.3666666666667, 
96.2833333333333)

我使用两个for 循环完成了这项工作。遍历seqb 的每个元素并测试它是否出现在大于seqa 的每个值但在10 秒限制内的时间。

matX <- matrix(nrow=length(seqa), ncol=length(seqb))

for(j in seq_along(seqb)){
  for(i in seq_along(seqa)){
    test1 <- seqb[j]>=seqa[i]
    test2 <- seqb[j]<=seqa[i]+10
    matX[i,j] <- sum(test1 + test2)
  }
}
matX

我将结果存储在矩阵中。您可以在第 1、2 和 3 列中看到 2 的值。

      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
 [1,]    1    1    1    1    1    1    1    1    1
 [2,]    1    1    1    1    1    1    1    1    1
 [3,]    2    2    1    1    1    1    1    1    1
 [4,]    1    1    1    1    1    1    1    1    1
 [5,]    1    1    1    1    1    1    1    1    1
 [6,]    1    1    1    1    1    1    1    1    1
 [7,]    1    1    1    1    1    1    1    1    1
 [8,]    1    1    1    1    1    1    1    1    1
 [9,]    1    1    1    1    1    1    1    1    1
[10,]    1    1    2    1    1    1    1    1    1
[11,]    1    1    2    1    1    1    1    1    1
[12,]    1    1    2    1    1    1    1    1    1
[13,]    1    1    1    1    1    1    1    1    1
[14,]    1    1    1    1    1    1    1    1    1
[15,]    1    1    1    1    1    1    1    1    1

out <- apply(matX, 2, function(x) any(x>=2))    
seqb[out]

# [1] 18.38333 18.38333 63.88333

这些值是在seqa 中至少一个值的 10 秒内出现的值。（前两个出现在 9.03333 的 10 秒内，第三个值 63.8333 出现在三个 seqa 值（55.1、56.78333、59.38333）的 10 秒内。

我正在尝试加快速度，因为我将对大约 2000 个元素的序列进行一些随机化。任何想法都非常感谢。

【问题讨论】：

标签： r

【解决方案1】：

这里有两个基本选项

seqa = c(4.53333333333333, 7.43333333333334, 9.03333333333333, 20.6166666666667, 
         20.6333333333333, 42.5666666666667, 48.3166666666667, 48.8, 49.75, 
         55.1, 56.7833333333333, 59.3833333333333, 110.15, 113.95, 114.6)

seqb = c(18.3833333333333, 18.3833333333333, 63.8833333333333, 72.3166666666667, 
         76.7166666666667, 85.2166666666667, 91.25, 91.3666666666667, 
         96.2833333333333)


## via alexis_laz
a <- function() seqb[seqa[findInterval(seqb, seqa)] + 10 >= seqb]
# [1] 18.38333 18.38333 63.88333


## f
(function() {
  la <- length(seqa)
  lb <- length(seqb)
  rr <- rep(seqb, each = la)
  m <- matrix(rep(seqa, length(seqb)) - rr, la)
  +(m < 0 & abs(m) <= 10)
})()

## g
o <- outer(seqa, seqb, `-`)
x <- +(o < 0 & abs(o) <= 10)

`dimnames<-`(x, list(round(seqa, 2), round(seqb, 2)))

#        18.38 18.38 63.88 72.32 76.72 85.22 91.25 91.37 96.28
# 4.53       0     0     0     0     0     0     0     0     0
# 7.43       0     0     0     0     0     0     0     0     0
# 9.03       1     1     0     0     0     0     0     0     0
# 20.62      0     0     0     0     0     0     0     0     0
# 20.63      0     0     0     0     0     0     0     0     0
# 42.57      0     0     0     0     0     0     0     0     0
# 48.32      0     0     0     0     0     0     0     0     0
# 48.8       0     0     0     0     0     0     0     0     0
# 49.75      0     0     0     0     0     0     0     0     0
# 55.1       0     0     1     0     0     0     0     0     0
# 56.78      0     0     1     0     0     0     0     0     0
# 59.38      0     0     1     0     0     0     0     0     0
# 110.15     0     0     0     0     0     0     0     0     0
# 113.95     0     0     0     0     0     0     0     0     0
# 114.6      0     0     0     0     0     0     0     0     0

我的破硬件上有一些长凳

library('microbenchmark')
seqa <- rep(seqa, 100)
seqb <- rep(seqb, 100)
microbenchmark(f(), g(), baseR(), DT(), unit = 'relative')
# Unit: relative
#      expr        min         lq       mean    median         uq       max neval  cld
#       f()   525.3178  374.23871  402.51609  386.4717  372.50657  496.6496   100   c 
#       g()   293.2158  223.21560  247.40211  241.3430  225.80202  443.5323   100  bc 
#   baseR() 13268.9357 9357.70517 8895.30834 9111.6828 8466.15623 6702.1735   100    d
#      DT()   136.1109   93.61985   96.88054   96.0771   95.03329  100.5602   100 ab  
#       a()     1.0000    1.00000    1.00000    1.0000    1.00000    1.0000   100 a

【讨论】：

如果“seqa”按看起来排序并且不需要中间矩阵，另一种方法-除非我遗漏了示例中不明显的内容-可能是seqb[seqa[findInterval(seqb, seqa)] + 10 >= seqb]以避免将一切与一切进行比较。
@alexis_laz 聪明！迄今为止最快的超级骗子
@alexis_laz 应该适用于所有不同长度 >1 的序列 - 尝试对其他一些不同长度的序列对，我收到以下消息：Warning message: In seqa[findInterval(seqb, seqa)] + dt >= seqb : longer object length is not a multiple of shorter object length
@jalapic 是 seqa 的最小值
如果您想将a() 与其他解决方案进行比较，它是否也应该构造矩阵。

【解决方案2】：

您可以尝试data.table 包中的foverlaps 函数。

library(data.table)

b <- data.table(seqb)
a <- data.table(seqa)
a[, end := seqa + 10]
setkey(a)
b[, end := seqb]

inds <- foverlaps(b, a,
                  by.x=c("seqb","end"), 
                  type="within",
                  mult="all",
                  which=TRUE # you can use nomatch=0L, but it doesn't change the final matrix
                 )
 #   xid yid
 #1:   1   3
 #2:   2   3
 #3:   3  10
 #4:   3  11
 #5:   3  12
 #6:   4  NA
 #7:   5  NA
 #8:   6  NA
 #9:   7  NA
#10:   8  NA
#11:   9  NA

这些索引现在可以用来创建你想要的矩阵。

mat <- matrix(1, nrow=length(seqa), ncol=length(seqb))
mat[cbind(inds$yid, inds$xid)] <- 2

这是在一个带有seqa 和seqb 硬代码的函数中：

DT <- function(){
    b <- data.table(seqb)
    a <- data.table(seqa)
    a[, end := seqa + 10]
    setkey(a)
    b[, end := seqb]

    inds <- foverlaps(b, a,
                      by.x=c("seqb","end"), 
                      type="within",
                      mult="all",
                      which=TRUE 
                     )

    mat <- matrix(1, nrow=length(seqa), ncol=length(seqb))
    mat[cbind(inds$yid, inds$xid)] <- 2
    mat
}

【讨论】：

【解决方案3】：

seqa = c(4.53333333333333, 7.43333333333334, 9.03333333333333, 20.6166666666667, 20.6333333333333, 42.5666666666667, 48.3166666666667, 48.8, 49.75, 55.1, 56.7833333333333, 59.3833333333333, 110.15, 113.95, 114.6)

seqb = c(18.3833333333333, 18.3833333333333, 63.8833333333333, 2.3166666666667, 76.7166666666667, 85.2166666666667, 91.25, 91.3666666666667, 96.2833333333333)

上面读取的数据。下面，我展示了我的方法，以及@jota 的方法。请注意，这是一个有点愚蠢的比较，因为数据太小了。 data.table 解决方案几乎可以肯定在处理更大的数据时要快得多。

library(microbenchmark)
library(data.table)

DT <- function(){
   b <- data.table(seqb)
   a <- data.table(seqa)
   a[, end := seqa + 10]
   setkey(a)
   b[, end := seqb]

   inds <- foverlaps(b, a,
                     by.x=c("seqb","end"), 
                     type="within",
                     mult="all",
                     which=TRUE 
                    )

   mat <- matrix(1, nrow=length(seqa), ncol=length(seqb))
   mat[cbind(inds$yid, inds$xid)] <- 2
   mat
}



baseR <- function(){
    out <- matrix(NA, ncol=length(seqb), nrow=length(seqa));
    for(i in 1:length(seqa)){
        out[i,] <- sapply(seqb, function(x){seqa[i] -10 < x  & x < seqa[i] +10})
    }
    out
}


microbenchmark(
    baseR(), DT()
)

以及微基准测试的结果（为了好玩）。

Unit: microseconds
    expr      min       lq     mean   median        uq      max neval
 baseR()  703.382  750.129  786.283  770.867  788.3085 1905.357   100
    DT() 7289.433 7415.906 7631.574 7503.236 7575.7345 8794.439   100

【讨论】：

如果您将sapply 中的函数更改为seqa[i] < x & x < seqa[i] +10，您将匹配jalapic 列出的输出。仅供参考，我稍微更改了我的 data.table 答案，因此您发布为 datTable() 的答案略有不同。

【解决方案4】：

您可以使用IRanges 包。

library(IRanges)

a.ir <- IRanges(round(seqa, 4)*1e4, round(seqa, 4)*1e4+10*1e4)
b.ir <- IRanges(round(seqb, 4)*1e4, round(seqb, 4)*1e4)

findOverlaps(b.ir, a.ir)
# Hits of length 5
# queryLength: 9
# subjectLength: 15
#   queryHits subjectHits 
#    <integer>   <integer> 
# 1         1           3 
# 2         2           3 
# 3         3          10 
# 4         3          11 
# 5         3          12 

seqb[unique(queryHits(findOverlaps(b.ir, a.ir)))]
# [1] 18.38333 18.38333 63.88333

【讨论】：

你可以直接使用seqb[countOverlaps(b.ir, a.ir) > 0]，而不是findOverlaps；它也可能更快。但是，老实说，我觉得随意取整浮点值并代表它们进行计算有点容易出错？
countOverlaps 也可以。由于 IRange 只接受整数，因此您必须将这些值更改为整数。 OP应该没问题。