1.1 Recursive solution

The first part asked explicitly for a recursive function. The question says

Write a recursive function in R that has input n and output baby_rabbits(n). You can test your function by checking that at the end of the first year there are 144 pairs of rabbits.

If you understand the idea of recursion, and you have seen the Fibonacci example, it is easy to adapt. There are two exit conditions: n<=1 and n<=3. In both of these cases the answer is easy. In the rest of the cases you need to use the function inside the function.

A typical correct answer looks like this:

baby_rabbits <- function(n) {
    if(n<=1) {
        return(1)
    } else if(n<=3){
        return(0)
    } else {
        return(baby_rabbits(n-1) + baby_rabbits(n-3))
    }
}

One wrong answer is this:

baby_rabbits <- function(n) {
  # Write your code here
  beybi <- rep(0,n)
  beybi[1] <- 1
  if(is.na(beybi[4])){}
  else {
    for (i in rep(4:length(beybi))) {
      beybi[i]  <-  sum(beybi[(i-3):(i-1)])-beybi[i-2]
    }
  }
  return(beybi[n])
}

This is not a recursive function. Moreover, it is over-complicated. For instance, it still has the comment # Write your code here. The comment should be deleted.

Then there is an if(is.na(beybi[4])){} question. Since the vector beibi is created with the value 0 on each element, then the question will always be false, and the else part will always be executed. So there is no need for this if(){}. Moreover, the code block between if() and else is empty. So it does nothing, and it should be deleted.

A little better version, still not recursive, would be

baby_rabbits <- function(n) {
  beybi <- rep(0,n)
  beybi[1] <- 1
  for (i in rep(4:length(beybi))) {
    beybi[i]  <-  sum(beybi[(i-3):(i-1)])-beybi[i-2]
  }
  return(beybi[n])
}

Still, the formula is weird. There is a sum, and then there is a minus. It is much easier to understand as this:

baby_rabbits <- function(n) {
  beybi <- rep(0,n)
  beybi[1] <- 1
  for (i in rep(4:length(beybi))) {
    beybi[i]  <-  beybi[i-3]+beybi[i-1]
  }
  return(beybi[n])
}

Why make things hard when they can be easy?

1.2 System-modeling solution

Again, a classical question. Basic dynamic systems modeling and simulation.

To solve this we can model the rabbits as a system. We have three kinds of rabbits: baby, adult and pregnant rabbits. In each month baby rabbits grow, so the number of adult rabbits increases and the number of baby rabbits decreases. The same happens with adult and pregnant rabbits. On the other side, each of pregnant rabbit produces a new baby rabbit, so the number of baby rabbits increases by the number of pregnant rabbits.

You only need to fill up the part of d_baby, d_adult and d_pregnant. As described in classes and in the complementary text, we start by multiplying the things going into the boxes. There are three boxes, so we have three formulas

growth <- growth_rate*baby[i-1]
mating <- mating_rate*adult[i-1]
birth  <- birth_rate*pregnant[i-1]

This is very easy, since each box has only one incoming arrow. It cannot be more easy. Literally.

Then we write the formulas for the circles. It is important to not forget to use the correct indices (that is, the i inside the brackets). Each incoming arrow is positive, each outgoing one is negative.

d_baby[i]     <- birth - growth
d_adult[i]    <- growth - mating
d_pregnant[i] <- mating

The two arrows of birth into pregnant have opposite signs so they cancel each other.

These two parts can be combined into one. This is the second way of writing the good solution

d_baby[i]     <- birth_rate*pregnant[i-1] - growth_rate*baby[i-1]
d_adult[i]    <- growth_rate*baby[i-1]    - mating_rate*adult[i-1]
d_pregnant[i] <- mating_rate*adult[i-1]

1.3 Change parameters to find new_rabbits

What happens if the pregnant rabbits produce more baby rabbits every month, but not all of the babies survive to be adults? For example, what happens if birth_rate is 2 and growth_rate is 0.5?

new_rabbits <- breeding(N=24, baby_ini=1, adult_ini=0, pregnant_ini=0,
                    growth_rate=0.5, mating_rate=1, birth_rate=2)
plot(baby ~ Month, data=rabbits)
points(baby ~ Month, data=new_rabbits, pch=2)
legend("topleft", legend=c("rabbits", "new_rabbits"), pch=1:2)
title(main="Two breedings of rabbits")

This question is very easy and most people answered it correctly. The main confusion was on how to label the axis, and how to make the legend. For example, some people tried this:

plot(rabbits$Month, rabbits$baby, col = 1)
points(new_rabbits$baby, col = 2)
legend("topleft", legend = c("rabbits             BLACK","new_rabbits    RED"))

We do not need to say col=1 or pch=1, since these are the default values.

1.4 Find the growth rate of new_rabbits

As with most unbounded growth processes, this system seems to grow exponentially. Write the code to fit a linear model for the column baby in new_rabbits and calculate the growth rate from the fitted coefficients. You may need to transform the data.

This is direct from what we learned in CMB1. The only thing to be aware is the exponential part. This means that we have to use log(baby) in the model, and then use exp() on the model’s coefficient.

model <- lm(log(baby) ~ Month, data = new_rabbits)
exp(coef(model))

Interestingly, we cannot fit the same model to data in rabbits, since some of the values in rabbits$baby are 0, and we get in trouble when we evaluate the logarithm of 0.