logo标签列表

寻找向量中的数值，这些数值落在另一个向量数值范围内。

r
4

我有两个序列，它们是以秒为单位的时间。 我希望知道序列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个元素的序列进行一些随机化。任何想法都非常感激。
4个回答
4

这里有两种基本选项

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
2
如果“seqa”按照它看起来的那样排序，并且如果不需要中间矩阵，另一种方法 - 除非我错过了在示例中不明显的东西 - 可以是seqb [seqa [findInterval（seqb，seqa）] +10> = seqb]，以避免将所有内容与所有内容进行比较。 - alexis_laz
1
@alexis_laz 聪明！而且绝对是最快的。 - rawr
@alexis_laz 这个方法适用于所有长度大于1的序列吗？当我尝试对一些不同长度的序列对进行操作时，出现了以下警告信息：Warning message: In seqa[findInterval(seqb, seqa)] + dt >= seqb : longer object length is not a multiple of shorter object length。 - jalapic
@jalapic是seqa的最小值<min seqb吗？如果不是，findinterval将会给你一些0。 - rawr
1
如果您想要将a()与其他解决方案进行比较，那么它也应该构建矩阵。 - Jota
显示剩余3条评论
1
你可以尝试使用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
}
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, 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()略有不同。 - Jota
0
你可以使用 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。这样做可能会更快。但说实话，我觉得任意舍入浮点值并代表它们进行计算有点容易出错。 - alexis_laz
countOverlaps也是可行的。由于IRanges仅接受整数，您必须将这些值更改为整数。对于OP来说应该没问题。 - Ven Yao
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接