修改后的引导答案

【问题标题】：Modified Bootstrapping修改后的引导
【发布时间】：2011-04-01 15:36:29
【问题描述】：

我有兴趣开发一种修改后的引导程序，该引导程序可以对长度为 x 的某个向量进行替换，但在停止采样之前必须满足许多标准。我试图计算种群增长率的 lambda 的置信区间，10000 次迭代，但在某些个体分组中，比如向量 13，很少有个体从组中生长出来。典型的自举会导致相当数量的实例，其中该向量不会发生增长，因此模型会崩溃。每个向量由一定数量的 1、2 和 3 组成，其中 1 表示留在一个组中，2 表示从一个组中生长出来，3 表示死亡。这是我到目前为止没有修改的内容，从时间上讲，这可能不是最好的方法，但我是 R 新手。

st13 <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,  
          1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,3,3)
#runs
n <- 10000
stage <- st13
stagestay <- vector()
stagemoved <- vector()
stagedead <- vector()
for(i in 1:n){
      index <- sample(stage, replace=T)
      stay <- ((length(index[index==1]))/(length(index)))
      moved <- ((length(index[index==2]))/(length(index)))
      stagestay <- rbind(stagestay,stay)
      stagemoved <- rbind(stagemoved,moved)
}

目前，此示例那么我的问题是：我可以通过什么方式修改示例函数以继续对这些数字进行采样，直到“索引”的长度至少与 st13 相同，并且直到“索引”中存在至少 1 个 2 的实例？

非常感谢，克里斯托弗·亨尼格硕士生密西西比大学牛津，MS，38677

【问题讨论】：

我有点困惑。这是您要修改的行：index <- sample(stage, replace=T)？如果是这样，这将始终产生一个向量，只要st13。那么真正的问题是更多关于如何使sample() 返回一个包含 1 个或多个值 > 1 的向量，然后再进行计算？

标签： r matrix modeling sample statistics-bootstrap

【解决方案1】：

更新： @lselzer 的回答提醒我，要求样本的长度至少与st13 一样长。我上面的代码只是不断采样，直到找到包含2 的引导示例。 @lselzer 的代码会增加样本，一次增加 1 个新索引，直到样本包含 2。这是非常低效的，因为您可能不得不多次调用sample()，直到您得到2。在示例中返回 2 之前，我的代码可能会重复很长时间。那么我们还能做得更好吗？

一种方法是使用对sample() 的一次调用来抽样一个带有替换的大样本。检查哪些是2s，并查看第一个length(st13) 条目中是否有2。如果有，则返回这些条目，如果没有，则在大样本中找到第一个 2 并返回所有条目，直到包含那个条目。如果没有2s，添加另一个大样本并重复。这是一些代码：

#runs
n <- 100 #00
stage <- st13
stagedead <- stagemoved <- stagestay <- Size <- vector()
sampSize <- 100 * (len <- length(stage)) ## sample size to try
for(i in seq_len(n)){
    ## take a large sample
    samp <- sample(stage, size = sampSize, replace = TRUE)
    ## check if there are any `2`s and which they are
    ## and if no 2s expand the sample
    while(length((twos <- which(samp == 2))) < 1) {
        samp <- c(samp, sample(stage, size = sampSize, replace = TRUE))
    }
    ## now we have a sample containing at least one 2
    ## so set index to the required set of elements
    if((min.two <- min(twos)) <= len) {
        index <- samp[seq_len(len)]
    } else {
        index <- samp[seq_len(min.two)]
    }
    stay <- length(index[index==1]) / length(index)
    moved <- length(index[index==2]) / length(index)
    stagestay[i] <- stay
    stagemoved[i] <- moved
    Size[i] <- length(index)
}

这是一个非常退化的向量，46 个条目中只有一个 2：

R> st14 <- sample(c(rep(1, 45), 2))
R> st14
 [1] 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[39] 1 1 1 1 1 1 1 1

如果我在其上使用上述循环而不是 st13，则在 100 次运行中的每一次运行中获得 2 所需的最小样本量如下：

R> Size
  [1]  65  46  46  46  75  46  46  57  46 106  46  46  46  66  46  46  46  46
 [19]  46  46  46  46  46 279  52  46  63  70  46  46  90 107  46  46  46  87
 [37] 130  46  46  46  46  46  46  60  46 167  46  46  46  71  77  46  46  84
 [55]  58  90 112  52  46  53  85  46  59 302 108  46  46  46  46  46 174  46
 [73] 165 103  46 110  46  80  46 166  46  46  46  65  46  46  46 286  71  46
 [91] 131  61  46  46 141  46  46  53  47  83

所以这表明我选择的sampSize (100 * length(stage)) 在这里有点矫枉过正，但由于我们使用的所有运算符都是矢量化的，我们可能不会因为初始样本过长而受到太多惩罚大小，我们当然不会招致任何额外的sample() 调用。

原文： 如果我理解正确，问题是sample() 可能根本不会返回任何2 指标。如果是这样，我们可以使用repeat 控制流构造继续采样直到它完成为止。

我已经相应地更改了您的代码，并对其进行了一些优化，因为您永远不会像以前那样在循环中增长对象。还有其他可以改进的方法，但我现在会坚持使用循环。解释如下。

st13 <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,  
          1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,3,3)
#runs
n <- 10000
stage <- st13
stagedead <- stagemoved <- stagestay <- vector()
for(i in seq_len(n)){
    repeat {
        index <- sample(stage, replace = TRUE)
        if(any(index == 2)) {
            break
        }
    }
    stay <- length(index[index==1]) / length(index)
    moved <- length(index[index==2]) / length(index)
    stagestay[i] <- stay
    stagemoved[i] <- moved
}

这是与您的 Q 相关的主要变化：

    repeat {
        index <- sample(stage, replace = TRUE)
        if(any(index == 2)) {
            break
        }
    }

这样做是重复包含在大括号中的代码，直到触发break 以将我们跳出repeat 循环。所以会发生什么是我们采取引导样本，然后检查任何样本是否包含索引2。如果有任何2s，那么我们就会中断并继续当前 for 循环迭代的其余部分。如果样本不包含任何2s，则不会触发中断，我们会再次循环获取另一个样本。这将一直发生，直到我们得到一个包含2 的样本。

【讨论】：

@lselzer 习惯 - 当然这里没有必要，但我之前用all.equal() 做了一些事情，你需要在isTRUE() 中换行，所以我认为手指会自动驾驶。将删除，谢谢。

【解决方案2】：

对于初学者，sample 有一个 size 参数，您可以使用它来匹配 st13 的长度。您问题的第二部分可以使用while 循环来解决。

st13 <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,  
          1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,2,2,3,3)
    #runs
    n <- 10000
    stage <- st13
    stagestay <- vector()
    stagemoved <- vector()
    stagedead <- vector()
    for(i in 1:n){
          index <- sample(stage, length(stage), replace=T)
          while(!any(index == 2)) {
            index <- c(index, sample(stage, 1, replace = T))
          }
          stay <- ((length(index[index==1]))/(length(index)))
          moved <- ((length(index[index==2]))/(length(index)))
          stagestay[i] <- stay
          stagemoved[i] <- moved
    }

在我写这篇文章时，Gavin 发布了他的答案，这与我的相似，但我添加了 size 参数以确保索引至少具有 st13 的长度

【讨论】：

size 位不需要。 ?sample有For ‘sample’ the default for ‘size’ is the number of items inferred from the first argument, so that ‘sample(x)’ generates a random permutation of the elements of ‘x’ (or ‘1:x’).所以我们根本不需要设置size，它是从stage的长度推断出来的。
@Gavin 是的，但是 OP 说索引应该至少是 length(st13)，并且他想知道如何继续采样，直到找到 2。所以我推断索引可以大于st13，但不能更小。现在我更仔细地阅读了你的代码，我发现它在每个重复循环中都替换了索引，所以每次都是length(index) == length(st13)。
这是正确的。索引可以大于长度（st13），并且在这些没有出现 2 的情况下是必要的，因为报废样本并重新洗牌直到我收到 2 会使我正在计算的值过度膨胀。感谢帮助和 cmets。很有帮助。
@lselzer 是的，我在写评论时忘记了这一点，因此我将其删除。我刚刚在我的代码中添加了一个更新，它显示了一种混合方法，我们采用比所需更大的样本并且只扩展它我们没有得到2s。这应该比每次扩展 1 个单位更有效，因为您可能会多次调用 sample()，直到您收到 2。
@user687814 这在我的代码中不会发生 - 我只是扔掉样本并继续这样做，直到我们得到一个包含 2 的样本。对于您显示的向量，我们几乎总是得到 2，但对于 2 较少的情况，代码可能会循环一段时间，直到找到 2。请参阅我更新的答案，了解与@lselzer 类似的代码，但没有可能多次调用sample() 以扩展样本大小。