将 for 循环方程中的最后一个变量存储到数组中答案

【问题标题】：Store last variable in a for loop equation into an array将 for 循环方程中的最后一个变量存储到数组中
【发布时间】：2016-11-14 18:51:23
【问题描述】：

我正在尝试计算虚构的蜥蜴种群规模灭绝的概率。为此，我在 30 年的时间里运行 100 次模拟的 for 循环，并查看每个模拟灭绝的概率。在我的 100 次模拟结束时，我需要绘制一个直方图，描绘 30 年间隔结束时的最终人口规模。我认为绘制直方图的最简单方法是创建一个不同的向量，并将每个模拟的最终种群大小存储到这个向量中 (pop)。但是，我不知道如何为此编写代码，也没有在网上找到解决我的困境的答案。

我正在使用以下代码：

tmax <- 31
runmax <- 100
Year <- 0:(tmax-1)
N <- numeric(tmax) %vector for the population size
N <- N + 1
epsilon <- numeric(tmax) 
rmax <- 0.87992 %maximum growth rate (a value previously calculated)
K <- 34.64252 %carrying capacity (a value previously calculated)
N[1] <- K
extinct <- 0

for(t in 2:tmax){
    sdr <- 0.9469428
    epsilon[t-1] <- rnorm(1,0,sdr) %this takes into account the random population stochasticity (random chance a population will go extinct)
    N[t] <- exp(rmax*(1-(N[t-1]/K))+epsilon[t-1])*N[t-1]
    if(N[t] < 1.0) {
        N[t] <- 0.0;break
    }
    pop=numeric(runmax)
    pop[1]=N[30]
}

extinct <- extinct + ifelse(N[tmax]<=1,1,0)

plot(Year,N,type='l',ylim=c(0,200))

for(i in 1:runmax){
    N <- numeric(tmax)
    N <- N+1
    N[1] <- K
    for(t in 2:tmax){
        sdr <- 0.9469428
        epsilon[t-1] <- rnorm(1,0,sdr) 
        N[t] <- exp(rmax*(1-(N[t-1]/K))+epsilon[t-1])*N[t-1]
        if(N[t] < 1.0) {
            N[t] <- 0.0
            break
        }
    for(w in 2:runmax){
        pop[w]<- N[30]
    }
}

extinct <- extinct + ifelse(N[tmax]<=1,1,0)
lines(Year,N,col=i)
}

所以在上面的代码中，pop 是我将人口存储在N[30] 的向量。然后的想法是使用hist(pop) 绘制直方图。

提前致谢！

【问题讨论】：

您可能应该添加一些上下文。现在我看到很多变量，我认为rmax和K和N和epsilon的概率密度
值得注意的是hist(pop)是蜥蜴灭绝时发出的声音。
在第一个循环中，您在每次循环中都写入pop -- 为什么？在第二个中，不清楚你在pop 变量中实际保存了什么，但它似乎每次都在重写结果（这里，将pop 的所有值设置为N[30]通过i 上的循环的时间。为什么不设置pop[i] <- N[30]（没有w 上的循环）？
@KenS。我回去并添加了一些编辑。这是否足够的上下文，或者你还需要什么？
@MarkPeterson 我不是 100% 确定我在做什么 - 我意识到每次循环通过时我都在写结果，但我不知道如何不这样做。您使用 pop[i]

标签： r algorithm

【解决方案1】：

你可以得到这样的矩阵结果：

pop=matrix(rep(0,runmax*tmax),ncol=tmax)
for(i in 1:runmax){
  N <- numeric(tmax)
  N <- N+1 # this can be removed
  N[1] <- K
  for(t in 2:tmax){
    sdr <- 0.9469428 # this could be placed outside the loops
    epsilon[t-1] <- rnorm(1,0,sdr)
    N[t] <- exp(rmax*(1-(N[t-1]/K))+epsilon[t-1])*N[t-1]
    if(N[t] < 1.0) {N[t] <- 0.0}
    pop[i,t]=N[t]
    if(N[t] ==0) {break}
  }

  extinct <- extinct + ifelse(N[tmax]<=1,1,0)
  lines(Year,N,col=i)
}
hist(pop[,tmax])  #simulation results for tmax

【讨论】：