# Stochastic model of density dependent population growth
# with variation in R0 & K possible
# Code also allows plotting of counts that contain
# observation errror, i.e., estimates
N0 <- 50
T <- 15
muR0 <- 0.2
sdR0 <- 0.0
muK <- 100
sdK <- 0
sdObs <- 0
N <- matrix(NA, T + 1, 1)
N[1] <- N0
Nobs <- matrix(NA, T + 1, 1)
# add observation error to 1st count
# and don't let it be <0
Nobs[1] <- max(0, N0 + rnorm(1, 0, sdObs))
for (t in 1:T) {
R0 <- rnorm(1, muR0, sdR0)
K <- rnorm(1, muK, sdK)
N[t + 1] <- N[t] * exp(R0 * (1 - N[t]/K))
# add observation error to count
# and don't let it be <0
Nobs[t + 1] <- N[t + 1] + rnorm(1, 0, sdObs)
}
Nmax <- max(max(Nobs), max(N))
Nmin <- min(min(Nobs), min(N))
# you might change the upper limits of the y-axis
# below for the more variable scenarios
plot(1:(T + 1), N, "o",
ylim = c(0, 130),
xlab = 'Year',
main = c('avg N =', round(mean(N[10:T]), 2)))
abline(h = muK, col = 'red', lwd = 3)
points(1:(T + 1), Nobs, "o", col = 'blue')