我如何计算3D阵列切片的表达式？

Question

我已经为24x24方阵（a）的每24个值计算了这个hc指数。它返回一个包含24个值的向量（hc）。

hc<-vector("numeric",24L)
for (i in 1:24) {
  hc[i]<- -sum((c(a[,i])/colSums(a)[i])*log((c(a[,i])/colSums(a)[i]), base = 2))
  }

我想为数组的1000个矩阵中的每一个计算这个，并且不知道如何继续（函数？，另一个嵌套的“for”语句？...）。并获得1000,24个长度的向量。一个用于阵列中的每个矩阵。我真的很感激任何帮助。谢谢！

Answer 1

如果您确实有一个数组，这是另一种方法：

# Simplified version of your function
f <- function(x) -sum(x / sum(x) * log2(x / sum(x)))

set.seed(1)
n <- 1000

matrices <- replicate(n, matrix(runif(24 ^ 2), nrow = 24))
str(matrices)
#>  num [1:24, 1:24, 1:1000] 0.266 0.372 0.573 0.908 0.202 ...

result <- apply(matrices, c(2, 3), f)
str(result)
#>  num [1:24, 1:1000] 4.36 4.36 4.37 4.36 4.34 ...

如果您的矩阵在列表中：

matrix_list <- lapply(seq_len(n), function(i) matrices[, , i])
list_result <- lapply(matrix_list, apply, 2, f)
str(head(list_result))
#> List of 6
#>  $ : num [1:24] 4.36 4.36 4.37 4.36 4.34 ...
#>  $ : num [1:24] 4.43 4.32 4.31 4.4 4.37 ...
#>  $ : num [1:24] 4.26 4.21 4.31 4.24 4.24 ...
#>  $ : num [1:24] 4.31 4.36 4.27 4.32 4.35 ...
#>  $ : num [1:24] 4.39 4.27 4.35 4.29 4.22 ...
#>  $ : num [1:24] 4.25 4.42 4.19 4.32 4.33 ...

由reprex package创建于2018-03-05（v0.2.0）。

Answer 2

首先，您的代码效率很低。您在每次迭代中计算矩阵中每列的总和两次。最好在开始时计算所有这些，或者只是在需要时计算总和。

以下是更像“R-like”的代码也应该更快：

mat <- matrix(runif(100, 1, 100), nrow = 24, ncol = 24)


mat_res <- sapply(seq_len(ncol(mat)), 
    function(i, mat) {
        col <- mat[, i]
        s <- sum(col)
        sc <- col / s
        l <- log2(col) / s
        - sum(sc * l)
    }, 
    mat = mat)

现在，将其应用于所有1000个矩阵？首先，列出1000个矩阵：

list_of_mats <- lapply(1:1000, 
    function(i) matrix(runif(100, 1, 100), nrow = 24, ncol = 24))

然后，在每个矩阵上应用一个函数：

lapply(list_of_mats, 
    function(mat) {
        ## Note that this process is the same as above!
        sapply(seq_len(ncol(mat)), 
            function(i, mat) {
                col <- mat[, i]
                s <- sum(col)
                sc <- col / s
                l <- log2(col) / s
                - sum(sc * l)
            }, 
            mat = mat)
    })