每月意味着申请多维数组

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

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

解决方案 1

您可以做的是单独计算一个日期向量，该向量标识与数组的每个 z-index 对应的日期。在 apply() call for each z-line, you can then call strftime() to extract the months for each such date, and group by that month value using tapply() to take monthly mean() 秒内。这是如何完成的：

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 天，但显然你应该确定使用您的真实数据实际涵盖的任何日期。
如您所见，以 1:2 作为边距调用 apply() 是正确的，因为这允许您在每条 z 线上独立操作，这样您就可以按月对该 z 线进行分组，并计算沿该 z 线的每个月的平均值。
不幸的是，apply() 有一个恼人的习惯，即以与人们通常预期不同的换位方式返回结果。对于二维用法，这通常通过简单调用 t(), but since we're working in three dimensions here, we need to call aperm() 来固定维度顺序来解决。
由于我选择的日期是从一月开始，然后按日历顺序推进到月份，结果中的方法最终将按日历月排序。 IOW，a2 中的 z-indexes 1:12 对应于 1 月至 12 月。如果您的日期不是从一月开始，那么这个解决方案应该仍然有效，但您必须注意结果中 z-index 和月份之间的对应关系。例如，我的 "proof of correctness" 代码假定索引 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() have to be fixed first to get the 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

每月意味着申请多维数组

monthly means with apply for multidimensional arrays

aggregate

r

apply

解决方案 1

解决方案 2

性能