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 animproved 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

recursivefunction 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:

```
<- 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,adultandpregnantrabbits. 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 and`growth_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`

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.

```
<- 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.