计算事先不知道长度的向量 - 我应该“增长”它吗？

Question

我需要计算一个向量的条目我事先不知道其长度。如何高效地做到这一点？

一个简单的解决方案是“增长”它：从一个小向量或空向量开始，然后连续添加新条目，直到达到停止标准。例如：

foo <- numeric(0)
while ( sum(foo) < 100 ) foo <- c(foo,runif(1))
length(foo)
# 195

但是，出于性能原因，R 中不赞成“增长”向量。

当然，我可以“分块增长”：预先分配一个“大小合适”的向量，填充它，当它满时将其长度加倍，最后将其缩小。但这感觉很容易出错，并且会导致代码不优雅。

有没有更好或规范的方法来做到这一点？ （在我的实际应用中，计算和停止条件当然要复杂一些。）

---

回复一些有用的cmets

即使你事先不知道长度，你知道它理论上可以有的最大长度吗？在这种情况下，我倾向于使用该长度初始化向量，然后在循环剪切 NA 或根据最新的索引值删除未使用的条目之后。

不，最大长度事先不知道。

随着向量的增长，您是否需要保留所有值？

是的，我愿意。

像rand_num <- runif(300); rand_num[cumsum(rand_num) < 100] 这样的东西怎么样？你选择一个足够大的向量，你知道这个条件很可能会满足？您当然可以检查它，如果不满足，则使用更大的数字。我已经测试到runif(10000) 它仍然比“增长”更快。

我的实际用例涉及动态计算，我不能简单地对其进行矢量化（否则我不会问）。

具体来说，为了近似负二项式随机变量的卷积，我需要计算Furman, 2007中定理2中整数随机变量$K$的概率质量，直到达到高累积概率。这些质量 $pr_k$ 涉及一些复杂的递归和。

Answer 1

我可以“分块增长”：预先分配一个“大小合适”的向量，填充它，当它满时将其长度加倍，最后将其缩小。但这感觉很容易出错，并且会导致代码不优雅。

听起来您指的是Collecting an unknown number of results in a loop 的公认答案。你有没有编码并尝试过？长度加倍的想法绰绰有余（见这个答案的结尾），因为长度会几何增长。我将在下面演示我的方法。

---

出于测试目的，请将您的代码包装在一个函数中。请注意我如何避免为每个 while 测试执行 sum(z)。

ref <- function (stop_sum, timing = TRUE) {
  set.seed(0)                            ## fix a seed to compare performance
  if (timing) t1 <- proc.time()[[3]]
  z <- numeric(0)
  sum_z <- 0
  while ( sum_z < stop_sum ) {
    z_i <- runif(1)
    z <- c(z, z_i)
    sum_z <- sum_z + z_i
    }
  if (timing) {
    t2 <- proc.time()[[3]]
    return(t2 - t1)                      ## return execution time
    } else {
    return(z)                            ## return result
    }
  }

分块对于降低串联的运营成本是必要的。

template <- function (chunk_size, stop_sum, timing = TRUE) {
  set.seed(0)                            ## fix a seed to compare performance
  if (timing) t1 <- proc.time()[[3]]
  z <- vector("list")                    ## store all segments in a list
  sum_z <- 0                             ## cumulative sum
  while ( sum_z < stop_sum ) {
    segmt <- numeric(chunk_size)         ## initialize a segment
    i <- 1
    while (i <= chunk_size) {
      z_i <- runif(1)                    ## call a function & get a value
      sum_z <- sum_z + z_i               ## update cumulative sum
      segmt[i] <- z_i                    ## fill in the segment
      if (sum_z >= stop_sum) break       ## ready to break at any time
      i <- i + 1
      }
    ## grow the list
    if (sum_z < stop_sum) z <- c(z, list(segmt))
    else z <- c(z, list(segmt[1:i]))
    }
  if (timing) {
    t2 <- proc.time()[[3]]
    return(t2 - t1)                      ## return execution time
    } else {
    return(unlist(z))                    ## return result
    }
  }

让我们先检查一下正确性。

z <- ref(1e+4, FALSE)
z1 <- template(5, 1e+4, FALSE)
z2 <- template(1000, 1e+4, FALSE)

range(z - z1)
#[1] 0 0

range(z - z2)
#[1] 0 0

接下来我们比较一下速度。

## reference implementation
t0 <- ref(1e+4, TRUE)

## unrolling implementation
trial_chunk_size <- seq(5, 1000, by = 5)
tm <- sapply(trial_chunk_size, template, stop_sum = 1e+4, timing = TRUE)

## visualize timing statistics
plot(trial_chunk_size, tm, type = "l", ylim = c(0, t0), col = 2, bty = "l")
abline(h = t0, lwd = 2)

看起来chunk_size = 200足够好，加速因子是

t0 / tm[trial_chunk_size == 200]
#[1] 16.90598

最后让我们看看用c，通过剖析来增长向量花费了多少时间。

Rprof("a.out")
z0 <- ref(1e+4, FALSE)
Rprof(NULL)
summaryRprof("a.out")$by.self
#        self.time self.pct total.time total.pct
#"c"          1.68    90.32       1.68     90.32
#"runif"      0.12     6.45       0.12      6.45
#"ref"        0.06     3.23       1.86    100.00

Rprof("b.out")
z1 <- template(200, 1e+4, FALSE)
Rprof(NULL)
summaryRprof("b.out")$by.self
#        self.time self.pct total.time total.pct
#"runif"      0.10    83.33       0.10     83.33
#"c"          0.02    16.67       0.02     16.67

---

自适应chunk_size 线性增长

ref 具有O(N * N) 操作复杂度，其中N 是最终向量的长度。 template 原则上具有O(M * M) 复杂性，其中M = N / chunk_size。为了达到线性复杂度O(N)，chunk_size 需要与N 一起增长，但线性增长就足够了：chunk_size <- chunk_size + 1。

template1 <- function (chunk_size, stop_sum, timing = TRUE) {
  set.seed(0)                            ## fix a seed to compare performance
  if (timing) t1 <- proc.time()[[3]]
  z <- vector("list")                    ## store all segments in a list
  sum_z <- 0                             ## cumulative sum
  while ( sum_z < stop_sum ) {
    segmt <- numeric(chunk_size)         ## initialize a segment
    i <- 1
    while (i <= chunk_size) {
      z_i <- runif(1)                    ## call a function & get a value
      sum_z <- sum_z + z_i               ## update cumulative sum
      segmt[i] <- z_i                    ## fill in the segment
      if (sum_z >= stop_sum) break       ## ready to break at any time
      i <- i + 1
      }
    ## grow the list
    if (sum_z < stop_sum) z <- c(z, list(segmt))
    else z <- c(z, list(segmt[1:i]))
    ## increase chunk_size
    chunk_size <- chunk_size + 1
    }
  ## remove this line if you want
  cat(sprintf("final chunk size = %d\n", chunk_size))
  if (timing) {
    t2 <- proc.time()[[3]]
    return(t2 - t1)                      ## return execution time
    } else {
    return(unlist(z))                    ## return result
    }
  }

快速测试验证我们已达到线性复杂度。

template1(200, 1e+4)
#final chunk size = 283
#[1] 0.103

template1(200, 1e+5)
#final chunk size = 664
#[1] 1.076

template1(200, 1e+6)
#final chunk size = 2012
#[1] 10.848

template1(200, 1e+7)
#final chunk size = 6330
#[1] 108.183