StochasticPopulationGrowth.R

#Stochastic geometric population growth rate

#Simulation settings
pgr = 1.05
var.pgr = 0.1
startPop = 10
nGen = 100
ntrials = 1000
pseudoExtinction = 1

####################################################################
#If you are unfamiliar with R, do not edit anything below this line!
####################################################################

#First randomly generate some lambda values
lambdas<-matrix(rnorm(ntrials*nGen, mean = pgr, sd = sqrt(var.pgr)),ncol=ntrials,nrow=nGen)

#Use a histogram to see what they look like
hist(lambdas,col="grey",main="")

#Now run the simulations to see what the resulting population growth looks like
trials = matrix(data = NA, nrow = nGen, ncol = ntrials)
for (j in 1:ntrials){
  popSize = startPop  
  for (i in 2: nGen){
  stoch.pgr = lambdas[i,j]
  popSize = append(popSize, popSize[i-1]*stoch.pgr)
}
trials[,j] = popSize
rm(popSize)
}

#Make a plot of the population trajectories
plot(1:nGen,log(seq(0.1,max(trials),length.out=nGen)),type = "n",axes=F,xlab ="Time", ylab = "N",ylim=log(c(0.1,100000)))
matlines(log(trials),col = "#FF234520",lty=1,lwd=3)
axis(1)
axis(2,at = log(c(0.1,1,10,100,1000,10000,100000)),
     label = c(0.1,1,10,100,1000,10000,100000))
abline(h=log(pseudoExtinction),lty=2)

#Calculate probability of (pseudo)extinction
minvals <- apply(trials,2,min)
nExtinct <- length(minvals[minvals<=pseudoExtinction])
nExtinct/ntrials