使用 R 替代带有“动态”变量的 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 和 counts 变量）。我的真实案例是通过更大的矩阵制作的，这应该在几十个领域完成，因此时间非常重要。 您知道任何更快的替代方案吗？ 非常感谢您的帮助！

编辑：我更正了 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')

Answer 2

我认为我通过 sapply 找到了一个很好的解决方案，包括 thr

f1<-function()
{   
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) } }
}
return(profile)
}

f2<-function()
{
profile=sapply(lev, function(i) {
            v0=var0[,,i];
            mV=sapply(1:indices, function(j) {mean(v0[octo[,,i] == j], na.rm = TRUE)})
            counts[,i]/np*mV
    })

profile[counts/np <= thr]=NA
profile<-matrix(profile, nz, indices, byrow = TRUE)
return(profile)
}

f3<-function()
{
profile=sapply(lev, function(i) {
            v0=var0[,,i];
            mV=sapply(1:indices, function(j) {if (counts[j,i]/np>thr) {mean(v0[octo[,,i] == j], na.rm = TRUE)} else {NA}})
            counts[,i]/np*mV
    })

profile<-matrix(profile, nz, indices, byrow = TRUE)
return(profile)
}

其实 f1() 是原版，f2() 是@parksw3 的版本，f3() 是我稍微改进的版本。

benchmark(f1(),f2(),f3(),replications=10)

test   replications elapsed relative user.self sys.self user.child  sys.child
1 f1()           10  27.382    1.411    27.375        0          0         0
2 f2()           10  35.195    1.814    35.186        0          0         0
3 f3()           10  19.403    1.000    19.392        0          0         0

这样，它总是比标准循环快。 data.table 可能更快，但它需要完全更改我目前无法执行的数据结构。希望这会有所帮助！

Answer 3

这在我的机器上似乎更快：

profile2 <- sapply(lev, function(i){
    v0 <- var0[,,i]
    mV <- sapply(1:indices, function(j){
        mean(v0[octo[,,i] == j], na.rm = TRUE)
    })
    counts[,i]/np*mV
})
profile2[counts/np > thr] <- NA
profile2<- t(profile2)

all.equal(profile, profile2)
## TRUE

我尝试将它们与microbenchmark 包进行比较，但需要相当长的时间...这是我与rbenchmark 包进行的快速比较

f1 <- function(){
    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)
            } 
        }
    }
}

f2 <- function(){
    prof <- sapply(lev, function(i){
        v0 <- var0[,,i]
        mV <- sapply(1:indices, function(j){
            mean(v0[octo[,,i] == j], na.rm = TRUE)
        })
        counts[,i]/np*mV
    })
    profile2[counts/np > thr] <- NA
    profile2<- t(profile2)
}

library(rbenchmark)
benchmark(f1(), f2(), replications = 10)

我将两个代码都放入一个函数并进行了测试。结果如下：

##   test replications elapsed relative user.self sys.self
## 1 f1()           10   89.03    1.342     85.15     1.72
## 2 f2()           10   66.34    1.000     61.50     0.75