R while循环与向量条件答案

【问题标题】：R while loop with vector conditionR while循环与向量条件
【发布时间】：2015-04-22 03:51:24
【问题描述】：

我想对使用 while 循环的函数进行矢量化。

原来的功能是

getParamsLeadtime <- function(leadtimeMean, in_tolerance, tolerance){
  searchShape=0
  quantil=0

  # iterates the parameters until the percentage of values is within the interval of tolerance
  while (quantil < in_tolerance){
    searchShape = searchShape+1
    quantil <- pgamma(leadtimeMean+tolerance,shape=searchShape,rate=searchShape/leadtimeMean) -
                  pgamma(leadtimeMean-tolerance,shape=searchShape,rate=searchShape/leadtimeMean) 
  }

  leadtimeShape <- searchShape
  leadtimeRate <- searchShape/leadtimeMean 

  return(c(leadtimeShape, leadtimeRate))
}

我想对该函数进行矢量化调用，以将其应用于数据框。目前我正在循环浏览它：

leadtimes <- data.frame()

for (a in seq(92:103)) {
  leadtimes <- rbind(leadtimes, getParamsLeadtime(a, .85,2))
}

当我尝试对函数进行向量化时，while 似乎不接受向量作为条件。出现以下警告：

Warning message:
In while (input["U"] < rep(tolerance, dim(input)[1])) { :
  the condition has length > 1 and only the first element will be used

这让我假设 while 不喜欢向量。你能告诉我如何对函数进行矢量化吗？

在旁注中，我想知道为什么生成的 leadtimes-data.frame 的列名似乎是值：

> leadtimes
   X1     X1.1
1   1 1.000000
2   1 0.500000
3   4 1.333333
4   8 2.000000
5  13 2.600000
6  19 3.166667
7  25 3.571429
8  33 4.125000
9  42 4.666667
10 52 5.200000
11 63 5.727273
12 74 6.166667

【问题讨论】：

你能提供一些示例数据吗？
如果将函数和 for 循环复制到 R 脚本中，它应该会运行。实际上，示例数据在 for 循环中。
对，你是……对不起！
你不能向量化这个函数。你总是需要一个循环。但是，您应该避免R inferno 的第二个圆圈。例如，您可以使用replicate。如果您深入研究统计数据，您可能会为您的函数找到替代算法，但我不确定。
pgamma 是矢量化的 - 您只需要确保 args 的长度匹配（即 length(shape) 应该等于 length(rate)）。

标签： r function while-loop dataframe vectorization

【解决方案1】：

这是一个非常高效的选项。

对于+tol 和-tol 的情况，在足够大的shp 序列上，我们将pgamma 的计算向量化。我们计算（矢量化）差异，并与in_tol 进行比较。大于in_tol 的向量的第一个元素的索引（减1，因为我们从0 开始）是shp 的最小值，它导致pgamma 大于in_tol。

f <- function(lead, in_tol, tol) {
  shp <- which(!(pgamma(lead + tol, 0:10000, (0:10000)/lead) - 
                 pgamma(lead - tol, 0:10000, (0:10000)/lead)) 
               < in_tol)[1] - 1
  rate <- shp/lead
  c(shp, rate)
}

然后，我们可以sapply 这个范围内的平均交货时间。

t(sapply(1:12, f, 0.85, 2))

##       [,1]     [,2]
##  [1,]    1 1.000000
##  [2,]    1 0.500000
##  [3,]    4 1.333333
##  [4,]    8 2.000000
##  [5,]   13 2.600000
##  [6,]   19 3.166667
##  [7,]   25 3.571429
##  [8,]   33 4.125000
##  [9,]   42 4.666667
## [10,]   52 5.200000
## [11,]   63 5.727273
## [12,]   74 6.166667

system.time(leadtimes <- sapply(1:103, f, 0.85, 2))

##   user  system elapsed 
##   1.28    0.00    1.30

您只需要确保为形状参数选择一个合理的上限（这里我选择了 10000，这非常慷慨）。请注意，如果您不选择足够高的上限，一些返回值将是NA。

【讨论】：