r - 使用 R 替代具有 "dynamic"变量的 for 循环

Question

我是 StackOverflow 的新手，尽管我使用 R 已经有一段时间了。我正在努力解决一个我无法在网站上找到任何答案的问题。如果我的追求不够准确，请纠正我。

我有两个 3d 阵列，在这个简化的情况下为 256x256x200。第一个是一个字段，第二个由索引组成，范围从 1 到 8。我想根据索引的值和计数计算每个垂直级别的平均值，即 200 个级别的字段的平均值对于每个索引(从 1 到 8)。仅当有足够的索引计数(即循环中的 if 条件)时才应执行此操作。我的输出必须是一个 8x200 的矩阵。

例如，我创建了两个随机数组。下面是我正在使用的基本代码:

nz=200
lev=1:nz
indices=8
var0=array(rnorm(256*256*nz),dim=c(256,256,nz))
#octo=array(sample(1:indices),dim=c(256,256,nz)) 
octo=array(sample(1:indices,size=256*256*nz,replace=T),dim=c(256,256,nz))
counts=apply(octo,3,function(x) table(factor(x,levels=1:indices)))
#thr=0.1
thr=0.125
np=length(var0[,1,1])*length(var0[1,,1])
profile=array(NA,dim=c(nz,indices))


t0=proc.time()
for (i in 1:indices)
{
    for (z in 1:length(lev)) 
    {
       if (counts[i,z]/np>thr) 
       {v0=var0[,,z];  profile[z,i]=counts[i,z]/np*mean(v0[octo[,,z]==i],na.rm=T)} 
    }
}
print(proc.time()-t0)

user  system elapsed 
5.169   0.001   5.170 

我尝试使用 apply 系列函数，但我无法以合理有效的方式将其写下来，因为我需要每次计算都考虑到一个改变其级别的“动态”变量(即 octo 和 计数 变量)。我的真实案例是通过更大的矩阵制作的，这应该在几十个领域完成，因此时间非常重要。
您知道任何更快的替代方案吗？
非常感谢您的帮助!

编辑:我更正了 octo 的原始定义，并调整了阈值 thr 。这样 if 条件就有意义了，因为它并不总是受到尊重。

Answer 1

这是一个 data.table避免循环和/或应用语句的 reshape 解决方案:

nz=200
lev=1:nz
indices=8
var0=array(rnorm(256*256*nz),dim=c(256,256,nz))
octo=array(sample(1:indices),dim=c(256,256,nz))
counts=apply(octo,3,function(x) table(factor(x,levels=1:indices)))
thr=0.1
np=length(var0[,1,1])*length(var0[1,,1])
profile=array(NA,dim=c(nz,indices))


# From here load data.table to do the manipulation
# reshape2 to convert back into a matrix at the end
library(data.table)
library(reshape2)

# Take the data long and convert to data.table
var01 <- setDT(melt(var0))
octo1 <- setDT(melt(octo))

# Join the data to get corresponding data
# EDIT, it currently works, but I think that's because all data is defined
# adding nomatch in case of missing data
octo1 <- octo1[var01, on = c('Var1','Var2','Var3'), nomatch = NA] 

# Make our calculation grouping by the vertical dimension and the value
profile <- octo1[,if(.N/np > thr) .N / np * mean(i.value, na.rm = TRUE) else NA, by = .(value,Var3)]

# Recast to matrix
profile <- acast(profile, value ~ Var3, mean, value.var = 'V1')