# Project sizes of population growing with a finite per capita growth rate that varies stochastically
# Assumes the finite rate is sampled from a binomial distribution every time step
# By R. Gomulkiewicz (8 Sep 2015)

initial.size <- 10 #fixed starting population size
reps <- 5 # number of replicates
steps <- 10 #number of generations (iterations) for each "replicate"

# Create function that computes the population size at the next time step assuming each individual has
# a binomially-distributed number of descendants, with mean 1+r and max number of offspring = max.off
# This is a binomial distribution with parameters N = max.off and p = (1+r)/max.off
# "rbinom(n,N,p)" is a built-in R function that generates a list of n binomial random variables with  parameters N and p
# "sum(l)" is a built-in R function that sums the elements in a list l
next.size.rand <- function(n, max.off,r) {sum(rbinom(n,max.off,(1+r)/max.off))}

r <- 0.01 #expected intrinsic rate of growth or decline (a fixed parameter)
max.off <- 2  #this is the maximum number of descendants per individual		

# create a (steps+1) x reps table to hold the population sizes for "reps" replicates
# each column contains the population sizes, including the starting size
nstoch <- matrix(0,steps+1,reps)

#store the initial population size in the first position (row) of each column
nstoch[1,] <- initial.size

#double loop that projects the population sizes for the replicates starting from initial.size
# results will be stored in the 1st column
for(j in 1:reps){
for(i in 1:steps){ nstoch[i+1,j] = next.size.rand(nstoch[i,j],max.off,r)}
}

#plot all three trajectories together
matplot(0:steps,nstoch,type="l",xlab="time step",ylab="population size")