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

Question

我正在尝试计算虚构蜥蜴种群规模的灭绝概率。为此，我在 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) 绘制直方图。

提前致谢!

Answer 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