1

我需要确定位于给定矩形内的一组 xy 坐标的分数。该矩形定义为边距坐标系边缘给定距离的区域（在这种情况下，坐标系的边界大致为 (-50, -20), (-50, 20), (50, 20 ), (50, -20)。另外，我希望能够在距边缘不同距离的矩形上测试结果。我的方法如下：

  # set initial limits to the coordinate system
  lim.xleft = -50
  lim.xright = 50
  lim.ybottom = -20
  lim.ytop = 20

frac.near.edge <- function(coord.pairs, tolerance){
  # set the coordinates of the rectangle of interest
  exclude.xleft = lim.xleft + tolerance
  exclude.xright = lim.xright - tolerance
  exclude.ybottom = lim.ybottom + tolerance
  exclude.ytop = lim.ytop - tolerance
  out <- vector()
  # loop through the pairs testing whether the point is inside the rectangle or outside
  for(i in 1:nrow(coord.pairs)){
    if(coord.pairs[i, 1] > exclude.xleft & coord.pairs[i, 1] < exclude.xright &  coord.pairs[i, 2] > exclude.ybottom & coord.pairs[i, 2] < exclude.ytop){
      out[i] <- "in"
    } else {
      out[i] <- "out"
    }
  }
  # return how many points were inside the rectangle and how many were outside
  return(table(out))
}

# try it out on something much bigger!
foo <- data.fram(x = runif(100), y = runif(100))
system.time(frac.near.edge(foo,  tolerance = 5))

这对于大型数据集来说非常慢（我的数据集包含大约 10^5 xy 对）。我怎样才能加快速度？绕过循环的方法？

4

2 回答 2

1

这可能更适合 SE 代码审查 ( https://codereview.stackexchange.com/questions/tagged/r )。我不知道这段编码是否有帮助，但你的问题真的不是关于编程，而是更多的代码改进。我还生成了一个更好的数据集，因为你们都生成了一个响应。

foo <- data.frame(x = sample(-100:100, 100, replace=TRUE), 
    y = sample(-100:100, 100, replace=TRUE))

xleft = -50
xright = 50
ybottom = -2
ytop = 20
foo$x >= xleft & foo$x <= xright  & foo$y >= ybottom & foo$y <= ytop
于 2012-04-30T21:17:34.437 回答
1 
exclude.xleft = lim.xleft + tolerance
exclude.xright = lim.xright - tolerance
exclude.ybottom = lim.ybottom + tolerance
exclude.ytop = lim.ytop - tolerance
out <- c("out", "in")[1+( findInterval(coord.pairs[ , 1], c(exclude.xleft, exclude.xright))==1 &
                       findInterval(coord.pairs[ , 2], c(exclude.ybottom,   exclude.ytop))==1)]

在 50K 的测试用例上，您的方法需要 0.01 秒而不是 19 秒：

coord.pairs<- cbind(rnorm(50000, 0,50), rnorm(1000,0,20)); tolerance=10
于 2012-04-30T21:29:54.347 回答