Claudia Takahashi Claudia Takahashi - 18 days ago 8
R Question

Store last variable in a for loop equation into an array

I'm trying to calculate the probability of extinction of a fictional lizard population size. To do this, I am running a for loop for 100 simulations over a period of 30 years, and seeing the probability of each simulation from going extinct. At the end of my 100 simulations, I need to plot a histogram depicting the final population size at the end of the 30 year interval. I figured that the easiest way to plot the histogram would be to create a different vector, and store the final population size of each simulation into this vector (

pop
). However, I have no idea how to code for this and have not found an answer online for my predicament.

I am using the following code:

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)
}


So in the above code,
pop
is the vector where I'm storing the population at
N[30]
. The idea is then to use
hist(pop)
to plot the histogram.

Thanks in advance!

Answer

You can get the results in a matrix like this:

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
Comments