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

Question

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

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

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

然而，在R中，由于性能原因，“增长”向量被拒绝。

当然，我可以“成片生长”：预先分配一个“大小很好”的向量，填充它，当它满的时候把它的长度增加一倍，最后把它切成大小。但是这感觉很容易出错，并且会导致代码不雅。

有更好的或规范的方法来做这件事吗？(在我的实际应用中，计算和停止标准当然要复杂一些。)

对一些有用的评论的答复

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

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

当向量增长时，是否需要保留所有的值？

是的，我知道。

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

我的实际用例涉及动态计算，我不能简单地将其向量化(否则我不会问)。

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

Answer 1

我可以“成批生长”：预先分配一个“大小好”的向量，填充它，当它满的时候把它的长度增加一倍，最后把它切成大小。但是这感觉很容易出错，并且会导致代码不雅。

听起来你指的是在循环中收集未知数目的结果的公认答案。你把它编好并试过了吗？长度加倍的想法已经足够了(见这个答案的结尾)，因为长度会按几何级数增长。我将在以下中演示我的方法。

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

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