我需要减少 R 中的 for 循环吗?如何减少?

Do I need to reduce for-loop in R and how?

我是 R 编码的新手,但我听说 R 中的循环比 Python 或 C 等其他语言慢得多。那么在 R 中编码时是否需要减少循环?

具体来说,在这个模拟代码中,我该如何提高自己的编码水平?

library(moments)
n <- c(5:20)
m <- c(1:10000)
skew <- c()
kurt <- c()
for(num in n){
  beta1 <- c()
  beta2 <- c()
  for(i in m){
    set.seed(num * 10000 + i)
    x <- rnorm(num, mean = 0, sd = 1)
    beta1 <- c(beta1, skewness(x))
    beta2 <- c(beta2, kurtosis(x) - 3)
  }
  skew <- c(skew, quantile(beta1, probs = c(0, 0.01, 0.1, 0.2, 0.5, 0.8, 0.9, 0.99, 1)))
  kurt <- c(kurt, quantile(beta2, probs = c(0, 0.01, 0.1, 0.2, 0.5, 0.8, 0.9, 0.99, 1)))
}

在 R 中不使用 for 循环的一个主要优点是利用其向量化。因此,在 Python 或 C 等语言中,您可以为向量的每个元素编写向量计算代码,而在 R 中,您可以方便地一次为整个向量编写计算代码(请参阅下面的编辑),还可以通过实际使用来减少计算时间快速底层 C、Fortran 等函数。

我会把你想为单个样本量进行的所有计算放入函数 statFUN 并将其放入 lapply 以遍历样本量向量 n.

对于分位数,我们可以使用 apply 或我推荐的 matrixStats::rowQuantiles,因为它更快。

set.seed() 在 运行 和 lapply 之前只需要一次,所有 res 结果都可以用那个种子重现。

n <- 5:20  ## different sample sizes
m <- 1e4   ## number of replications in each iteration
probs <- c(0, 0.01, 0.1, 0.2, 0.5, 0.8, 0.9, 0.99, 1)

library(moments)
library(matrixStats)

statFUN <- function(i, num) {
  r <- replicate(i, {
    x <- rnorm(num, mean=0, sd=1)
    c(kurt=kurtosis(x) - 3, skew=skewness(x))
  })
  # t(apply(r, 1, quantile, probs=probs))  ## using base R
  rowQuantiles(r, probs=probs)  ## using matrixStats
}

set.seed(42)
res <- lapply(n, statFUN, m)

结果

result 是每个样本大小的峰度分位数和偏度分位数的列表。

res
# [[1]]
#               0%          1%         10%         20%
# kurt -0.04710729 -0.04658709 -0.04190536 -0.03670343
# skew -0.03045563 -0.02969417 -0.02284104 -0.01522645
#              50%          80%           90%         99%
# kurt -0.03388803 -0.006250622  1.068998e-03 0.007656657
# skew -0.01028591 -0.006132523 -5.883157e-05 0.005407491
#             100%
# kurt 0.008388619
# skew 0.006014860
# 
# [[2]]
#               0%          1%         10%         20%
# kurt -0.09089922 -0.08859363 -0.06784329 -0.04478737
# skew -0.03252828 -0.03165837 -0.02382918 -0.01513009
#               50%          80%        90%        99%
# kurt -0.023634727 -0.005277533 0.01038904 0.02448896
# skew  0.003433589  0.017711708 0.01947178 0.02105585
#            100%
# kurt 0.02605562
# skew 0.02123186
#
# [...]

哪里

length(res)
# [1] 16

编辑

这里有一个小例子可以更好地说明 R 中向量化的实际含义。虽然在大多数编程语言中,两个向量的加法是按元素编码的,但在 R 中,向量的加法可以直接编码(即在矢量化的方式)。

a <- 1:9
b <- rev(a)

## element wise addition of vectors a and b
s1 <- c()
for (i in seq(a)) {
  s1[i] <- a[i] + b[i]
}
s1
# [1] 10 10 10 10 10 10 10 10 10

## direct addition of vectors a and b (i.e. vectorized)
s2 <- a + b
s2
# [1] 10 10 10 10 10 10 10 10 10

我们可以研究 *apply 系列而不是 for 循环。但是,大多数情况下仍然隐藏着 for 循环。 (要查看功能代码类型,例如 lapply 不带括号或任何内容。)

您可能想阅读例如那些很棒的问答:

  • Grouping functions (tapply, by, aggregate) and the *apply family
  • Why are loops slow in R?
  • How to see the source code of R .Internal or .Primitive function?

注:向量化其实只是R的语言特性,所谓的“向量化函数”内部往往使用C、Fortran等代码,在其中你还是在最后找到 for 循环,但是使用一种更快的语言。例如,当我们使用 sum().

时调用 source code of summary.c