Solving the Random Pielou Logistic Equation with the Random Variable Transformation Technique: Theory and Applications

TL;DR

This paper introduces a probabilistic approach to the Pielou logistic equation by randomizing it and applying the transformation of random variables, providing a full stochastic description and demonstrating its application through examples.

Contribution

It develops a novel probabilistic framework for the randomized Pielou logistic model using the transformation of random variables, including real data application.

Findings

Derived probability density functions for the stochastic process and steady state.

Provided a full probabilistic description of the model.

Demonstrated applications with numerical and real data examples.

Abstract

The study of the dynamics of the size of a population via mathematical modelling is a problem of interest and widely studied. Traditionally, continuous deterministic methods based on differential equations have been used to deal with this problem. However discrete versions of some models are also available and sometimes more adequate. In this paper, we randomize the Pielou logistic equation in order to include the inherent uncertainty in modelling. Taking advantage of the method of transformation of random variables, we provide a full probabilistic description to the randomized Pielou logistic model via the computation of the probability density functions of the solution stochastic process, the steady state and the time until a certain level of population is reached. The theoretical results are illustrated by means of two examples, the first one consists of a numerical experiment and the second one shows an application to study the diffusion of a technology using real data.

Figures

Tables

Table 1. Table 1: Number of mobile lines in Spain, x n subscript 𝑥 𝑛 x_{n} , at year n 𝑛 n during the range of years 1999–2015. Example 3.2 . Source: Spanish Commission of Markets and Competition (’Comisión Nacional de los Mercados y la Competencia’) [ 18 ] .

Year	$n$	$x_n$	Year	$n$	$x_n$
$1999$	$0$	$\mathrm{15\hspace{0.17em}003\hspace{0.17em}708}$	$2008$	$9$	$\mathrm{49\hspace{0.17em}623\hspace{0.17em}339}$
$2000$	$1$	$\mathrm{24\hspace{0.17em}265\hspace{0.17em}059}$	$2009$	$10$	$\mathrm{51\hspace{0.17em}052\hspace{0.17em}693}$
$2001$	$2$	$\mathrm{29\hspace{0.17em}655\hspace{0.17em}729}$	$2010$	$11$	$\mathrm{51\hspace{0.17em}389\hspace{0.17em}417}$
$2002$	$3$	$\mathrm{33\hspace{0.17em}530\hspace{0.17em}997}$	$2011$	$12$	$\mathrm{52\hspace{0.17em}590\hspace{0.17em}507}$
$2003$	$4$	$\mathrm{37\hspace{0.17em}219\hspace{0.17em}839}$	$2012$	$13$	$\mathrm{50\hspace{0.17em}665\hspace{0.17em}099}$
$2004$	$5$	$\mathrm{38\hspace{0.17em}622\hspace{0.17em}582}$	$2013$	$14$	$\mathrm{50\hspace{0.17em}158\hspace{0.17em}689}$
$2005$	$6$	$\mathrm{42\hspace{0.17em}693\hspace{0.17em}832}$	$2014$	$15$	$\mathrm{50\hspace{0.17em}806\hspace{0.17em}251}$
$2006$	$7$	$\mathrm{45\hspace{0.17em}675\hspace{0.17em}855}$	$2015$	$16$	$\mathrm{51\hspace{0.17em}067\hspace{0.17em}569}$
$2007$	$8$	$\mathrm{48\hspace{0.17em}422\hspace{0.17em}470}$	–	–	–

Equations

$$x^{\prime }(t)=x(t)(\alpha -\beta x(t)),\alpha ,\beta >0,$$

$$x_{n+1}=\frac{a{x}_n}{1+b{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,\\ \\ \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,\\ \\ \displaystyle\frac{1}{\frac{1}{c}+bn}&&&\text{if }a=1.\end{array}\right.$$

$${\mu }_X(n)=\mathbb{E}[X_n]={\int }_{-\infty }^{\infty }x{f}_1^X(x;n)\mathtt{d}x,$$

$${\sigma }_X^2(n)=\mathbb{V}[X_n]={\int }_{-\infty }^{\infty }{x}^2{f}_1^X(x;n)\mathtt{d}x-{\left({\mu }_X(n)\right)}^2,$$

$$\mathbb{E}[(X_n{)}^k]={\int }_{-\infty }^{\infty }{x}^k{f}_1^X(x;n)\mathtt{d}x,k=0,1,2,\dots$$

$$\mathbb{P}[a\le X_n\le b]={\int }_a^b{f}_1^X(x;n)\mathtt{d}x.$$

$$\mathbb{P}\left[|X_n-{\mu }_X(n)|\ge \lambda \right]\le \frac{({\sigma }_X(n){)}^2}{{\lambda }^2},\lambda >0,$$

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

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

$$\left\{\begin{array}[]{l}\displaystyle X_{n+1}=\frac{AX_{n}}{1+BX_{n}},\qquad n=0,1,2,\ldots\\ X_{0}=C,\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},\\ y_{2}&=&r_{2}\left(c,a,b\right)&=&a,\\ 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)},\\ a&=&s_{2}\left(y_{1},y_{2},y_{3}\right)&=&y_{2},\\ 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|,$$

$$f_{Y_1,Y_2,Y_3}\left(y_1,y_2,y_3\right)=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)\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|.$$

$$f_1^X(x;n)={\iint }_{\mathcal{D}(A,B)}{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|dadb,$$

$$\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},\\ y_{2}&=&r_{2}\left(c,a,b\right)&=&a,\\ y_{3}&=&r_{3}\left(c,a,b\right)&=&c.\end{array}$$

$$X_{\infty }=\underset{n\to \infty }{\lim }{X}_n=\underset{n\to \infty }{\lim }\frac{A^n(A-1)}{B{A}^n+\frac{1}{C}(A-1)-B}=\frac{A-1}{B}.$$

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

$$f_{X_{\infty }}(x)={\iint }_{\mathcal{D}(C,B)}{f}_{C,A,B}\left(c,xb+1,b\right)\left|b\right|dcdb,$$

$$\hat{X}=X_n=\frac{A^n(A-1)}{B{A}^n+\frac{1}{C}(A-1)-B}.$$

$$N=\frac{\log \left(\frac{\hat{X}(CB-A+1)}{C(B\hat{X}-A+1)}\right)}{\log (A)}.$$

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

$$f_N(n)={\iint }_{\mathcal{D}(A,B)}{f}_{C,A,B}\left(\frac{\hat{x}(1-a)}{\hat{x}b(a^n-1)+a^n(1-a)},a,b\right)\frac{\left|\hat{x}(a-1)a^n(1-a+\hat{x}b)\log (a)\right|}{{\left(\hat{x}b(a^n-1)+a^n(1-a)\right)}^2}dadb,$$

$$\underset{{\mu }_A,{\mu }_B,{\mu }_C,{\sigma }_A,{\sigma }_B,{\sigma }_C}{\min }\sum _{n=0}^{16}{\left(x_n-\mathbb{E}\left[X_n({\mu }_A,{\mu }_B,{\mu }_C,{\sigma }_A,{\sigma }_B,{\sigma }_C)\right]\right)}^2$$

$$\mathbb{E}\left[X_n({\mu }_A,{\mu }_B,{\mu }_C,{\sigma }_A,{\sigma }_B,{\sigma }_C)\right]={\int }_{-\infty }^{\infty }x{f}_1^X(x;n)dx.$$

$$\begin{array}[]{ccccc}\mu_{A}=1.4912,&&\mu_{B}=0.095109,&&\mu_{C}=1.76917,\\ \sigma_{A}=0.00531,&&\sigma_{B}=0.0025587,&&\sigma_{C}=0.0050285.\end{array}$$

Full text

††journal: XXX

1 Introduction

Some interesting problems in population dynamics are to model the changes in the size and to quantify the composition of the population along the time[1]. Mathematical models of populations can be used to accurately describe changes in a population and, more importantly, to predict future changes.

Population growth can be described by several mathematical models. The simplest model is the Malthusian one. It can predict an exponential increase in the population with time, which is unrealistic for long time. But it can be an adequate model in short-time studies. A typical application of Malthusian model is the study of the evolution of bacteria in a laboratory.

A more realistic and popular continous model of population growth is the logistic o Verhulst-Pearl equation

[TABLE]

where $x(t)$ is the size of the population at time $t$, $\alpha$ is the rate of population growth if the resources were unlimited and the individuals did not affect one another, and $-\beta x^{2}(t)$ represents the negative effect on the growth due to crowdedness and limited resources. This model was introduced by Pierre Verhulst [2, 3]. Although originally this equation was introduced to model populations, it has also been used to model different problems [4].

Most models that describe population dynamics are continuous. But many times only data is available for discrete times, so it is interesting to have the solution for those moments. In this case, a difference equation instead of a differential equation can be more appropriate to model this situation. For the sake of clarity, as time variable can be considered discrete or continuous, depending upon the context, hereinafter we will called period in the former case and time in the latter.

The discrete version of the logistic equation is known as Pielou logistic equation and it is stated as

[TABLE]

where $a>1$ and $b>0$[5, 6, 7]. By letting $z_{n}=1/x_{n}$, the nonlinear difference Equation (2) is transformed into the following linear difference equation

[TABLE]

Taking the initial condition $z_{0}=1/c$, the solution of Equation (3) is

[TABLE]

Thus, the solution of Pielou logistic equation (2) 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 RVs instead of deterministic values. Notice that hereinafter capital letters are used to denote a RV.

There are some recent interesting contributions concerning continuous random modelling in population dynamics[8, 9]. Discrete random models have not been widely studied in the literature. Recently, authors of this paper have made several contributions related to Markov models[10, 11] and linear difference equations[12].

The main objective of this paper is to construct a randomized version of the Pielou logistic equation. To the best of our knowledge, this problem has not been considered in the extant literature yet, but randomizing parameters of a model is a technique used in other contexts.

Solving a random difference equation means not only to calculate the exact solution of the stochastic process but also its main statistical functions. The computation of the 1-PDF, $f_{1}^{X}(x;n)$, allows us to have a complete statistical description of the solution. From the 1-PDF we can calculate easily the mean, variance and other higher statistical moments, by means of the following expressions, respectively

[TABLE]

[TABLE]

[TABLE]

Furthermore significant information such as the probability of the solution lies within a set of interest can be determined from the 1-PDF

[TABLE]

This improves the computation of rough bounds, like the one derived via Chebyshev’s inequality[13]

[TABLE]

usually applied in practice.

RVT method is a powerful technique that has been recently used by the authors to compute the 1-PDF of the solution of some differential and difference equations [14, 15, 16, 12]. The RVT technique permits to compute the PDF of a RV which results from mapping another RV whose PDF is known. The multidimensional version of the RVT technique is stated in Theorem 1.

Theorem 1** (RVT Multidimensional version[13, 15] )**

Let $\mathbf{U}=(U_{1},\ldots,U_{m})^{\top}$ and $\mathbf{V}=(V_{1},\ldots,V_{m})^{\top}$ be two $m$-dimensional absolutely continuous random vectors. Let $\mathbf{r}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ 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 probability density function of vector $\mathbf{U}$, and $\mathbf{s}=\mathbf{r}^{-1}=(s_{1}(v_{1},\ldots,v_{m}),\ldots,s_{m}(v_{1},\ldots,v_{m}))^{\top}$ represents the inverse mapping of $\mathbf{r}=(r_{1}(u_{1},\ldots,u_{m}),\ldots,r_{m}(u_{1},\ldots,u_{m}))^{\top}$, the joint probability density function of vector $\mathbf{V}$ is given by

[TABLE]

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

[TABLE]

Remark 1

In the context of solving random difference equations, when Theoerem 1 is applied the choice of random vector $\mathbf{V}$ (that defines the mapping $\mathbf{r}$) often can be made in different ways. The convenience of choosing one mapping against another depends heavily upon the feasibility of obtaining the inverse mapping $\mathbf{s}$ as well as the involved computations to obtain the 1-PDF of the solution stochastic process. Anyway, the 1-PDF obtained is equivalent independently of the mapping initialy chosen. This issue will be illustrated later (see Remark 2).

This paper is organized as follows. In Section 2, the Pielou logistic equation (2) is randomized. Then, the solution, the steady state and the time until a given population level is achieved are determined. Also a full probabilistic description of these magnitudes is provided via their 1-PDF and PDFs, respectively. In Section 3, the theoretical findings are illustrated by examples. Finally, conclusions are drawn in Section 4.

2 Randomized Pielou logistic equation

Random Pielou logistic equation can be written as

[TABLE]

where all the input parameters $A$,$B$ and $C$ are assumed to be absolutely continuous RVs defined on a common complete probability space ($\Omega,\mathcal{F},\mathbb{P}$). As a natural extension of its deterministic counterpart[5], we assume $\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$. For the sake of generality, hereinafter we will assume that $A,B,C$ are dependent RVs whose joint PDF is $f_{C,A,B}(c,a,b)$.

The main goal of this section is to obtain the 1-PDF of the solution of the random Pielou logistic equation (8), say $f_{1}^{X}(x;n)$. In this section we will also determine the PDF of another interesting quantities in dealing with the Pielou equation. Specifically, the steady state and the time until a given proportion is reached will be studied from a probabilistic standpoint.

2.1 1-PDF of the solution of the randomized Pielou equation

Inspired in the deterministic theory (2)–(4), by introducing the change of variable $Z_{n}=1/X_{n}$ in Equation (8), it is linearized and then solved. As $A$ is an absolutely continuous RV, then $\mathbb{P}\left[\{\omega\in\Omega:\,A(\omega)=1\}\right]=0$, for all event $\omega\in\Omega$. As a consequence, we obtain that the solution of random Pielou logistic problem (8) is given by

[TABLE]

and this solution is well-defined from a probabilistic point of view.

In order to compute the 1-PDF of (9), we can define the transformation $\mathbf{r}:\mathbb{R}^{3}\longrightarrow\mathbb{R}^{3}$,

[TABLE]

whose inverse mapping $\mathbf{s}=\mathbf{r}^{-1}$ is given by

[TABLE]

and the absolute value of the Jacobian of the inverse mapping, $\mathbf{s}$, is

[TABLE]

which is different from zero if $y_{2}=a>1$, which holds by hypothesis.

Applying RVT technique (Theorem 1) for an arbitrary but fixed value of $n$, the PDF of the random vector $(Y_{1},Y_{2},Y_{3})$ defined by mapping $\mathbf{r}$ is given by

[TABLE]

Finally, the 1-PDF of the solution of the randomized Pielou logistic model, $X_{n}$ is obtained by marginalizing expression (10) with respect to $A$ and $B$, being $n$ arbitrary,

[TABLE]

where $\mathcal{D}(A,B)$ stands for the domain of the random vector $(A,B)$.

Remark 2

Notice that another possibility to choose an adequate mapping $\mathbf{r}$ when applying Theorem 1 would be

[TABLE]

This choice would lead to an equivalent expression to $f_{1}^{X}(x;n)$ given in (11).

2.2 PDF of the steady state of the randomized Pielou equation

A key magnitude in population dynamics is to determine its long term behaviour. Mathematically it can be obtained calculating the limit as $n$ tends to infinite of the solution $X_{n}$. So, the steady state of the randomized Pielou equation, $X_{\infty}$, is

[TABLE]

Notice that, as $\mathbb{P}\left[\{\omega\in\Omega:\,A(\omega)>1\}\right]=1$ and $\mathbb{P}\left[\{\omega\in\Omega:\,B(\omega)>0\}\right]=1$ for all event $\omega\in\Omega$ , $X_{\infty}$ is a well-defined positive RV as is required for the size of a population.

By applying again the RVT technique (Theorem 1), considering the mapping $\mathbf{r}:\mathbb{R}^{3}\longrightarrow\mathbb{R}^{3}$,

[TABLE]

and, after some technical computations, we obtain the PDF of the RV steady state,

[TABLE]

where $\mathcal{D}(C,B)$ stands for the domain of the random vector $(C,B)$.

2.3 PDF of time until a given population size is reached

Another interesting question in dealing with discrete population models is to determine the distribution of the time where the size of the population reaches a certain specific value, $\hat{X}$,

[TABLE]

Again the RVT technique is very useful to answer this interesting question. In order to determine the PDF of the time where the size of the population reaches a certain specific value, say $f_{N}(n)$, we first isolate $n$ from expression (14), and as the obtained expression is a RV, we use capital letter notation to denote it, $N$,

[TABLE]

Defining an appropriate mapping, for example, $\mathbf{r}:\mathbb{R}^{3}\longrightarrow\mathbb{R}^{3}$,

[TABLE]

and applying the RVT method (Theorem 1), one can obtain the PDF of $N$,

[TABLE]

where $\mathcal{D}(A,B)$ stands for the domain of the random vector $(A,B)$.

3 Examples

In this section we will present two examples. In Subsection 3.1, we show a numerical example aimed to illustrate the theoretical results obtained in Section 2. In Subsection 3.2, a second example is presented concerning modelling. Although Pielou equation is usually applied to model population dynamics, it can model other interesting situations as the diffusion of a technology. In particular, we describe the dynamics of the number of mobile lines in Spain during the range of years 1999–2015 using real data by means of the randomized Pielou equation. We assume that $A$, $B$ and $C$ are independent RVs.

3.1 Numerical example

We consider the randomized Pielou model (8) and we assume the following probability independent distributions for the inputs. For the initial condition we choose a truncated normal distribution, ${C}\sim\text{N}_{]0,1[}(0.5;0.05)$, and for $A$ and $B$ uniform and beta distributions ${A}\sim\text{Un}([1.1,2])$ and ${B}\sim\text{Be}(2;3)$, respectively. We have chosen these distributions to illustrate the capability of our approach to deal with different probability distributions, although other distribution can be used too. All the computations have been carried out using the software Mathematica®[17].

We have calculated the 1-PDF of the solution stochastic process $X_{n}$, $f_{1}^{X}(x;n)$, given by (11). Also we have determined the PDF of the steady state, $f_{X_{\infty}}(x)$, given by (13). In Figure 1 it is plotted $f_{1}^{X}(x;n)$ for different values of $n$ and $f_{X_{\infty}}(x)$. We can observe that $f_{1}^{X}(x;n)$ tends to $f_{X_{\infty}}(x)$ as $n$ increases. For the sake of clarity, in Figure 2 it is plotted $f_{1}^{X}(x;15)$ and $f_{X_{\infty}}(x)$, and we can observe that both PDFs are similar.

In Figure 3 it is represented the expectation (left) and the standard deviation (right) of $X_{n}$. One can observe that they tend to the expectation and the standard deviation of the steady state, respectively.

3.2 Application: modelling the diffusion of a technology

This subsection is addressed to apply the theoretical results obtained in Section 2 to model a problem using real data. In particular, we consider the number of mobile lines in Spain including prepaid or postpaid contracts. We assume that this situation can be modelled by the random Pielou logistic equation (8). For this proposal we consider data provided by the Spanish Commission of Markets and Competition (’Comisión Nacional de los Mercados y la Competencia’)[18]. In Table 1 the total number of mobile lines in Spain in each year during the range of years 1999–2015 are collected. For computations, data and outputs are expressed in ten millions of mobile lines. So, all the results presented in this section are expressed in these units.

Since the model parameters $A$, $B$ and $C$ do not have a physical meaning and involve very complex factors that determine the dynamics of the number of mobile phones in Spain, we take advantage of the Central Limit Theorem in Probability to approximate them via Gaussian random variables. Therefore, let us take Gaussian distributions for the random inputs $A\sim\text{N}(\mu_{A};\sigma_{A})$, $B\sim\text{N}(\mu_{B};\sigma_{B})$ and $C\sim\text{N}(\mu_{C};\sigma_{C})$, where the means and standard deviations are determined by adjusting the real data $x_{n}$ to the theoretical expectation of random Pielou logistic model. For this proposal, we find a solution of the following optimization problem that consists of minimizing the square error between real data $x_{n}$ and punctual predictions via the expectation ($\mathbb{E}\left[X_{n}(\mu_{A},\mu_{B},\mu_{C},\sigma_{A},\sigma_{B},\sigma_{C})\right]$),

[TABLE]

where

[TABLE]

Using the command NMinimize of Mathematica® software for optimization, the following parameters are obtained

[TABLE]

In Figure 4 the number of mobile lines in Spain during the range of years 1999–2015, $x_{n}$, obtained from Table 1 is represented by blue points. To calculate the expectation and confidence intervals, first we calculate $f_{1}^{X}(x;n)$ with parameters (16). The expectation of $X_{n}$ is calculated by expression (5) and it has been plotted by a solid line. We can observe a good fitting between $x_{n}$ and the expected values $\mathbb{E}[X_{n}]$. In the same graphical representation, the 75% and 99% confidence intervals are plotted in dash-dotted lines. These confidence intervals have been computed in the following way. Firstly a value of the period $\hat{n}\geq 1$ and $\alpha\in(0,1)$ are fixed, and secondly $z_{1}=z_{1}(\hat{n})$ and $z_{2}=z_{2}(\hat{n})$ are determined such that

[TABLE]

Then, $(1-\alpha)\times 100\%$-confidence interval is specified by

[TABLE]

As we can observe most of the data lies in the confidence intervals. Only two points are outside the corresponding confidence intervals. This is a good result taking into account the uncertainties in the market of mobile lines due to people preferences and economic oscillations.

A key feature of the previous methodology with respect to other available approaches is that it permits to construct confidence intervals considering the real probability distribution of the model output rather than using asymptotic approximations, which usually rely in the Gaussian distribution via rules of the type mean plus/minus two standard deviations. Notice that using the methodology proposed in this paper, and assuming appropriate distributions for the input data (in our case we have assumed that $A$, $B$ and $C$ are Gaussian), we are able to compute the theoretical probability distribution of the model output (which in our case is not Gaussian), and from it, to determine precise confidence intervals using a prefixed but arbitrary confidence level.

In Figure 5 the 1-PDF of $X_{n}$, $f_{1}^{X}(x;n)$, at different fixed periods $n$ are plotted. We can observe that these 1-PDFs tend, as $n$ increases, to the PDF of the steady state, $f_{X_{\infty}}(x)$, calculated by expression (13) and represented in light blue colour.

In order to better visualize the tendency of $f_{1}^{X}(x;n)$ to $f_{X_{\infty}}$, as $n$ increases, in Figure 6 we have plotted $f_{1}^{X}(x;16)$ together with $f_{X_{\infty}}(x)$. We can observe the mentioned tendency.

An interesting question that arises in modelling the dynamics of the number of mobile lines is to know the period $n$ when there will be a certain number of mobile phones. This allows to predict the period when mobile telecommunication infrastructures might be necessary to expand, for example. The PDF of this time is calculated by (15). In Figure 6 it is plotted the distribution of time until a given proportion of the population has got line mobile. For example, if we observe the PDF for $\hat{x}=5$, we can derive that there were about $50\,000\,000$ line mobiles in Spain at year 2008 ($n=9$).

4 Conclusions

In this paper a randomization of the Pielou logistic model to study the dynamics of a population has been considered. Important goals in this type of equations are to compute its solution, to study the steady state and to compute the time until a certain level of population is reached. A full probabilistic description of these randomized magnitudes has been provided. This has been achieved taking advantage of the RVT technique.

RVT technique allows us to calculate the PDF from a particular mapping. In each case studied throughout this paper, appropriate mapping have been chosen to reduce the computational cost.

Theoretical findings have been illustrated via two examples. The first one consists in some numerical experiments to illustrate our results from a mathematical standpoint. While the second one shows a real application to modelling the dynamics of the number of mobile lines in Spain during the range of years 1999–2015. We can conclude that the obtained results are satisfactory.

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.