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

Question

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}$$

What about the resulting distribution?

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.

Answer 1

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