This page is a method.
The page identifier is Op_en7816 |
|---|
| Moderator:Jouni (see all) |
|
Give your opinion to the peer rating of the content of this page. |
| Upload data
|
How to model exponential growth in a population?
library(OpasnetUtils) objects.latest(Op_enXXXX) # growth2 growth2(100, 30, 0.1, 0.05, 10) |
# This is code Op_enXXXX/initiate on page [[Population growth]]
library(OpasnetUtils)
growth2 <-function(
N, # size of population in the beginning
tmax, # end of follow-up time
r, # growth rate (fraction per time unit)
sdr, # standard deviation of r
it # number of iterations
){
popsize <- data.frame(time=c(0:tmax)) # the data.frame for the final data with tmax rows
for(i in 1:it){ # doing all this it times
pop.size <- c(N, rep(NA, tmax))
for(t in 1:tmax){ # this is the same bit as before
nt <- ifelse(t==1, N, nt1)
nt1 <- nt*exp(rnorm(1, mean=r, sd=sdr))
pop.size[t+1] <- nt1
}
popsize <- cbind(popsize, pop.size) # adding the data of each iteration to the same data frame
}
popsize$mean.size <- rowMeans(popsize[,c(2:it)]) # mean sizes for times from iterations
quantiles <- data.frame()
# with a loop calculating quantiles for each row, aka. each time.
for(p in 1:tmax+1){
quantiles <- rbind(quantiles, quantile(popsize[p,c(2:it+1)], probs=c(0.025,0.975)))
}
quantiles <- stack(quantiles) #flipping the quantile
quantiles[c((tmax+1):(2*tmax)),"values"] <- rev(quantiles[c((tmax+1):(2*tmax)),"values"])
quantiles$time <- c(1:tmax, tmax:1) # adding the time as x-axis points
plot(x=popsize$time, y=popsize$mean.size, type="n",
xlab = "Time", ylab="Population size", main="Population growth over time",
ylim = c(0, max(quantiles$value))) #y lim so that the upper quantile is also visible
polygon(quantiles$time, quantiles$values, col="palegreen") # quantile lines
lines(popsize$time, popsize$mean.size) # pop size line on top
}
objects.store(growth2)
cat("Function growth2 stored.\n")
|
