每月意味着申请多维数组

Question

我想计算多维数组的 3-D 平均值。由于这个维度应该是时间，我想计算每月的平均值。为此，我尝试使用 apply，但我不确定问题出在哪里。假设我的数据如下：

       #Creating a sample  
       m <-array(1:12, dim=c(20,4,36))
       #number of months
       months <- seq(1:12)
       #Compute the mean over each month (dimension of the result should be [20,4,12]
       monmean <- apply(m,1:2,function(x) for(i in 1:12) mean(x[,,months==i],na.rm=TRUE))

有什么想法吗？？ 提前致谢

Answer 1

我想我明白你在追求什么。这实际上比看起来要复杂一些，因为几个月不是固定的时间段。它们的天数会有所不同，并且由于闰年，2 月会因年份而异。因此，简单的常规逻辑或数字索引向量不足以精确计算此结果。您需要考虑阵列的 z 维度所涵盖的确切日期。

解决方案 1

您可以做的是单独计算一个日期向量，该向量标识与数组的每个 z-index 对应的日期。在每个 z 线的 apply() 调用中，您可以调用 strftime() 以提取每个此类日期的月份，并使用 tapply() 按该月份值分组以获取每月 mean()s。可以这样做：

set.seed(1);
R <- 48;
C <- 39;
Z <- 3653;
N <- R*C*Z;
a1 <- array(rnorm(N,10,2),c(R,C,Z));
dates <- seq(as.Date('2000-01-01'),as.Date('2009-12-31'),1);
a2 <- aperm(apply(a1,1:2,function(x) tapply(x,strftime(dates,'%m'),mean)),c(2,3,1));

这是一个演示，展示了一些具体的正确性证明：

for (r in sample(1:nrow(a2),2)) for (c in sample(1:ncol(a2),2)) for (m in sample(1:dim(a2)[3],2)) cat(sprintf('[%02d,%02d,%3s] %f %f\n',r,c,month.abb[m],mean(a1[r,c,strftime(dates,'%m')==sprintf('%02d',m)]),a2[r,c,m]));
## [14,05,Aug] 10.030313 10.030313
## [14,05,Apr] 10.200982 10.200982
## [14,25,Jan] 9.957879 9.957879
## [14,25,Apr] 10.185447 10.185447
## [26,34,Oct] 10.056931 10.056931
## [26,34,Nov] 9.876327 9.876327
## [26,17,Apr] 10.005423 10.005423
## [26,17,Sep] 10.009785 10.009785

备注

• 我随机选择了 2000-01-01 到 2009-12-31 的日期范围，因为它涵盖了 10 年期间（由于闰年）正好有 3653 天，但显然您应该确保使用您的真实数据实际涵盖的日期。
• 如您所见，调用apply() 并以1:2 作为边距，您走在正确的轨道上，因为这允许您在每条z 线上独立操作，这样您就可以按以下方式对z 线进行分组月并计算沿该 z 线的每个月的平均值。
• 不幸的是，apply() 有一个令人讨厌的习惯，即以不同于人们通常预期的换位方式返回结果。对于二维的使用，这通常通过简单的调用t() 来解决，但由于我们在这里处理三维，我们需要调用aperm() 来修复维度顺序。
• 由于我选择的日期从一月开始并按日历顺序逐月推进，因此结果中的方法最终将按日历月排序。 IOW，a2 中的 z-indexes 1:12 对应于 1 月至 12 月。如果您的日期不是从一月开始，那么这个解决方案应该仍然有效，但您必须注意结果中 z-indexes 和月份之间的对应关系。例如，我的“正确性证明”代码假定索引 1:12 对应于 Jan-Dec 月份，但如果月份在输入数组中以不同的顺序出现，那将是不正确的。

---

解决方案 2

在写这个答案时，我实际上想到了一个稍微不同的，并且可以说稍微好一点的解决方案。您可以只调用一次tapply()，然后按行分组，然后按列分组，最后按月分组。不幸的是，tapply() 似乎没有被设计为自然循环其组向量以覆盖输入向量，因此我们必须使用精心设计的对 rep() 的调用（使用 each 和 times 参数）自己循环它们仔细——我想tapply() 实际上甚至不知道如何正确处理我们的输入数据），但除此之外，它相当简单：

a3 <- tapply(a1,list(rep(1:R,C*Z),rep(1:C,each=R,times=Z),rep(strftime(dates,'%m'),each=R*C)),mean);

这是结果与我的第一种方法相同的证明（必须先修复dimnames() 才能使identical() 调用正常工作，但这很简单）：

dimnames(a3) <- dimnames(a2);
identical(a3,a2);
## [1] TRUE

---

性能

以下是使用system.time() 进行的一些基本性能测试，以了解第二种解决方案的优越性：

first <- function() a2 <- aperm(apply(a1,1:2,function(x) tapply(x,strftime(dates,'%m'),mean)),c(2,3,1));
second <- function() a3 <- tapply(a1,list(rep(1:R,C*Z),rep(1:C,each=R,times=Z),rep(strftime(dates,'%m'),each=R*C)),mean);
system.time({ first() });
##    user  system elapsed
##   3.672   0.015   3.719
system.time({ first() });
##    user  system elapsed
##   3.672   0.016   3.720
system.time({ second() });
##    user  system elapsed
##   1.797   0.344   2.135
system.time({ second() });
##    user  system elapsed
##   1.719   0.391   2.124