Storing sum of values produced in a while-loop-CodePudding

I am plotting results from a simulation with multiple loops (bothfor and while) - consider the below simplified, reproducible example plotting the number of bunnies produced in a family for 100 simulated bunny families (each family starting with one bunny) and how long it takes for the family to go extinct. I would like to add a line in the plot representing the median number of bunnies from all simulated families and a band representing the interquartile range.

I specifically need help in determining where within the loops structure I can store the final family size of each simulated family, which is the sum of the initial bunnies and all subsequent bunnies_produced from each generation.

Below I am providing the current figure that I have (left) and the figure that I want (right; with made up values of the median and IQR).

I have tried initiating a variable totalbunnies and adding totalbunnies or totalbunnies[i] in various places, but it never correctly calculates the sum for each family (see commented out lines).

I do not need help with the actual visualization, I can add in the lines/bands, just need specific help on where to store the final family sizes for each simulation.

# Initiate plot
plot(NA, xlim = c(0, 10), ylim = c(0, 25), 
     xlab="Duration", ylab = "Number of Bunnies", frame = FALSE)

totalbunnies <- c() # initiate to store total number of bunnies in each sim

# Run simulations
set.seed(0506)
for(i in 1:100) {
  bunnies <- 1 
  t <- rep(0, 1) 
  duration <- t 
  while(bunnies > 0) {
    bunnies_produced <- rnbinom(bunnies, mu = 0.5, size = 0.25) 
    t.new <- numeric()
    for(j in 1:length(bunnies_produced)) {
      t.new <- c(t.new, t[j]   rgamma(bunnies_produced[j], shape = 0.25, rate = 0.5))
      #totalbunnies <- sum(totalbunnies, bunnies_produced)
    }
    bunnies <- length(t.new) 
    t <- t.new
    duration <- c(duration, t.new)
    # totalbunnies[i] <- sum(bunnies, bunnies_produced)
    # totalbunnies <- sum(totalbunnies, bunnies_produced)
  } 
  #totalbunnies[I] <- max(duration)
  lines(sort(duration), 1:length(duration), col = "maroon", lwd = 1) 
  points(max(duration), length(duration), col = "maroon", pch = 16)
}

CodePudding user response：

The sum of the total number of bunnies is found by adding totalbunnies[i] <- length(duration) instead of max(duration) just before the lines() argument