The first question is a variation of the Fibonacci sequence. Just read the first paragraph.
Nine hundred years ago the Italian businessman Leonardo Bonacci (known as Fibonacci) created a famous model to predict the number of rabbits to breed. We want to make an improved model, by assuming the following hypotheses:
Like the original sequence, this sequence can be calculated using a recursive function.
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 outputbaby_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:
<- function(n) {
baby_rabbits 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:
<- function(n) {
baby_rabbits # Write your code here
<- rep(0,n)
beybi 1] <- 1
beybi[if(is.na(beybi[4])){}
else {
for (i in rep(4:length(beybi))) {
<- sum(beybi[(i-3):(i-1)])-beybi[i-2]
beybi[i]
}
}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
<- function(n) {
baby_rabbits <- rep(0,n)
beybi 1] <- 1
beybi[for (i in rep(4:length(beybi))) {
<- sum(beybi[(i-3):(i-1)])-beybi[i-2]
beybi[i]
}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:
<- function(n) {
baby_rabbits <- rep(0,n)
beybi 1] <- 1
beybi[for (i in rep(4:length(beybi))) {
<- beybi[i-3]+beybi[i-1]
beybi[i]
}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_rate*baby[i-1]
growth <- mating_rate*adult[i-1]
mating <- birth_rate*pregnant[i-1] birth
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.
<- birth - growth
d_baby[i] <- growth - mating
d_adult[i] <- mating d_pregnant[i]
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
<- birth_rate*pregnant[i-1] - growth_rate*baby[i-1]
d_baby[i] <- growth_rate*baby[i-1] - mating_rate*adult[i-1]
d_adult[i] <- mating_rate*adult[i-1] d_pregnant[i]
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 andgrowth_rate
is 0.5?
<- breeding(N=24, baby_ini=1, adult_ini=0, pregnant_ini=0,
new_rabbits 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
innew_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.
<- lm(log(baby) ~ Month, data = new_rabbits)
model 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.