Claudia Takahashi - 7 months ago 39

R Question

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`

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`

`N[30]`

`hist(pop)`

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
```

Source (Stackoverflow)