向量的阶乘

Factorial of a vector

作为新手,我尝试定义自己的函数来计算阶乘。我已经设法构建了完美适用于数字的函数。

fact1 = function(x){
    a=1 
    for(i in 1:x){
        a = a*i
    }
    return(a)
}   

factorial = function(x){
    ifelse(x>=0 & round(x) == x , fact1(as.integer(x)),"NA")
}

但是,我怎样才能改进它,使其可以输入一个向量并计算每个元素的阶乘?

使用lapply函数

lapply(c(1,2,3),factorial)
[[1]]
[1] 1

[[2]]
[1] 2

[[3]]
[1] 6

R Documentation for lapply function

添加到上面的 lapply 注释,您还可以使用 vapplysapply 到 return 向量而不是列表:

vapply(c(1, 2, 3),
       factorial, 
       FUN.VALUE = numeric(1))

[1] 1 2 6

这似乎是 Vectorize 的完美案例:只需在 factorial 函数的定义周围使用 Vectorize 即可使其在其输入上向量化。

fact1 = function(x){
  a=1 
  for(i in 1:x){
    a = a*i
  }
  return(a)
}   

factorial = Vectorize(function(x){
  ifelse(x>=0 & round(x) == x , fact1(as.integer(x)),"NA")
})

factorial(c(1,2,3))
#> [1] 1 2 6

问题的答案似乎有点复杂。 阶乘已经是一个存在的函数,如果你有一些数据,你可以简单地将它放入函数中,它就这样被矢量化了。如果要将负数定义为 return 0,也可以使用逻辑语句将其合并。请注意,我使用的是下面的内置函数 factorial 而不是问题中的函数。

dat <- round(runif(1000, -10, 10))
dat_over_zero <- dat > 0 
fact_vector <- numeric(1000)
fact_vector <- factorial(dat[dat_over_zero])

现在,如果您只是创建一个练习来学习,您可以使用相同的想法非常简单地向量化函数,避免不必要的 for 循环。只需使用一个循环并在此循环中迭代向量中的每个元素。

R_factorial <- function(x){
  if(!is.numeric(x) || length(dim(x)))
    stop("X must be a numeric vector!")
  #create an output vector
  output <- numeric(NROW(x))
  #set initial value
  output[x >= 1] <- 1
  output[x < 1] <- NA
  #Find the max factor (using only integer values, not gamma approximations)
  mx <- max(round(x))
  #Increment each output by multiplying the next factor (only on those which needs to be incremented) 
  for(i in seq(2, mx)){
    output[x >= i] <- output[x >= i] * i
  }
  #return output
  output
}

注意几点:

  1. 首先使用 output <- numeric(length) 分配整个向量,其中长度是输出的数量(例如这里的 length(x) 或更一般的 NROW(x))。
  2. 对 none 数值使用 R 常量 NA 而不是 "NA"。第一个被识别为数字,而后者将在字符向量中更改您的向量。

现在备选答案建议 lapply 或 vapply。这或多或少与遍历向量中的每个值并对每个值使用函数相同。因此,向量化函数通常是一种缓慢(但非常可读!)的方法。但是,如果可以避免这种情况,您通常可以获得速度提升。 For loops 和 apply 不一定不好,但与向量化函数相比,它通常要慢得多。请参阅 this Whosebug page,它以非常容易理解的方式解释了原因。 另一种替代方法是使用建议的 Vectorize 函数。这是一个快速而肮脏的解决方案。根据我的经验,它通常比执行一个简单的循环要慢,而且它可能会对多参数函数产生一些意想不到的副作用。它不一定是坏事,因为底层代码的可读性通常会有所提高。


速度比较

现在,与替代答案相比,矢量化版本要快得多。使用 microbenchmark 包中的 microbenchmark 函数,我们可以看到究竟快了多少。下面显示了多少(注意这里我在问题描述中使用阶乘函数):

microbenchmark::microbenchmark(R_factorial = R_factorial(x),
                               Vapply = vapply(x,
                                              factorial, 
                                              FUN.VALUE = numeric(1)),
                               Lapply = lapply(x, factorial),
                               Vfactorial = Vfactorial(x))
Unit: microseconds
        expr       min        lq      mean    median       uq       max neval
 R_factorial   186.525   197.287  232.2394  212.9565  241.464   395.706   100
      Vapply  2209.982  2354.596 3004.9264 2428.7905 3842.265  6165.144   100
      Lapply  2182.041  2299.092 2584.3881 2374.9855 2430.867  5061.852   100
Vfactorial(x) 2381.027 2505.4395 2842.9820 2595.3040 2669.310  5920.094   100

正如你所见,R_factorial 与 vapply 或 lapply (2428.8 / 212.96 = 11.4) 相比大约快 11 - 12 倍。这是相当大的速度提升。可以进行其他改进以进一步加快速度(例如,使用阶乘近似算法、Rcpp 和其他选项),但对于此示例,它可能就足够了。

您还可以使用类型安全 purrr::map_dbl-函数:

purrr::map_dbl(c(1,2,3), fact1)

[1] 1 2 6