代码之家 › 专栏 › 技术社区 › Maiasaura

当数据不连续时,如何随机抽取离y值最近的n个值?

statistics-bootstrap time-series r

Maiasaura · 技术社区 · 15 年前

我有一个数据集,其中包括物种列表、它们的数量以及调查开始时的天数。由于许多天没有取样,所以这一天是不连续的。例如,在第5天、第6天、第9天、第10天、第15天、第34天、第39天等等,可能会有鸟类被计算在内。我把最早的日期定为0天。

实例数据:

species     counts      day
Blue tit    234         0
Blue tit    24          5
Blue tit    45          6
Blue tit    32          9
Blue tit    6           10
Blue tit    98          15
Blue tit    40          34
Blue tit    57          39
Blue tit    81          43
..................

我需要引导这些数据,并获取一个结果数据集,我在其中指定何时开始、进行的时间间隔和要采样的点数。

例如:假设我随机选择第5天作为开始日,间隔为30,要采样的行数为2。这意味着我将从5开始,加上30,在35天左右寻找2行(但不是第35天本身)。在这种情况下,我将抓住第34和第39天的两行。

接下来我加30到35,找65左右的两点。漂洗,重复直到数据集结束。

我编写了这个函数来进行抽样,但它有缺陷(见下文):

resample <- function(x, ...) x[sample.int(length(x), ...)]
 locate_points<- function(dataz,l,n) #l is the interval, n is # points to sample. This is called by another function that specifies start time among other info.
{
   tlength=0
   i=1
    while(tlength<n)   
    {
        low=l-i
        high=l+i
        if(low<=min(dataz$day)) { low=min(dataz$day) }
        if(high>=max(dataz$day)) { high=max(dataz$day) }
        test=resample(dataz$day[dataz$day>low & dataz$day<high & dataz$day!=l])
          tlength=length(test)
         i=i+1
      } 
  test=sort(test)
  k=test[1:n]
 return (k)
 }

我需要帮助的两个问题:

虽然我的函数没有返回所需的点数,但它不是围绕我的搜索值居中的。这是有道理的,因为随着我的范围越来越广,我得到了更多的分数,当我对这些分数进行排序并选择第一个n时,它们往往不是低值。
第二,我该如何得到实际的行?现在,我有另一个函数来定位这些行,使用 which 然后 rbind 把这些排在一起。似乎应该有更好的方法。

谢谢!

2 回复 | 直到 11 年前

Charles 15 年前

像这样的怎么样:

day = 1:1000

search = seq(from=5, to=max(day), by=30)
x = sort(setdiff(day, search))
pos = match(x[unlist(lapply(findInterval(search, x), seq, len=2))], day)

day[pos]

要从data.frame中获取行,只需将其子集:

rows = data[pos, ]

这可能比unlist/lapply/seq组合稍微干净一点:

pos = match(x[outer(c(0, 1), findInterval(search, x), `+`)], day)

还要注意的是,如果你想要一个更大的窗口(比如说4),这只是一个返回一点的问题:

pos = match(x[outer(-1:2, findInterval(search, x), `+`)], day)

Joris Meys 15 年前

喜欢查尔斯的解决方案,它对n=2的情况非常有效。唉,它不能延伸到更大的窗户。它仍然存在OP所描述的问题:对于较大的窗口,选择并不是围绕搜索值进行的。考虑到n是偶数,我在很大程度上基于查尔斯的想法提出了以下解决方案。

函数控制边框。如果有100天,下一个中间点是最后一天的第二天,那么一个4的窗口将意味着您选择索引101,它给出 NA . 此函数移动窗口,使所有选定的索引都位于原始数据中。这也有副作用,取决于启动值( st (长度) l 和窗口( n )开始和结束的值被选择两次的可能性更大。长度应该至少是窗口大小的两倍。

函数的输出是引导样本的索引。它可以用作 pos 向量和数据帧上的charles变量。

bboot <- function(day,st,l,n){
  mid <- seq(st,max(day),by=l)
  x <-sort(setdiff(day,mid))
  lx <- length(x)

  id <- sapply(mid,
          function(y){
            m <- match(T,x>y)
            seq(
              from=min( lx-n, max(1,m+(-n/2)) ),
              to=min( lx, max(n,m+(n/2-1)) )
            )
          }
        )

  pos <- match(x[id],day)
  return(pos)
}

然后

>   day <- sample(1:100,50)
> sample.rownr <- bboot(day,10,20,6)
> sort(day)
 [1]  3  4  5  7  9 10 13 15 16 18 19 21 22 24 25 26 27 28 29 
[20] 30 31 32 35 36 38 40 45 49 51 52 54 55 58 59 62 65 69 72 73
[40] 74 80 84 87 88 91 92 94 97 98 99
> day[sample.rownr]
 [1]  5  7  9 13 15 16 27 28 29 31 32 35 40 45 49 51 52 54 62 
[20] 65 69 72 73 74 84 87 88 91 92 94
>

编辑:关于时间序列的引导,您应该 CRAN taskview on time series 尤其是关于重新采样的部分。对于不规则时间序列, zoo 软件包还提供了许多其他的功能,可以派上用场。