A probabilistic analysis of a Beverton-Holt type discrete model: Theoretical and computing analysis

TL;DR

This paper introduces a randomized Beverton-Holt discrete model, deriving its probabilistic properties and steady states, supported by numerical examples to illustrate the theoretical findings.

Contribution

It provides a full probabilistic analysis of a stochastic Beverton-Holt model, including density functions and steady states, which is a novel extension of the classical model.

Findings

Derived the first and second probability density functions of the model

Identified the random steady state of the stochastic process

Validated theoretical results with numerical examples

Abstract

In this paper a randomized version of the Beverton-Holt type discrete model is proposed. Its solution stochastic process and the random steady state are determined. Its first probability density function and second probability density function are obtained by means of the random variable transformation method, providing a full probabilistic description of the solution. Finally, several numerical examples are shown.

Figures

Equations

$$x^{\prime }(t)=f(x(t),t),$$

$$x_{n+1}=\frac{\nu K_n{x}_n}{K_n+(\nu -1)x_n},n\ge 0,$$

$$x_{n+1}=\frac{a{x}_n}{1+b{x}_n},$$

$$x_{n+1}=\frac{\alpha /\beta }{1+\left({\mathrm{e}}^{-\alpha (n+1)}/(\gamma \beta )\right)}=\frac{{\mathrm{e}}^{\alpha }\alpha /\beta }{{\mathrm{e}}^{\alpha }+\left({\mathrm{e}}^{-\alpha n}/(\gamma \beta )\right)}=\frac{{\mathrm{e}}^{\alpha }\alpha /\beta }{1+\left({\mathrm{e}}^{-\alpha n}/(\gamma \beta )\right)+{\mathrm{e}}^{\alpha }-1}.$$

$$\begin{array}[]{ccl}x_{n+1}&=&\dfrac{\mathrm{e}^{\alpha}\alpha/\beta}{1+\mathrm{e}^{-\alpha n}/(\gamma\beta)}\dfrac{1+\mathrm{e}^{-\alpha n}/(\gamma\beta)}{1+\left(\mathrm{e}^{-\alpha n}/(\gamma\beta)\right)+\mathrm{e}^{\alpha}-1}\\ \\ &=&\dfrac{\mathrm{e}^{\alpha}\alpha/\beta}{1+\mathrm{e}^{-\alpha n}/(\gamma\beta)}\dfrac{1}{1+\dfrac{\mathrm{e}^{\alpha}-1}{1+\mathrm{e}^{-\alpha n}/(\gamma\beta)}}\\ \\ &=&\mathrm{e}^{\alpha}\dfrac{\alpha/\beta}{1+\mathrm{e}^{-\alpha n}/(\gamma\beta)}\dfrac{1}{1+\beta/\alpha\left(\mathrm{e}^{{\color[rgb]{0,0,1}\alpha}}-1\right)\dfrac{\alpha/\beta}{1+(\mathrm{e}^{-\alpha n})/(\gamma\beta)}}\\ \\ &=&\mathrm{e}^{\alpha}x_{n}\dfrac{1}{1+\beta/\alpha\left(\mathrm{e}^{{\color[rgb]{0,0,1}\alpha}}-1\right)x_{n}}\\ \\ &=&\dfrac{\mathrm{e}^{\alpha}x_{n}}{1+\dfrac{\beta}{\alpha}\left(\mathrm{e}^{{\color[rgb]{0,0,1}\alpha}}-1\right)x_{n}}.\end{array}$$

$$x_{n+1}=\dfrac{\mathrm{e}^{\alpha}x_{n}}{1+\dfrac{\beta}{\alpha}\left(\mathrm{e}^{{\color[rgb]{0,0,1}\alpha}}-1\right)x_{n}}.$$

$$z_{n+1}=\frac{1}{a}{z}_n+\frac{b}{a},$$

$$z_{n}=\left\{\begin{array}[]{llll}\displaystyle\left[\frac{1}{c}-\frac{b}{a-1}\right]a^{-n}+\frac{b}{a-1},&&&\text{if }a\neq 1,\\[17.07182pt] \displaystyle\frac{1}{c}+bn,&&&\text{if }a=1.\end{array}\right.$$

$$x_{n}=\left\{\begin{array}[]{llll}\displaystyle\frac{a^{n}(a-1)}{ba^{n}+\frac{1}{c}(a-1)-b},&&&\text{if }a\neq 1,\\[17.07182pt] \displaystyle\frac{1}{\frac{1}{c}+bn}&&&\text{if }a=1.\end{array}\right.$$

$$\lim_{n\rightarrow\infty}x_{n}=\left\{\begin{array}[]{llll}\displaystyle\frac{a-1}{b},&&&\text{if }a\neq 1,\\[17.07182pt] \displaystyle 0&&&\text{if }a=1.\end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle X_{n+1}=\frac{AX_{n}}{1+BX_{n}},\qquad n=0,1,2,\ldots,\\[11.38092pt] X_{0}=C,\end{array}\right.$$

$$\mathbb{E}\left[X_n\right]={\int }_{\mathbb{R}}x{f}_1^X(x,n)\mathrm{d}x,$$

$${\sigma }^2\left[X_n\right]={\int }_{\mathbb{R}}{x}^2{f}_1^X(x,n)\mathrm{d}x-{\left(\mathbb{E}\left[X_n\right]\right)}^2.$$

$${\Gamma }_X(n_1,n_2)=\mathbb{E}\left[X_{n_1}{X}_{n_2}\right]={\int }_{{\mathbb{R}}^2}{x}_1{x}_2{f}_2^X(x_1,n_1;x_2,n_2)\mathrm{d}{x}_1\mathrm{d}{x}_2.$$

$$C_X(n_1,n_2)={\Gamma }_X(n_1,n_2)-\mathbb{E}\left[X_{n_1}\right]\mathbb{E}\left[X_{n_2}\right].$$

$$f_{\mathbf{V}}(\mathbf{v})=f_{\mathbf{U}}\left(\mathbf{s}(\mathbf{v})\right)\left|J_n\right|,$$

$$J_{n}=\det\left(\frac{\partial\mathbf{s}^{\top}}{\partial\mathbf{v}}\right)=\det\left(\begin{array}[]{ccc}\dfrac{\partial s_{1}(v_{1},\ldots,v_{n})}{\partial v_{1}}&\cdots&\dfrac{\partial s_{n}(v_{1},\ldots,v_{n})}{\partial v_{1}}\\ \vdots&\ddots&\vdots\\ \dfrac{\partial s_{1}(v_{1},\ldots,v_{n})}{\partial v_{n}}&\cdots&\dfrac{\partial s_{n}(v_{1},\ldots,v_{n})}{\partial v_{n}}\\ \end{array}\right).$$

$$X_n=\frac{A^n(A-1)}{B{A}^n+\frac{1}{C}(A-1)-B},n=0,1,\dots$$

$$\begin{array}[]{ccccl}y_{1}&=&r_{1}\left(c,a,b\right)&=&\displaystyle\frac{a^{n}(a-1)}{ba^{n}+\frac{1}{c}(a-1)-b},\\[5.69046pt] y_{2}&=&r_{2}\left(c,a,b\right)&=&a,\\[5.69046pt] y_{3}&=&r_{3}\left(c,a,b\right)&=&b.\end{array}$$

$$\begin{array}[]{lllll}c&=&s_{1}\left(y_{1},y_{2},y_{3}\right)&=&\displaystyle\frac{y_{1}\left(y_{2}-1\right)}{y_{2}^{n}\left(y_{2}-1\right)-y_{3}y_{1}\left(y_{2}^{n}-1\right)},\\[5.69046pt] a&=&s_{2}\left(y_{1},y_{2},y_{3}\right)&=&y_{2},\\[5.69046pt] b&=&s_{3}\left(y_{1},y_{2},y_{3}\right)&=&y_{3},\end{array}$$

$$\left|J_3\right|=\left|\frac{\partial s_1}{\partial y_1}\right|=\left|\frac{{\left(y_2-1\right)}^2{y}_2^n}{{\left(y_2^n\left(y_2-1\right)-y_3{y}_1\left(y_2^n-1\right)\right)}^2}\right|,$$

$$\begin{array}[]{ccl}f_{Y_{1},Y_{2},Y_{3}}\left(y_{1},y_{2},y_{3}\right)&=&\displaystyle f_{C,A,B}\left(\frac{y_{1}\left(y_{2}-1\right)}{y_{2}^{n}\left(y_{2}-1\right)-y_{3}y_{1}\left(y_{2}^{n}-1\right)},y_{2},y_{3}\right)\\[11.38092pt] &&\times\displaystyle\left|\frac{\left(y_{2}-1\right)^{2}y_{2}^{n}}{\left(y_{2}^{n}\left(y_{2}-1\right)-y_{3}y_{1}\left(y_{2}^{n}-1\right)\right)^{2}}\right|.\end{array}$$

$$f_1^X(x,n)={\iint }_{{\mathbb{R}}^2}{f}_{C,A,B}\left(\frac{x\left(a-1\right)}{a^n\left(a-1\right)-bx\left(a^n-1\right)},a,b\right)\left|\frac{{\left(a-1\right)}^2{a}^n}{{\left(a^n\left(a-1\right)-bx\left(a^n-1\right)\right)}^2}\right|\mathrm{d}a\mathrm{d}b.$$

$$\begin{array}[]{ccccl}y_{1}&=&r_{1}\left(c,a,b\right)&=&\displaystyle\frac{a^{n_{1}}(a-1)}{ba^{n_{1}}+\frac{1}{c}(a-1)-b},\\[14.22636pt] y_{2}&=&r_{2}\left(c,a,b\right)&=&a,\\[5.69046pt] y_{3}&=&r_{3}\left(c,a,b\right)&=&\displaystyle\frac{a^{n_{2}}(a-1)}{ba^{n_{2}}+\frac{1}{c}(a-1)-b}.\end{array}$$

$$\begin{array}[]{lllll}c&=&s_{1}\left(y_{1},y_{2},y_{3}\right)&=&\displaystyle\frac{y_{1}y_{3}\left(y_{2}^{n_{1}}-y_{2}^{n_{2}}\right)}{y_{1}y_{2}^{n_{2}}\left(-1+y_{2}^{n_{1}}\right)-y_{3}y_{2}^{n_{1}}\left(-1+y_{2}^{n_{2}}\right)},\\[14.22636pt] a&=&s_{2}\left(y_{1},y_{2},y_{3}\right)&=&y_{2},\\[5.69046pt] b&=&s_{3}\left(y_{1},y_{2},y_{3}\right)&=&\displaystyle\frac{(-1+y_{2})\left(y_{3}y_{2}^{n_{1}}-y_{1}y_{2}^{n_{2}}\right)}{y_{1}y_{3}\left(y_{2}^{n_{1}}-y_{2}^{n_{2}}\right)},\end{array}$$

$$\left|J_3\right|=\left|\frac{\partial s_1}{\partial y_1}\frac{\partial s_3}{\partial y_3}-\frac{\partial s_1}{\partial y_3}\frac{\partial s_3}{\partial y_1}\right|=\left|\frac{(-1+y_2)y_2^{n_1+n_2}(y_2^{n_1}-y_2^{n_2})}{(y_1{y}_2^{n_2}(-1+y_2^{n_1})-y_3{y}_2^{n_1}(-1+y_2^{n_2}){)}^2}\right|,$$

$$\begin{array}[]{lcl}f_{Y_{1},Y_{2},Y_{3}}(y_{1},y_{2},y_{3})&=&\displaystyle f_{C,A,B}\left(\frac{y_{1}y_{3}\left(y_{2}^{n_{1}}-y_{2}^{n_{2}}\right)}{y_{1}y_{2}^{n_{2}}\left(-1+y_{2}^{n_{1}}\right)-y_{3}y_{2}^{n_{1}}\left(-1+y_{2}^{n_{2}}\right)},y_{2},\frac{(-1+y_{2})\left(y_{3}y_{2}^{n_{1}}-y_{1}y_{2}^{n_{2}}\right)}{y_{1}y_{3}\left(y_{2}^{n_{1}}-y_{2}^{n_{2}}\right)}\right)\\ \\ &&\displaystyle\times\left|\frac{(-1+y_{2})y_{2}^{n_{1}+n_{2}}(y_{2}^{n_{1}}-y_{2}^{n_{2}})}{\left(y_{1}y_{2}^{n_{2}}(-1+y_{2}^{n_{1}})-y_{3}y_{2}^{n_{1}}(-1+y_{2}^{n_{2}})\right)^{2}}\right|\end{array}$$

$$\begin{array}[]{lcl}f_{2}^{X}(x_{1},n_{1};x_{2},n_{2})&=&\displaystyle\int_{\mathbb{R}}f_{C,A,B}\left(\frac{x_{1}x_{2}\left(a^{n_{1}}-a^{n_{2}}\right)}{x_{1}a^{n_{2}}\left(-1+a^{n_{1}}\right)-x_{2}a^{n_{1}}\left(-1+a^{n_{2}}\right)},a,\frac{(-1+a)\left(x_{2}a^{n_{1}}-x_{1}a^{n_{2}}\right)}{x_{1}x_{2}\left(a^{n_{1}}-a^{n_{2}}\right)}\right)\\ \\ &&\qquad\displaystyle\times\left|\frac{(-1+a)a^{n_{1}+n_{2}}(a^{n_{1}}-a^{n_{2}})}{\left(x_{1}a^{n_{2}}(-1+a^{n_{1}})-x_{2}a^{n_{1}}(-1+a^{n_{2}})\right)^{2}}\right|\mathrm{d}a.\end{array}$$

$$X_{\infty }=\frac{A-1}{B}.$$

$$\begin{array}[]{ccccl}y_{1}&=&r_{1}\left(c,a,b\right)&=&\displaystyle\frac{a-1}{b},\\[5.69046pt] y_{2}&=&r_{2}\left(c,a,b\right)&=&c,\\[5.69046pt] y_{3}&=&r_{3}\left(c,a,b\right)&=&b,\end{array}$$

Full text

††journal: XXX

1 Introduction

One of the most interesting problems in population dynamics is modelling the changes in the size of the population over the time [1]. The breeding population of a species has been mathematically described using both continuous and discrete models. In the former case, the corresponding models are mainly based on differential equations while difference equations are the most used mathematical tools for the latter approach.

Explicit first-order differential equations are formulated in the form

[TABLE]

where the unknown, $x(t)$, represents the number or percentage of individuals of the population under study, while the right-hand size term, $f(x(t),t)$, is a specific function depending on the proposed model. In the the simplest model, usually termed Malthusian or exponential model, $f(x(t),t)=\alpha x(t)$. This model is satisfactory as long as the population is not too large with respect to its environment (resources). Otherwise, the accuracy of the model fails and, when $\alpha>0$, the Malthusian model predicts an exponential increase in the population as time goes on, which is unrealistic because of resources, like food, are always limited. This drawback is overcame by the most popular continuous model of the growth of population, usually referred to as the logistic o Verhulst-Pearl model formulated via a nonlinear differential equation where $f(x(t),t)=x(t)(\alpha-\beta x(t))$, being $\alpha,\beta>0$. In this model $\alpha$ represents the rate of the growth of the population if the resources were unlimited and the individuals did not affect one another, and the term $-\beta x^{2}(t)$ can be interpreted as the negative effect on the growth due to crowdedness and limited resources. This model was first published by Pierre Verhulst [2, 3]. This equation has been extensively used to model different problems in a variety of realms including Social Sciences, Biology, Ecology, etc. [4, 5], showing its great flexibility in applications. It is well-known that the general solution of logistic model is given by $x(t)=(\alpha/\beta)/(1+\left(\mathrm{e}^{-\alpha t}/(\gamma\beta)\right))$, where the parameter $\gamma$ is a real number that can be determined once an initial condition (the initial population) has been set.

Many times only data is available for discrete times, so it is interesting to have the solution for these time instants. In that case, the application of discrete models is natural and difference equations emerge. In this context, a very useful discrete equation in ecology is the non-autonomous Beverton-Holt equation[6]

[TABLE]

where $\nu\geq 1$, $K_{n}>0$ and $x_{0}>0$. The Beverton-Holt equation is of great importance in the fishery industry concerning the grow and exploitation of some species. The sequence $K_{n}$ is the carrying capacity and $\nu$ is the inherent growth rate. Sometimes the carrying capacity can be considered constant[6], it is $K_{n}=K$. This leads to the following difference equation

[TABLE]

where $a=\nu\geq 1$ and $b=\frac{1}{K}(\nu-1)=\frac{1}{K}(a-1)\geq 0$.

An interesting property is that Beverton-Holt equation (2) can be derived from logistic equation, as we will show below. Denoting an arbitrary time instant by $n$, and using the standard notation $x_{n}$ for the solution of the logistic model at $n$, we clearly have $x_{n}=(\alpha/\beta)/(1+\mathrm{e}^{-\alpha n}/(\gamma\beta))$. Now, let us observe that

[TABLE]

From this latter expression, we can derive a difference equation associated with $x_{n}$. Indeed, let us observe that

[TABLE]

Therefore,

[TABLE]

Observe that this model has the structure of Eq. (2), where $a=\mathrm{e}^{\alpha}>1$ and $b=\frac{\beta}{\alpha}\left(\mathrm{e}^{\alpha}-1\right)>0$, that belongs to the so-called Beverton-Holt type (sometimes referred to as Pielou logistic equation [7]).

Beverton-Holt type equation (2) is a nonlinear difference equation. By letting $x_{n}=1/z_{n}$ it is transformed into the linear equation

[TABLE]

whose solution, taking the initial condition $z_{0}=1/c$, is

[TABLE]

Or equivalently, in terms of unknown $x_{n}$ can be written as

[TABLE]

An important magnitude is the steady state, which is given by

[TABLE]

In real problems coefficients and initial conditions are not usually known exactly. This may be due to measurement errors or the inherent complexity associated to their own nature. So, it seems more realistic to consider that parameters and initial conditions are random variables instead of deterministic magnitudes.

Denoting random variables by capital letters, $A=\mathrm{e}^{\alpha}>1$ and $B=\frac{\beta}{\alpha}\left(\mathrm{e}^{\alpha}-1\right)>0$, the random Beverton-Holt model (3) can be written as

[TABLE]

where all the input parameters $A$,$B$ and $C$ are assumed to be absolutely continuous random variables defined on a common complete probability space ($\Omega,\mathcal{F},\mathbb{P}$). Hereinafter, we assume that $\mathbb{P}\left[\{\omega\in\Omega:\,A(\omega)>1\}\right]=1$, $\mathbb{P}\left[\{\omega\in\Omega:\,B(\omega)>0\}\right]=1$ and $\mathbb{P}\left[\{\omega\in\Omega:\,C(\omega)>0\}\right]=1$. Moreover, for the sake of completeness, we will suppose that $A,B,C$ are dependent random variables whose joint PDF is $f_{C,A,B}(c,a,b)$.

Solving a random difference equation means not only to calculate its exact solution stochastic process but also to determine its main statistical functions like the mean and the variance Additionally, an important goal is the computation of the 1-PDF, say $f^{X}_{1}(x,n)$, since from it we can compute all the one-dimensional statistical moments of any order. In particular, from the 1-PDF we can straightforwardly determine the mean function

[TABLE]

and the variance function

[TABLE]

The computation of the 1-PDF allows us to obtain a complete statistical description of the solution stochastic process at every time instant or period, say $n$. In general, a more challenging problem is to obtain the rest $n$-dimensional PDFs of the solution stochastic process for $n\geq 2$ (the so-called fidis: finite-dimensional distributions) because it usually involves complex computations. These functions provide us important probabilistic information. For example, the 2-PDF, $f^{X}_{2}(x_{1},n_{1};x_{2},n_{2})$, allows us to obtain a complete probabilistic description of the solution stochastic process at every arbitrary pair of periods, say $n_{1}$ and $n_{2}$. In particular, from the 2-PDF we can compute the correlation function. The correlation function gives a measure of linear statistical interdependence between $X_{n_{1}}$ and $X_{n_{2}}$

[TABLE]

Furthermore, $\Gamma_{X}(n_{1},n_{2})$ together with the mean function, permits us the computation of the covariance function

[TABLE]

The so-called RVT method is a powerful method that has been recently used to determine the PDF of the solution of some relevant differential and difference equations [8, 9, 10, 11]. The RVT technique permits to compute the PDF of a random vector which results from mapping another random vector whose PDF is known. For the sake of completeness, now we state the multidimensional version of the RVT technique that will be extensively applied throughout this paper.

Theorem 1

(Multidimensional version, [12, pp. 24–25]). Let $\mathbf{U}=(U_{1},\ldots,U_{n})^{\top}$ and $\mathbf{V}=(V_{1},\ldots,V_{n})^{\top}$ be two $n$-dimensional absolutely continuous random vectors. Let $\mathbf{r}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ be a one-to-one deterministic transformation of $\mathbf{U}$ into $\mathbf{V}$, i.e., $\mathbf{V}=\mathbf{r}(\mathbf{U})$. Assume that $\mathbf{r}$ is continuous in $\mathbf{U}$ and has continuous partial derivatives with respect to $\mathbf{U}$. Then, if $f_{\mathbf{U}}(\mathbf{u})$ denotes the joint PDF of vector $\mathbf{U}$, and $\mathbf{s}=\mathbf{r}^{-1}=(s_{1}(v_{1},\ldots,v_{n}),\ldots,s_{n}(v_{1},\ldots,v_{n}))^{\top}$ represents the inverse mapping of $\mathbf{r}=(r_{1}(u_{1},\ldots,u_{n}),\ldots,r_{n}(u_{1},\ldots,u_{n}))^{\top}$, the joint PDF of vector $\mathbf{V}$ is given by

[TABLE]

where $\left|J_{n}\right|$ is the absolute value of the Jacobian, which is defined by

[TABLE]

The layout of the paper is as follows. In Section 2 we compute explicit expressions of the 1-PDF (Subsection2.1) and the 2-PDF (Subsection2.2) of the solution stochastic process of Beverton-Holt model in terms of the joint PDF of the random inputs. In Section 3, the PDF of the equilibrium of Beverton-Holt model is determined. Section 4 is addressed to illustrate of our theoretical results via some numerical experiments assuming a wide range of statistical distributions for the input model parameters. Conclusions are drawn in Section 5.

2 Solving the randomized Beverton-Holt model

The goal of this section is to obtain the two first PDFs (1-PDF and 2-PDF) of the solution of (7), say $f_{1}^{X}(x,n)$ and $f_{2}^{X}(x_{1},n_{1};x_{2},n_{2})$, respectively. As $A$ is an absolutely continuous random variable, then $\mathbb{P}\left[\{\omega\in\Omega:\,A(\omega)=1\}\right]=0$, for all event $\omega\in\Omega$. As a consequence, taking into account that the solution of deterministic problem (3) (with $a=\mathrm{e}^{\alpha}>1$ and $b=\frac{\beta}{\alpha}\left(\mathrm{e}^{\alpha}-1\right)>0$) is given by (5), we obtain that the solution of random Beverton-Holt model (7) is given by

[TABLE]

and this solution is well-defined.

2.1 Computing the 1-PDF of the solution stochastic process

As previously indicated, in this section we will obtain the 1-PDF of (12) using the RVT method. To do this, given an arbitrary but fixed period $n$, we will apply Theorem 1 for the following choice of mapping $\mathbf{r}$

[TABLE]

So, the inverse mapping $\mathbf{s}$ of mapping $\mathbf{r}$ is given by

[TABLE]

and the absolute value of its Jacobian is

[TABLE]

which is different from zero since $y_{2}=a=\mathrm{e}^{\alpha}>1$. Observe that $0\neq y_{2}\neq 1$ since $A$ is an absolutely continuous RV.

Applying Theorem 1, the PDF of the random vector $(Y_{1},Y_{2},Y_{3})$ defined by mapping $\mathbf{r}$ is

[TABLE]

Now, marginalizing expression (13) with respect to $A$ and $B$ and being $n$ arbitrary, we obtain the 1-PDF of $X_{n}$

[TABLE]

2.2 Computing the 2-PDF of the solution stochastic process

To compute the 2-PDF of the solution stochastic process $X_{n}$, an analogous reasoning exhibited in the previous subsection will be applied. Let us consider two periods of time $n_{1},n_{2}\geq 0$ fixed, being $n_{1}\neq n_{2}$, and then we apply Theorem 1 choosing the mapping $\mathbf{r}$ as follows

[TABLE]

The inverse mapping $\mathbf{s}$ of $\mathbf{r}$ is given by

[TABLE]

which absolute value of Jacobian

[TABLE]

is different from zero since $y_{2}>1$ and $0\neq y_{2}\neq 1$, as it was justified previously. Applying the RVT technique, Theorem 1

[TABLE]

Now, marginalizing with respect to $A$ and taking $n_{1}$ and $n_{2}$ arbitrary, the 2-PDF is

[TABLE]

3 PDF of the steady state

An important issue in dealing with Beverton-Holt model is to determine the steady state. From the deterministic theory it is known that its steady state is given by (6). As $\mathbb{P}\left[\{\omega\in\Omega:\,A(\omega)=1\}\right]=0$, for all event $\omega\in\Omega$, the steady state of Beverton-Holt model (7) is given by

[TABLE]

This section is devoted to compute the PDF of (16). To achieve this objective, we define the mapping $\mathbf{r}$ by

[TABLE]

and using RVT technique (Theorem 1), it is easy to check that the PDF corresponding to $f_{Y_{1},Y_{2},Y_{3}}(y_{1},y_{2},y_{3})$ is given by

[TABLE]

Finally, the PDF of the steady state (16) is obtained marginalizing last expression, obtaining

[TABLE]

We will check via the example exhibited in next section that $f_{1}^{X}(x,n)$ given by (14) converges to $f_{X_{\infty}}(x)$ as $n$ increases.

4 Numerical experiments

In this section we illustrate our theoretical results via several numerical examples. Computations have been carried out using the software Mathematica®.

Example 1

This section is addressed to illustrate our previous theoretical results through an example. With this aim, we need to choose the distributions of random variables $C$, $A$ and $B$. We will assume that these random variables are independent. Therefore, its joint PDF can be factorized as the product of each individual PDF, $f_{C,A,B}(c,a,b)=f_{C}(c)f_{A}(a)f_{B}(b)$. As $X$ must be positive we have chosen for $C$ a uniform distribution in the interval $[0,1]$. For $A$ we have taken a uniform distribution in the interval $[1.1,2]$ and for $B$ a uniform distribution in the interval $[0.1,1]$

In Figure 1, the 1-PDF, $f_{1}^{X}(x,n)$, of the solution stochastic process, $X_{n}$, for different values of $n\in\{1,2,3,5,10,20\}$ and the PDF, $f_{X_{\infty}}(x)$, of the equilibrium point, $X_{\infty}$, are shown. As it has been previously pointed out, in this graphical representation one can observe that $f_{1}^{X}(x,n)$ tends to $f_{X_{\infty}}(x)$ when $n$ increases. In Figure 2, we have plotted separately the PDFs $f_{1}^{X}(x;20)$ and $f_{X_{\infty}}(x)$ to better highlight this behaviour. In this figure we can observe that both PDFs match.

In the left part of Figure 3, the mean of the solution stochastic process, $X_{n}$, and the threshold computed by the mean of the equilibrium have been represented. As it occurs with the 1-PDF, now we can observe that $\mathbb{E}\left[X_{n}\right]$ tends to $\mathbb{E}\left[X_{\infty}\right]$ as $n$ increases. In the right part a similar plot is presented for the standard deviation.

Figure 1. 1-PDF, $f_{1}^{X}(x;n)$, of the solution, $X_{n}$, with $n\in\{1,2,3,5,10,20\}$ and the PDF, $f_{X_{\infty}}(x)$, of the equilibrium or steady state $X_{\infty}$. Example 1.

Figure 2. 1-PDF, $f_{1}^{X}(x;n)$, of the solution, $X_{n}$ with $n=20$, and PDF, $f_{X_{\infty}}(x)$, of steady state $X_{\infty}$. Example 1.

Figure 3. Left: Points represent the mean of $X_{n}$, $\mathbb{E}\left[X_{n}\right]$, for different $n\in\{1,2,\dots,50\}$. Solid line represents the mean, $\mathbb{E}\left[X_{\infty}\right]$, of the steady state $X_{\infty}$. Right: Points represent the standard deviation of $X_{n}$, $\sigma\left[X_{n}\right]$, for different $n\in\{1,2,\dots,50\}$. Solid line represents the standard deviation, $\sigma\left[X_{\infty}\right]$, of the steady state $X_{\infty}$. Example 1.

In Figure 4, the 2-PDF of the solution stochastic process $X_{n}$ has been plotted in two cases: (1) $(n_{1},n_{2})=(1,2)$ and $(n_{1},n_{2})=(2,1)$. Finally, in Figure 5, the covariance surface, $C_{X}(n_{1},n_{2})$ of the solution stochastic process $X_{n}$ has been represented. This important deterministic function has been computed taking advantage of the explicit expressions of the 2-PDF $f_{2}^{X}(x_{1},n_{1};x_{2},n_{2})$ given in (15) together with expressions (8)–(9).

Figure 4. 2-PDF, $f_{2}^{X}(x_{1},n_{1};x_{2},n_{2})$, of the solution, $X_{n}$. Left: $n_{1}=1$ and $n_{2}=2$. Right: $n_{1}=2$ and $n_{2}=1$. Example 1.

Figure 5. Covariance function, $C_{X}(n_{1},n_{2})$, of the solution, $X_{n}$, for the values of $(n_{1},n_{2})\in[0,20]\times[0,20]$. Example 1.

Example 2

In contrast to Example 1, where we have considered that the involved input data are statistically independent, in this example we will illustrate our theoretical findings assuming that $A$, $B$ and $C$ are dependent random variables. Specifically, hereinafter we will assume that the random vector $(C,A,B)$ has a truncated Gaussian distribution $(C,A,B)\sim\text{N}_{T}(\mu,\Sigma)$, where $T=[0,1]\times[1.1,2]\times[0,1]$ is the domain of truncation, and the parameters $\mu$ and $\Sigma$, that represent the mean vector and the variance-covariance matrix, respectively, are given by

[TABLE]

In Figure 6, the 1-PDF, $f_{1}^{X}(x,n)$, of the solution stochastic process, $X_{n}$, for different values of $n\in\{1,2,4,8,16\}$ and the PDF, $f_{X_{\infty}}(x)$, of the equilibrium point, $X_{\infty}$, are shown. As it has been previously pointed out, in this graphical representation one can observe that $f_{1}^{X}(x,n)$ tends to $f_{X_{\infty}}(x)$ as $n$ increases. In Figure 7, we have plotted separately the PDFs $f_{1}^{X}(x,16)$ and $f_{X_{\infty}}(x)$ to show better this behaviour. In this figure we can observe that both PDFs match.

In Figure 8 (left), both the mean of the solution stochastic process, $X_{n}$, and the threshold computed by the mean of the equilibrium have been represented. As it occurs with the 1-PDF, now we can observe that $\mathbb{E}\left[X_{n}\right]$ tends to $\mathbb{E}\left[X_{\infty}\right]$ as $n$ increases. In Figure 8 (right), a similar plot is shown for the standard deviation.

Figure 6. 1-PDF, $f_{1}^{X}(x,n)$, of the solution, $X_{n}$, with $n\in\{1,2,4,8,16\}$ and the PDF, $f_{X_{\infty}}(x)$, of the equilibrium or steady state $X_{\infty}$. Example 2.

Figure 7. 1-PDF, $f_{1}^{X}(x,n)$, of the solution, $X_{n}$, with $n=16$, and PDF, $f_{X_{\infty}}(x)$, of the equilibrium or steady state $X_{\infty}$. Example 2.

Figure 8. Left: Points represent the mean, $\mathbb{E}\left[X_{n}\right]$, of $X_{n}$ for $n\in\{1,2,\dots,30\}$. Solid line represents the mean, $\mathbb{E}\left[X_{\infty}\right]$, of the steady state, $X_{\infty}$. Right: Points represent the standard deviation, $\sigma\left[X_{n}\right]$, of $X_{n}$ for $n\in\{1,2,\dots,30\}$. Solid line represents the standard deviation, $\sigma\left[X_{\infty}\right]$, of the steady state $X_{\infty}$. Example 2.

In Figure 9, the 2-PDF of the solution stochastic process, $X_{n}$, has been plotted in two cases: (1) $(n_{1},n_{2})=(1,2)$ and $(n_{1},n_{2})=(2,1)$. Finally, in Figure (10), the covariance surface, $C_{X}(n_{1},n_{2})$, of the solution stochastic process, $X_{n}$, has been represented on the square $(n_{1},n_{2})\in[0,16]\times[0,16]$. This important deterministic function has been computed taking advantage of the explicit expression of the 2-PDF, $f_{2}^{X}(x_{1},n_{1};x_{2},n_{2})$, given in (15) together with expressions (8)–(9).

Figure 9. 2-PDF, $f_{2}^{X}(x_{1},n_{1};x_{2},n_{2})$, of the solution, $X_{n}$. Left: $n_{1}=1$ and $n_{2}=2$. Right: $n_{1}=2$ and $n_{2}=1$. Example 2.

Figure 10. Covariance function, $C_{X}(n_{1},n_{2})$, of the solution, $X_{n}$, for the values of $(n_{1},n_{2})\in[0,16]\times[0,16]$. Example 2.

5 Conclusions

In this paper we have randomized the Beverton-Holt model. Then, we have provided a full probabilistic description of its solution stochastic process under very general assumptions on random input data. That description has been made through the 1-PDF and the 2-PDF of the discrete solution stochastic process. A full probabilistic description through its PDF is also given to the steady state. Finally, some numerical examples illustrating our theoretical results have been shown.

Acknowledgements

This work has been partially supported by the Ministerio de Economía y Competitividad grant MTM2017-89664-P. Ana Navarro Quiles acknowledges the postdoctoral contract financed by DyCon project funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 694126-DYCON).

Conflict of Interest Statement

The authors declare that there is no conflict of interests regarding the publication of this article.