# "Difficult to please" rabbits and Fibonacci Numbers - A probabilistic variation

Fibonacci numbers were shown by the eponymous mathematician to be the solution of an admittedly idealized rabbit population growth problem.

Rabbits become fertile one month after their birth, after which they immediately, and successfully, biblically sleep with a mate: their gestation lasts one month as well and they give birth to a male and a female.

The number of rabbits over time is then found to be represented by Fibonacci numbers $$F_n = F_{n-1} + F_{n-2}$$ where the index stands for generation number.

I thought about a probabilistic variation on the theme, whereby not all couples mate successfully. Instead, they are rather picky and a constant proportion $p$ of couples does not generate progeny during each given generation. What will be the number, rather the distribution of couple of rabbits over time?

As far as the expected value, I believe this should follow the recursion $$F_n = F_{n-1} + p F_{n-2}$$

I think one could write $$F_n = F_{n-1} + \mathcal{B} (p; F_{n-2})$$ where $\mathcal{B} (p;n)$ stands for the binomial distribution for $n$ trials, with probability success $p$. Here I am unable to simplify the ensuing expression any further and would like to see if and how it could be done.

• Is the mathematician eponymous? Or are the number eponymous after the mathematician? And, if you're going to introduce randomness, why don't you include the number of rabbits born to each happy couple? Jul 4, 2017 at 16:26
• The mathematician in question being Leonardo of Pisa Jul 4, 2017 at 16:28
• @Scott, your variation is also interesting, but I have to learn to crawl first. With regards to your first point, I embellished my question with a hysteron-proteron, that you certainly did not fail to notice. Jul 4, 2017 at 16:28
• You are asking for the distribution of $F_n$ if $$F_n=F_{n-1}+\sum_{k=1}^{F_{n-2}}X_{k,n}$$ for every $n\geqslant2$, where (presumably) $F_0=F_1=1$ and (most certainly) $(X_{k,n})_{k,n}$ is i.i.d. Bernoulli $p$. Then $1\leqslant F_n\leqslant F_n^0$ almost surely, where $F^0_n$ is the usual $n$th Fibonacci number, $(F_n)$ is almost surely nondecreasing and unbounded, and each $E(F_n)$ can be computed recursively. But which kind of supplementary information do you hope for? I see no reason to expect a simple formula for the distribution of every $F_n$, if ever this was your aim...
– Did
Jul 4, 2017 at 18:07
• @Did, I am no hoping for any particular additional information. Would juts like to know what is the maximum that can be reasonably achieved. Possible aims, without reference to how plausible they are, a closed form for the distribution after $n$ generations, a formula for the variance, maybe a statement on the asymptotic convergence to a certain distribution, etc. Jul 5, 2017 at 7:32

A rigorous definition of the random process $(F_n)_{n\geqslant0}$ described in the question is that $F_0=F_1=1$ and, for every $n\geqslant0$, $$F_{n+2}=F_{n+1}+\sum_{k=1}^{F_n}Z_{n,k}$$ where the doubly-indexed family $(Z_{n,k})_{n\geqslant0,k\geqslant1}$ is i.i.d. Bernoulli with $P(Z_{n,k}=1)=p$ and $P(Z_{n,k}=0)=1-p$.
In turn, this random recursion can be rewritten under the form of the bivariate branching process $$Y_n=\begin{pmatrix}Y^1_n\\ Y^2_n\end{pmatrix}=\begin{pmatrix}F_{n+1}\\ F_n\end{pmatrix}$$ starting from $$Y_0=\begin{pmatrix}1\\1\end{pmatrix}$$ with reproduction mechanism $$Y^1_{n+1}=\sum_{k=1}^{Y^1_n}1+\sum_{k=1}^{Y^2_n}Z_{n,k}\qquad Y^2_{n+1}=\sum_{k=1}^{Y^1_n}1+\sum_{k=1}^{Y^2_n}0$$ Thus, the mean reproduction matrix of $(Y_n)$ is $$\begin{pmatrix}1&p\\1&0\end{pmatrix}$$ whose eigenvalues are $$\lambda=\frac{1+\varrho}2\qquad\mu=\frac{1-\varrho}2$$ where $$\varrho=\sqrt{4p+1}$$ with eigenvectors $$U_\lambda=\begin{pmatrix}\lambda\\1\end{pmatrix}\qquad U_\mu=\begin{pmatrix}\mu\\1\end{pmatrix}$$ A consequence is that, for each $\nu$ in $\{\lambda,\mu\}$, $$M_n^\nu=\frac{F_{n+1}+(\nu-1)F_n}{\nu^{n+1}}$$ defines a martingale $(M_n^\nu)_{n\geqslant0}$ starting from $M_0^\nu=1$. Now, $\lambda>1$ hence $M^\lambda$ is a positive martingale, bounded in $L^2$, in particular, $M^\lambda_n\to W_\lambda$ almost surely, for some nonnegative random variable $W_\lambda$ such that $E(W_\lambda)=1$. From there, follows the almost sure limit $$\lim_{n\to\infty}\frac{F_n}{\lambda^n}= \frac{\lambda}{2\lambda-1}W_\lambda=\frac{\lambda}{\varrho}W_\lambda$$ This convergence is rather sharp since, using the second eigenvalue $\mu$ with $-1<\mu<0$, one can show that, almost surely, $$\varrho\,F_n=\lambda^{n+1}W_\lambda-\mu^{n+1}W_\mu+\mu^nG_n$$ for some second random variable $W_\mu$ such that $E(W_\mu)=1$ and some sequence $(G_n)_{n\geqslant0}$ such that $E(G_n)=0$ and $G_n\to0$ almost surely.
Note finally that this approach also provides the (much easier) result that $$E(F_n)=\frac{\lambda^{n+1}-\mu^{n+1}}\varrho$$