2

I am trying to simulate populations in R, but my code currently includes nested for loops, which I would like to replace with apply loops so my code will run faster, though I'm not quite sure how. I've already read some existing topics about this, but I think my situation is different because I'm not just projecting growth, but reproduction and survivorship.

I have two vectors containing probabilities of survival and reproduction. I'll simulate some values for these.

St = rep(c(0.82, 0.89, 0.93), 3)
Bt = rep(c(0.13, 0.24, 0.17), 3)

What I want to do now is create a data frame with variables year, survivors, births, population size, and population growth. Year is trivial. The way I'm simulating the others is by killing off some of the starting population with whatever the survival probability is for that year, then applying the probability of reproduction to the survivors to get the number of births. Population size and growth are calculated from these. I have done this with two for loops nested inside a while loop.

count <- 0
years <- 9
population <- 200
survivors <- vector()
births <- vector()
growth <- vector()

while (count < years) {
count <- count + 1
for (i in St) {
  survivors[count] <- rbinom(n=1, size=population[count], prob=i)
}
for (i in Bt) {
  births[count] <- rbinom(n=1, size=survivors[count], prob=i)
}
growth[count] <- (population[count] - survivors[count]) + births[count]
population[count+1] <- population[count] + growth[count]
}
population <- population[2:length(population)]
year <- 1:years
data <- data.frame(year, births, survivors, population, growth)

Is there any way to generate this simulated time series using apply loops instead? I'm kind of stumped because the result of the births for loop in each year depends on the outcome of the survivors for loop.

Community
  • 1
  • 1
ecogrammer
  • 73
  • 8
  • you probably can't use `apply` to eliminate the `while` loop, but you can speed things up a bit by eliminating the `for` loops: try e.g. `survivors <- rbinom(n=length(survivors), size=population, count=seq(St))` – Ben Bolker Mar 13 '14 at 17:13
  • `apply` is not the first place I'd look for speed improvements. http://stackoverflow.com/questions/2908822/speed-up-the-loop-operation-in-r/8474941#8474941 If you really need speed, particularly for a simulation exercise, I'd consider `Rcpp`. – Ari B. Friedman Mar 13 '14 at 17:18
  • The problem with eliminating the for loops is that `population` in `rbinom` is no longer sensitive to changes in population size due to survivorship and reproduction. It just uses the starting population size, 200, and I get runaway growth. The advantage of using `population[count]` is that it represents the population size at time t - 1. – ecogrammer Mar 13 '14 at 17:40
  • You should have a look at package `simecol` for simulating models like yours. Make a single simulation for each component of St and Bt and combine the results if required after the running all variants into a single dataframe. I don't know if that would be faster but it would make it easier to understand and modify. You could also consider putting some or all of the calculations in a single function and then user the `compiler` package to compile the function to byte code, which could give some speed up. – Bhas Mar 13 '14 at 20:13

0 Answers0