Population growth
Jump to navigation
Jump to search
Moderator:Jouni (see all) |
|
Upload data
|
Question
How to model exponential growth in a population?
Moderator:Jouni (see all) |
|
Upload data
|
How to model exponential growth in a population?
library(OpasnetUtils) objects.latest("Op_en7816", code_name = "initiate") # growth2 growth2(100, 30, 0.1, 0.05, 10) |
# This is code Op_en7816/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") |