R下标超出范围错误-for循环自定义函数答案

【问题标题】：R subscript out of bound error - for loop over custom functionR下标超出范围错误-for循环自定义函数
【发布时间】：2018-03-20 20:04:36
【问题描述】：

我编写了一个函数来模拟遗传漂变，我想在 t（世代数）的各种值上循环它，但是每当我这样做时，我都会收到以下错误：

locifreq<-runif(49, .4, 0.8)

gen<-1:100
for (i in 1:length(gen)){
  pop[i]<-lapply(locifreq,wright.fisher,3000,200,gen[i])
}
Error in `[<-`(`*tmp*`, i, j, value = rbinom(1, 2 * N, prob = k[i - 1,  : 
  subscript out of bounds

我认为这是因为我的函数无法生成适当的矩阵，或者它无法访问列表中的矩阵（尽管我可能完全错了！），但是我不确定如何解决这个问题，因为每次尝试都会导致错误的下标数量。

我的函数代码如下

wright.fisher<-function(p,Ne,nsim,t){
  N <-Ne/2
  NA1 <- 2*N*p 
  NA2 <- 2*N*(1-p)
  k <- matrix(0, nrow = t, ncol = nsim)
  k[1,] <- rep(NA1, nsim) 
  for (j in 1:nsim) {
    for (i in 2:t) {
      k[i, j] <- rbinom(1, 2*N, prob = k[i-1, j] / (2*N))  
    }
  }
  k <- as.matrix(k/(2*N)) 
  t(k)
}

有谁知道如何解决这个问题？

【问题讨论】：

在第一个循环中，i=1，您正在传递 gen[1]=t=1，并且 k 被定义为 1 行和 nsim 列，但您正在尝试为 k 赋值[2,1] 不存在。

标签： r loops for-loop

【解决方案1】：

您的这部分代码重现了您的问题：

wright.fisher(locifreq[1], 3000, 200, 1)

导致下标越界错误的特定步骤是循环的第一步中第 9 行中k[i,j] 的赋值。此时，i 为 2，j 为 1。您通常可以通过将对象分配给超出其范围的索引来扩展对象，但您不能以需要更改维数。（这就是 k[2,1] 在那个时候所做的，因为它是一个 1 行 200 列的矩阵）。

k <- matrix(1:200, 1)
k[2,1] <- 5
# Error in `[<-`(`*tmp*`, 2, 1, value = 5) : subscript out of bounds

您可以通过让函数初始化矩阵以获得正确的维数来解决此问题。对函数的第 5 行进行更改

k <- matrix(0, nrow = max(2, t), ncol = nsim)

我选择了max(2, t)，因为这是嵌套 for 循环的要求：for (i in 2:t)。它完成了，解决了您帖子中的问题，但请确认它会产生您想要的行为。

此外，您的代码中有许多不是特别有效的习语，但这是另一个问题。

如何在 R 中调试

如果你想知道如何调试这样的东西……我正在考虑写一篇关于如何使用浏览器的教程（这里经常出现）。在函数体的顶部添加一行 browser()。然后，当您尝试运行浏览器时，它会在浏览器中将其拉起，这将允许您一次单步执行一个语句。输入help 以查看如何导航。在 RStudio 中，密切关注 Environment 选项卡以查看变量包含的值。请注意，您可以在浏览器中计算任何 R 表达式。这将帮助您检查您认为导致问题的原因。

【讨论】：

这工作得比较好，但不幸的是，当我循环该函数时，它不会遍历每个级别的 locifreq。理想情况下，结果应该给我一个包含 gen 列表数量的列表 - 每个列表包含 49 个矩阵
@MolEcologist29 这是一个不同的问题:)（而不是如何修复下标越界错误）。该问题已解决，但正如我所说，您需要在解决此错误后确认，您的函数为您提供了您想要的行为。您是否尝试过使用浏览器来尝试回答您的新问题？如果您想将新问题作为附加问题发布，我相信我们可以提供帮助。不过，我鼓励您试用浏览器！

【解决方案2】：

您的 Wright-Fisher 模拟看起来有点复杂，但我不知道您在模拟什么，所以也许这是正确的。我将使用一个更简单的版本，并在给定频率和有效种群规模的情况下对下一代的基因频率进行采样。

wf <- function(f, Ne) {
    rbinom(1, 2*Ne, prob = f) / (2*Ne)
}

真正的问题是表以及如何索引到它。如果wf 在给定f 和Ne 的情况下计算一个频率，并且您希望no_sims 模拟no_gen 代，您可以构建一个no_gens x no_sims 表并将其填满，使得行对应于代和列对应于独立的模拟。它可能看起来像这样：

no_sims <- 5
no_gens <- 4
Ne <- 10000

locifreq <- runif(5, .4, 0.8)
sims <- matrix(NA, ncol = no_sims, nrow = no_gens)
sims[1,] <- runif(5, .4, 0.8) # first generation
for (gen in 2:no_gens) {
    sims[gen,] <- sapply(sims[gen-1,], FUN = wf, Ne = Ne)
}

这些（小）参数的结果如下所示：

> sims
          [,1]      [,2]      [,3]      [,4]      [,5]
[1,] 0.5948922 0.7185469 0.6290239 0.4303951 0.5701607
[2,] 0.5987000 0.7108500 0.6254000 0.4270500 0.5754000
[3,] 0.6020000 0.7103500 0.6260000 0.4320500 0.5723500
[4,] 0.5982500 0.7110500 0.6276000 0.4332500 0.5751000

【讨论】：