我有一个包含两组天数范围的数据集(一天被编码为一年中的天数).对于每一行,我想计算这些范围每月相应的天数.
在我的示例数据中,列’deb’和’fin’是每行中第一个子范围的开始和结束日期,’deb2’和’fin2’是第二个子范围的限制.
d <- data.frame(deb = c(1, 32, 90, 91), fin = c(31, 59, 91, 91),
deb2 = c(50, 0, 0, 0), fin2 = c(60, 0, 0, 0))
d
# deb fin deb2 fin2
#1 1 31 50 60
#2 32 59 0 0
#3 90 91 0 0
#4 91 91 0 0
例如,对于第1行,第一个范围(从“deb”到“fin”)从第1天到第31天,第二个范围从第50天到第60天.
在计算了两个范围的每月天数之后,我希望得到类似的结果:
# jan feb mar
#[1,] 31 10 1
#[2,] 0 28 0
#[3,] 0 0 2
#[4,] 0 0 1
(NAs而不是零不是问题)
我尝试了几个解决方案,如下面的三个(第三个’g3’是最快的),并尝试使用tidyverse,显示速度较慢.我想知道是否有最快的选择,因为在现实生活中我有很多行.问题似乎在于从范围到月份参考列表的转换,但也可能是计算方式.
f1<-function(deb,fin,deb2,fin2,...) {
f<-factor(c(deb:fin,deb2:fin2))
levels(f)<-list(jan=1:31,feb=32:59,mar=60:91)
table(f)
}
g1 <- function() do.call(rbind,d %>% pmap(f1))
K <- vector(10,mode="character")
K[1:31] <- "jan"; K[32:59] <- "feb"; K[60:91] <- "mar"
f2 <- Vectorize(function(deb,fin,deb2,fin2) table(c(K[deb:fin],K[deb2:fin2])))
g2 <- function() do.call(bind_rows,f2(d$deb,d$fin,d$deb2,d$fin2))
L <- K
names(L) <- 1:91
f3 <- Vectorize(function(deb,fin,deb2,fin2) c(L[deb:fin],L[deb2:fin2]))
g3 <- function() {
as.matrix(do.call(bind_rows,f3(d$deb,d$fin,d$deb2,d$fin2))) -> m
z <- unlist(map(list("jan","feb","mar"),
function(y) apply(m,1,function(x) sum(x==y,na.rm=TRUE))))
dim(z)<-c(nrow(d),3)
z
}
更新
一些基准测试如下.我在试验中添加了Chinsson12的解决方案,该解决方案在优雅的解决方案中表现良好.
firstOfMths <- seq(as.Date("2018-01-01"), as.Date("2019-01-01"), by="month")
daysPerMth <- c(1L, cumsum(as.integer(diff(firstOfMths))))
chinsoon12 <- function()
t(apply(d, 1, function(x)
table(cut(c(x["deb"]:x["fin"],x["deb2"]:x["fin2"]), daysPerMth, labels=month.abb, include.lowest=TRUE, right=TRUE))
))
N <- 500
d<-data.frame(deb=rep(c(1,32,90,91),N),fin=rep(c(31,59,91,91),N),deb2=rep(c(50,0,0,0),N),fin2=rep(c(60,0,0,0),N))
microbenchmark(g1(),g2(),g3(),chinsoon12())
#Unit: milliseconds
# expr min lq mean median uq max neval
# g1() 571.3890 615.1020 649.7619 639.6632 662.4808 976.9566 100
# g2() 306.7141 341.3056 360.9687 353.1227 373.8194 505.0882 100
# g3() 282.2767 304.4331 320.4908 314.2377 325.8846 543.4680 100
# chinsoon12() 429.7627 469.6998 500.6289 488.5176 512.0520 729.0995 100
最佳答案 使用findInterval,Map和table:
# create breaks to be used in findInterval
b <- <- as.numeric(format(seq(as.Date("2018-01-01"), as.Date("2018-12-31"), by = "month"), "%j"))
# use Map to expand the day of year ranges by row
# use findInterval to convert day of year to month number
# use the month numbers to index month.abb
l <- Map(function(from, to, from2, to2) month.abb[findInterval(c(from:to, from2:to2), b)], d$deb, d$fin, d$deb2, d$fin2)
# create a row index
i <- rep(1:nrow(d), lengths(l))
# use table to get a contigency table of row indices and months
table(i, factor(unlist(l), levels = month.abb))
# i Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
# 1 31 10 1 0 0 0 0 0 0 0 0 0
# 2 0 28 0 0 0 0 0 0 0 0 0 0
# 3 0 0 1 1 0 0 0 0 0 0 0 0
# 4 0 0 0 1 0 0 0 0 0 0 0 0
在较大的数据集上看起来比g3()更快(d < – d [rep(1:nrow(d),1e4),]).