计算事先不知道长度的矢量 - 我应该"增长"吗？

Question

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

一个简单的解决方案是"增长"它:从一个小的或空的向量开始,并连续添加新的条目,直到达到停止标准.例如:

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

然而,出于性能原因,R中的"增长"载体是不受欢迎的.

当然,我可以"以块的形式增长":预先分配一个"大小合适"的矢量,填充它,当它满时加倍它的长度,最后将其缩小到大小.但这感觉容易出错,并且会产生不优雅的代码.

有没有更好或规范的方法来做到这一点？(在我的实际应用中,当然,计算和停止标准有点复杂.)

---

回复一些有用的评论

即使您事先不知道长度,您是否知道它理论上可能具有的最大长度？在这种情况下,我倾向于使用该长度初始化向量,并且在循环切割NA之后或基于最新的索引值移除未使用的条目.

不,事先不知道最大长度.

随着向量的增长,你需要保留所有的值吗？

是的,我愿意.

那么rand_num <- runif(300); rand_num[cumsum(rand_num) < 100]你选择一个足够大的向量,你知道条件满足的概率很高吗？你当然可以检查一下,如果不符合则使用更大的数字.我已经测试过,直到runif(10000)它仍然比"增长"更快.

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

具体来说,为了近似负二项式随机变量的卷积,我需要计算2007年Furman中定理2中整数随机变量$ K $的概率质量,直到高累积概率.这些质量$ 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通过分析来增加 vector 花费了多少时间。

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