用 apply/maply/etc 替换循环。加快速度

Question

我有以下使用嵌套循环的代码。这使用样本随机数据和样本（相对于实际应用来说较小）数字

N

、

Taumax

和

Tmax

。

N <- 10000
Taumax <- 50
Tmax <- 100

set.seed(42)

X1 <- matrix(rnorm(N * (Tmax+1)), N)
X2 <- matrix(rnorm(N * (Tmax+1)), N)
X3 <- matrix(rnorm(N * (Tmax+1)), N)
Phi <- matrix(rnorm(Taumax * (Tmax+1)), Taumax)
Psi <- matrix(rnorm(Taumax * 3), Taumax)

y <- array(0.0, dim = c(N, Taumax, (Tmax + 1)))
for (t in 1:(Tmax + 1)) {
  for (tau in 1:Taumax) {
    y[, tau, t] <-
      exp(-(
        Phi[tau, t] + Psi[tau, 1] * X1[, t] +
          Psi[tau, 2] * X2[, t] +
          Psi[tau, 3] * X3[, t]
      ) / tau) - 1
  }
}

我很好奇是否可以通过使用例如

apply

、

mapply

或某种其他类型的向量化（或矩阵代数，例如那件事）。我以前用过

apply

，但距离熟练还差很远。

可能会带来麻烦的一件事是循环索引之一 (

tau

) 作为最终计算中的一个因素。

（还有

X1

、

X2

、

X3

的东西，它们一起也可以是一个数组

X

。但是有没有一种方法可以加快速度？）

Answer 1

您可以用

outer()

函数替换内部循环，这应该可以提高速度。

N <- 10000
Taumax <- 50
Tmax <- 100

set.seed(42)

X1 <- matrix(rnorm(N * (Tmax+1)), N)
X2 <- matrix(rnorm(N * (Tmax+1)), N)
X3 <- matrix(rnorm(N * (Tmax+1)), N)
Phi <- matrix(rnorm(Taumax * (Tmax+1)), Taumax)
Psi <- matrix(rnorm(Taumax * 3), Taumax)

y <- array(0.0, dim = c(N, Taumax, (Tmax + 1)))
for (t in 1:(Tmax + 1)) {
   for (tau in 1:Taumax) {
      y[, tau, t] <-
         exp(-(
            Phi[tau, t] + Psi[tau, 1] * X1[, t] +
               Psi[tau, 2] * X2[, t] +
               Psi[tau, 3] * X3[, t]
         ) / tau) - 1
   }
}

#Replace the inner loop with the outer function: %o%
z <- array(0.0, dim = c(N, Taumax, (Tmax + 1)))
for (t in 1:(Tmax + 1)) {
   temp <- Phi[,t] + Psi[, 1] %o% X1[, t] + Psi[, 2] %o% X2[, t] + Psi[, 3] %o% X3[, t]
   z[, , t] <- t(exp(-temp /c(1:Taumax)) -1)
}

#compare results
identical(y, z)