# Set-up parameters for stochastic exponential growth N0 <- 20 yrs <- 300 mu.Lambda <- 1.06 sd.Lambda <- sqrt(0.05) sd.Lambda # calculate variance and mean for stochastic 'r' or 'r.bar' r.bar.var <- (sd.Lambda ^ 2) / (mu.Lambda ^ 2) r.bar <- log(mu.Lambda) - r.bar.var / 2 r.bar sd.r.bar <- sqrt(r.bar.var) sd.r.bar # plot the expected trajectory if there were no temporal variation # working on log scale for population size plot(1:(yrs+1), log(N0 * mu.Lambda ^ (0:yrs)), 'l', col='blue', ylim=c(0, 20), xlab="Time", ylab="ln N") # add a line for the expected stochastic trajectory # again working on log scale for population size lines(1:(yrs+1), log(N0 * exp(r.bar) ^ (0:yrs)), 'l', col='green') #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%# # Work with simulated stochastic trajectories # #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%# # choose how many stochastic trajectories to run sims <- 5 # create object to use for storing random growth rates rand.r <- matrix(data = NA, nrow = yrs, ncol = sims) # create object to use for storing population sizes N <- matrix(data = NA, nrow = yrs+1, ncol = sims) # store starting value in 1st row of each column N[1, ] <- N0 # generate dynamics and store values for (j in 1:sims) { for (i in 1:yrs) { # generate a random growth rate rand.r[i, j] = rnorm(1, r.bar, sd.r.bar) # apply growth rate to population for 1 time step N[i+1, j] = N[i, j] * exp(rand.r[i, j]) } # add simulated trajectory for each simulation # We'll slow things down to better see each simulation trajectory Sys.sleep(0.25) lines(1:(yrs+1), log(N[, j]), 'l', col='red') } # Re-plot the expected lines in thicker lines over red lines lines(1:(yrs+1), log(N0 * mu.Lambda ^ (0:yrs)), 'l', col = 'blue', lwd = 3) lines(1:(yrs+1), log(N0 * exp(r.bar) ^ (0:yrs)), 'l', col = 'green', lwd = 3) # calculate realized growth rates from simulations RGR <- (N[(yrs+1), ] / N[1, ]) ^ (1/yrs) RGR mean(RGR) exp(r.bar) # calculate the sum of 'r' values for each simulation (i.e., by column) # and check that these can be used along with N0 # to get final population size r.sums <- apply(rand.r, 2, sum) alt.N <- exp(log(N0) + r.sums) diffs <- alt.N - N[(yrs + 1), ] # should be ~0 summary(diffs) # check that all are ~0