假设我有一个100棵树的年龄向量.然后,我将这些树木老化5年、10年、15年和20年,以创建今年和future 四个5年规划期的树龄矩阵.

但随后,我决定砍伐其中一些树木(每个规划期仅10棵),记录在T/F值矩阵中,其中T是收获的,F不是(树木不能收获两次).

age.vec <- sample(x = 1:150, size = 100, replace = T) # create our trees
age.mat <- cbind(age.vec, age.vec+5, age.vec + 10, age.vec + 15, age.vec + 20) # grow them up
x.mat <- matrix(data = F, nrow = 100, ncol = 5) # create the empty harvest matrix
x.mat[cbind(sample(1:100, size = 50), rep(1:5, each = 10))] <-  T # 10 trees/year harvested

因此,当年收获的树木年龄变为零:

然后,我想在接下来的时间里再次对收获的树木进行老化.E、 g.如果一棵树在第一个计划期收获,在第二个计划期(5年后),我希望树的年龄为5岁,然后在第三个计划期(10年后),我希望树的年龄为10岁.我已在以下for循环中成功实现了这一点:

for (i in 2:5){ # we don't need to calculate over the first year
  age.mat[,i]<-age.mat[,i-1]+5L # add 5 to previous year
  age.mat[x.mat[,i],i] <- 0L # reset age of harvested trees to zero
}

然而,这是一种笨重而缓慢的工作方式.有没有办法更快地实现这一点(即没有for循环)?它也在函数中实现,这意味着使用"apply"实际上会减慢速度,因此需要直接对其进行矢量化.这是我迭代了数千次的东西,所以速度至关重要!

非常感谢.

推荐答案

乔恩·斯普林(JonS普林)的答案中,t(apply的替代品是matrixStats::rowCumsums.

library(matrixStats)

n <- 1e4L
n10 <- n/10L
age.mat <- outer(sample(150, n, TRUE), seq(0, 20, 5), "+")
x.mat <- matrix(FALSE, n, 5) # create the empty harvest matrix
# sample harvests so that no tree is harvested twice
x.mat[matrix(c(sample(n, n/2L), sample(n10:(6L*n10 - 1L)) %/% n10), n/2L)] <- TRUE

f1 <- function(age, x) {
  age[x[,1],] <- 0
  for (i in 2:5){ # we don't need to calculate over the first year
    age[,i] <- age[,i - 1] + 5L # add 5 to previous year
    age[x[,i], i] <- 0L # reset age of harvested trees to zero
  }
  age
}

f2 <- function(age, x) {
  age - rowCumsums(x*age)
}

microbenchmark::microbenchmark(f1 = f1(age.mat, x.mat),
                               f2 = f2(age.mat, x.mat),
                               check = "equal")
#> Unit: microseconds
#>  expr   min    lq     mean median     uq     max neval
#>    f1 294.4 530.2 1023.450  566.6 629.35 33222.8   100
#>    f2 135.2 263.6  334.622  284.2 307.15  4343.6   100

R相关问答推荐

显示矩阵的上半部分或下半部分时正确放置重要相关性

基于两个现有列创建新列

有没有方法将paste 0功能与列表结合起来?

更改网格的crs以匹配简单要素点对象的crs

使用Shiny组合和显示复制和粘贴的数据

使用scale_x_continuous复制ggplot 2中的离散x轴

为什么以及如何修复Mapview不显示所有点并且st_buffer合并一些区域R?

在R中查找每个组不同时间段的总天数

向gggplot 2中的数据和轴标签添加大写和星号

然后根据不同的列值有条件地执行函数

如何在R中合并两个基准点?

将包含卷的底部25%的组拆分为2行

将二进制数据库转换为频率表

方法::slotName如何处理非类、非字符的参数?

调换行/列并将第一行(原始数据帧的第一列)提升为标题的Tidyr类似功能?

R -如何分配夜间GPS数据(即跨越午夜的数据)相同的开始日期?

解析嵌套程度极高的地理数据

使用ggplot2中的sec_axis()调整次轴

ggplot斜体轴刻度标签中的单个字符-以前的帖子建议不工作

如何获取R chromote中的当前URL?