Generalized Beta Models and Population Growth: So Many Routes to Chaos

Abstract

Logistic and Gompertz growth equations are the usual choice to model sustainable growth and immoderate growth causing depletion of resources, respectively. Observing that the logistic distribution is geo-max-stable and the Gompertz function is proportional to the Gumbel max-stable distribution, we investigate other models proportional to either geo-max-stable distributions (log-logistic and backward log-logistic) or to other max-stable distributions (Fréchet or max-Weibull). We show that the former arise when in the hyper-logistic Blumberg equation, connected to the Beta $(p,q)$ function, we use fractional exponents $p-1=1\mp 1/\alpha$ and $q-1=1\pm 1/\alpha$, and the latter when in the hyper-Gompertz-Turner equation, the exponents of the logarithmic factor are real and eventually fractional. The use of a BetaBoop function establishes interesting connections to Probability Theory, Riemann–Liouville’s fractional integrals, higher-order monotonicity and convexity and generalized unimodality, and the logistic map paradigm inspires the investigation of the dynamics of the hyper-logistic and hyper-Gompertz maps.

1. Introduction

Forecasting the size of populations is essential for resource management, either in demographic studies, for instance, applied to the dynamics of overlapping generations and the equilibrium needed to preserve social security (de la Croix and Michel [1]; Michel et Wigniolle [2]), forestry (Yang et al. [3]; Payandeh and Wang [4]), Life Sciences and Medicine (Laird [5]; Laird et al. [6]; Norton et al. [7]; Bajzer [8]; Waliszewski et al. [9,10]; Waliszewski and Konarski [11,12,13]; Molski and Konarski, [14], and many others), and growth modeling (Tjørve and Tjørve [15], and references therein).

Most natural populations, after an initial period of exponential expansion followed by an almost linear growth, reach a stable sustainable growth status, for which the logistic function provides a good fit. There are, however, populations whose expansion is immoderate, for instance, the growth of cancer tumor cells, or the spread of epidemics, with the depletion of resources eventually causing ultimate extinction. The most widely used mathematical models to depict those frameworks are the Verhulst logistic and the Gompertz (Gumbel) population dynamics equations.

The main purpose of this review is to bring together different areas of Mathematics to explore the relevance of some forms of extreme growth. After showing that the widely used logistic and Gompertz models, sketched below, are connected to geo-max-stable (GMS) and to max-stable (MS) extreme value (EV) distributions, we use fractional exponents in the extensions due to Blumberg [16] and to Turner et al. [17,18] to explicit other forms of extreme growth and the dynamics of their associated hyper-logistic and hyper-Gompertz maps.

Important issues to achieve that main purpose are briefly stated next. In Section 2, it is provided an overview of BetaBoop functions,

$$G_{p,q,\mathsf{\Phi },\mathsf{\Psi }}\left(x\right)=x^{p-1}{(1-x)}^{q-1}{\left[-ln(1-x)\right]}^{\mathsf{\Phi }-1}{\left(-lnx\right)}^{\mathsf{\Psi }-1}{\mathbb{I}}_{(0,1)}\left(x\right),$$

(1)

for $p,q,\mathsf{\Phi },\mathsf{\Psi }>0$ such that ${\int }_0^1{G}_{p,q,\mathsf{\Phi },\mathsf{\Psi }}\left(x\right)\mathrm{d}x<\infty$ (see Appendix A). We discuss in detail the expansion role played by the factors $x^{p-1}$ and ${[-ln(1-x)]}^{\mathsf{\Phi }-1}$, the compression role played by the factors ${(1-x)}^{q-1}$ and ${(-lnx)}^{\mathsf{\Psi }-1}$, and the connection of $G_{p,2,1,1}$ and of $G_{p,1,1,2}$ to power laws, which are relevant for understanding the connection of hyper-logistic and hyper-Gompertz differential equations (DEs) to several extreme growth models.

BetaBoop functions and BetaBoop random variables (RVs) were initially introduced in Brilhante et al. [19], and their relevance in the investigation of forms of extreme growth has been carried out in Brilhante et al. [20,21] and in Gomes et al. [22].

In Section 3, a description of the unifying role of BetaBoop functions in important issues of fractional calculus, such as the Riemann–Liouville integral, generalized monotonicity, convexity and unimodality, and issues in Probability Theory such as the existence of non-normal/non-Cauchy additive stable laws, is performed. In Section 4, a preliminary discussion of logistic and of Gompertz growth and of the nonlinear maps arising from the Verhulst and from the Gompertz DEs, and a brief sketch of Blumberg hyper-logistic and of Turner hyper-Gompertz growth models is put forward. Section 5 is a brief sketch of max-stability and of geo-max-stability, highlighting the fact that the classical logistic growth and the Gompertz growth are closely related to probabilistic EV models. In Section 6, a general theory of growth with EV contours is discussed. Section 1.3 is a more detailed description of the review organization, and the introductory paragraph of Section 6 further highlights the purposes and most relevant topics discussed.

1.1. From Fibonacci Unbound Growth to Verhulst Sustainable Logistic Growth and Gompertz Extreme Growth

Fibonacci (c. 1170–c. 1250) in his Liber Abaci presented the first mathematical model for population growth. From rather crude assumptions, his model for the progeny of an initial couple of rabbits is the sequence of Fibonacci’s numbers ${\left\{F_n\right\}}_{n\in \mathbb{N}}$ starting with $\left\{0,1\right\}$ such that $F_{n+2}=F_n+F_{n+1}$. Using generating functions, it is easy to obtain Binet’s formula $F_n=\frac{{\left(1+\sqrt{5}\right)}^n-{\left(1-\sqrt{5}\right)}^n}{2^n\sqrt{5}}$.

The main drawback of Fibonacci’s model is the fact that the population will ultimately be infinite (curiously, von Foerster et al. [23] predicted that on Doomsday: Friday, 13 November, A.D. 2026, the human population will become infinite).

Denote by $N\left(t\right)$ a population’s size at time t. Imposing some natural regularity conditions on $N\left(t\right)$, namely that ${\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}N\left(t\right)=\sum _{k=1}^{\infty }{A}_k{\left[N\left(t\right)\right]}^k}$, Verhulst [24,25,26] used the second-order approximation $\frac{\mathrm{d}}{\mathrm{d}t}N\left(t\right)=A_{1}N\left(t\right)+A_2{\left[N\left(t\right)\right]}^2$, with $A_1>0$ and $A_2<0$, which can be rewritten as

$$\mathrm{Verhulst}\mathrm{equation}:\frac{\mathrm{d}}{\mathrm{d}t}N\left(t\right)=rN\left(t\right)\left[1-\frac{N\left(t\right)}{K}\right].$$

(2)

When modeling natural breeding populations, $r>0$ is the Malthusian instantaneous growth rate parameter, and the carrying capacity $K>0$ is the equilibrium limit size of the population. This logistic Verhulst equation is more realistic to model sustainable growth: as long as abundant resources allow, there is an initial period of exponential growth, which is followed by more moderate approximately linear growth. Finally, the limitation of natural resources (or success of predators or of competing species) enforces very moderate growth, maintaining the population size and consumption of resources in sustainable equilibrium.

It has been a shock when May [27] published a short note on the very complicated dynamics of the nonlinear map $x_{n+1}=r{x}_n(1-x_n)$ resulting from the logistic DE by using Euler’s algorithm to obtain a lattice approximating the difference equation to the DE in (2). The far-fetched idea of solving a straightforward quadratic equation using the fixed point method has been the root of explosive theoretical developments in the fields of self-similarity, fractals, chaos, and inspired far-reaching consequences in the applications of Mathematics to many scientific fields. Although some skeptics such as Dubois [28] claimed that “the chaos emerging from the so-called chaos map is [not] due to […] fundamental biological properties”, Waliszewski and co-workers and many others insisted on biological examples of fractals, even within our body, such as the growth of cancer tumors, the lungs’ fractal geometry, or the fractal circulatory branching network, where blood vessels branch and branch down to the width of a capillary, see Lorthois and Cassot [29]. Consequently, a sketch of the dynamics of maps tied to extremal distributions has been considered.

Not all growth is sustainable, and in the Gompertz model (Laird [5]; Laird et al. [6]; Norton et al. [7]; Bajzer [8]; Waliszewski et al. [9,10]; Waliszewski and Konarski [11,12,13]; Molski and Konarski, [14]; Tjørve and Tjørve [15])

$$\mathrm{Gompertz}\mathrm{DE}:\frac{\mathrm{d}}{\mathrm{d}t}N\left(t\right)=rN\left(t\right)\left[-ln\left(\frac{N\left(t\right)}{K}\right)\right],$$

(3)

the retroaction factor $-ln\left(\frac{N\left(t\right)}{K}\right)$ is weaker than the retroaction factor $1-\frac{N\left(t\right)}{K}$ in the Verhulst equation. In fact, as ${\displaystyle -ln\left(\frac{N\left(t\right)}{K}\right)=\sum _{n=1}^{\infty }\frac{{\left[1-\frac{N\left(t\right)}{K}\right]}^n}{n}}$, we can consider that the Verhulst model is a first-order approximation of the Gompertz model, partially eliminating moderation exerted by the retroaction factor $-ln\left(\frac{N\left(t\right)}{K}\right)$. As a consequence, population growth is steeper, and the Gompertz model is appropriate for non-sustainable growth, namely cells in cancer tumors.

Observe that the normalized Gompertz function is the Gumbel distribution of maxima; on the other hand, the normalized logistic function is the logistic distribution, the geo-maxima distribution in geometrically thinned sequences, see Section 5.

Arguably, in evolutionary terms, populations tend to grow as much as resources in their ecological niches and competition allow, although it is expectable that the stationary saturation level $\frac{N\left(t\right)}{K}$ is attained at different speeds. Our aim is to describe models having other extremal characteristics, namely finding the exponents in the natural generalizations of the Verhulst and of the Gompertz Equations (2) and (3), namely in the Blumberg [16] hyper-logistic model

$$\mathrm{Blumberg}\mathrm{DE}:\frac{\mathrm{d}}{\mathrm{d}\mathrm{t}}N\left(t\right)=r{\left[N\left(t\right)\right]}^{p-1}{\left[1-\frac{N\left(t\right)}{K}\right]}^{q-1},$$

(4)

in the Turner et al. hyper-Gompertz model

$$\mathrm{Turner}\mathrm{DE}:\frac{\mathrm{d}}{\mathrm{d}t}N\left(t\right)=r{\left[N\left(t\right)\right]}^{p-1}{\left[-ln\left(\frac{N\left(t\right)}{K}\right)\right]}^{\mathsf{\Psi }-1}$$

(5)

and in the modified hyper-Gompertz model considered in Brilhante et al. [20]

$$\mathrm{modified}\mathrm{Turner}\mathrm{DE}:\frac{\mathrm{d}}{\mathrm{d}t}N\left(t\right)=r{\left[1-N\left(t\right)\right]}^{q-1}{\left[-ln\left(1-\frac{N\left(t\right)}{K}\right)\right]}^{\mathsf{\Phi }-1},$$

(6)

whose solutions are proportional to EV probability distributions.

In Equations (2)–(6), the factors ${\left[N\left(t\right)\right]}^{p-1}$ and ${\left[-ln\left(1-\frac{N\left(t\right)}{K}\right)\right]}^{\mathsf{\Phi }-1}$ intervene in growth rate expansion, while the retroaction factors ${\left[1-N\left(t\right)\right]}^{q-1}$ and ${\left[-ln\left(\frac{N\left(t\right)}{K}\right)\right]}^{\mathsf{\Psi }-1}$ account for moderation of growth rate due to resource limitations. This will be analyzed in depth via their connection to the BetaBoop family of functions defined in (1) and the expansions in power laws in Section 2.1.

1.2. BetaBoop Functions

The unifying tools in our research is the family of BetaBoop functions in (1) and the RVs $X_{p,q,\mathsf{\Phi },\mathsf{\Psi }}⌢\mathrm{BetaBoop}(p,q,\mathsf{\Phi },\mathsf{\Psi })$ with probability density function (PDF)

$$f_{p,q,\mathsf{\Phi },\mathsf{\Psi }}\left(x\right)=\frac{G_{p,q,\mathsf{\Phi },\mathsf{\Psi }}\left(x\right)}{{\int }_0^1{G}_{p,q,\mathsf{\Phi },\mathsf{\Psi }}\left(x\right)\mathrm{d}x}$$

(7)

for $p,q,\mathsf{\Phi },\mathsf{\Psi }$ such that ${\int }_0^1{G}_{p,q,\mathsf{\Phi },\mathsf{\Psi }}\left(t\right)\mathrm{d}t<\infty$, see a sufficient condition in Appendix A.

Why BetaBoop? The immediate reason is that the $\mathrm{Beta}(p,q)$ family of PDFs

$$f_{p,q}\left(x\right)=f_{X_{p,q,1,1}}\left(x\right)=\frac{1}{B(p,q)}{x}^{p-1}{(1-x)}^{q-1}{\mathbb{I}}_{(0,1)}\left(x\right),p,q>0,$$

(8)

where $B(p,q)={\int }_0^1{x}^{p-1}{(1-x)}^{q-1}\mathrm{d}x,p,q>0$, is the Beta function or Euler’s integral of the first kind, which can be expressed in terms of the Gamma function or Euler’s integral of the second kind, $\mathsf{\Gamma }\left(\alpha \right)={\int }_0^{\infty }{x}^{\alpha -1}{\mathrm{e}}^{-x}\mathrm{d}x,\alpha >0$, namely $B(p,q)=\frac{\mathsf{\Gamma }\left(p\right)\mathsf{\Gamma }\left(q\right)}{\mathsf{\Gamma }(p+q)}$, is a subfamily of the broader family of BetaBoop PDFs in (7) ($f_{p,q}\left(x\right)=f_{p,q,1,1}\left(x\right)$), and the ultimate and more serious reason is that we love Betty Boop. As vox populi says, “good girls go to Heaven, Betty Boop goes everywhere”, so we expect BetaBoops to lead us everywhere—namely, fractional calculus and generalized monotonicity, convexity and unimodality in Section 3, the logistic and the Gompertz population models and associated maps in Section 4; Richards [30], Blumberg [16], and Turner et al. [18] population models, the hyper-logistic and the hyper-Gompertz maps, and population models connected to EV and geo-EV laws in Section 6. In fact, the logistic parabola $G_{2,2,1,1}\left(x\right)$ is an inspiring connection between the logistic map and BetaBoop PDFs (in this case, the middle $U_{2:3}$ order statistic (OS) of three independent standard Uniform RVs). For general information on the Gamma and on the Beta functions, see Whittaker and Watson [31] and Erdélyi et al. [32].

1.3. Review Organization

An appropriate starting point is the standard uniform RV, $U⌢\mathrm{Uniform}(0,1)$, with PDF $f_U\left(x\right)=f_{1,1,1,1}\left(x\right)=1{\mathbb{I}}_{(0,1)}\left(x\right)$. As it models ideal equiprobability, it does not seem useful for taking decisions on competing risks, and a very naïve view is that it is solely a simple pedagogical example for illustrating basic concepts on continuous univariate distributions theory. A deeper insight, however, shows that it is central in Probability Theory, namely since Borel [33] used the fact that

$${\displaystyle U\stackrel{\mathrm{d}}{=}\sum _{k=1}^{\infty }\frac{X_k}{2^k},\mathrm{where}{X}_k⌢\mathrm{Bernoulli}\left(\frac{1}{2}\right),\mathrm{i}.\mathrm{e}.,X_k=\left\{\begin{array}{cc}0& 1\\ 1/2& 1/2\end{array}\right.,}$$

to provide a rigorous construction of continuous probability, and under very mild conditions, denoting $F_X$, the cumulative distribution function (CDF) of the RV X, $F_X\left(X\right)⌢\mathrm{Uniform}(0,1)$, and therefore, ${\displaystyle F_X^{\leftarrow }\left(U\right)\stackrel{\mathrm{d}}{=}X⌢F_X}$, where $F_X^{\leftarrow }\left(y\right)=inf\{x:F_X\left(x\right)\ge y\}$ is the generalized inverse function of $F_X$. This has been a landmark for X-random numbers generation and computational statistics, and as $1-U⌢\mathrm{Uniform}(0,1)$, antithetic random numbers can be simultaneously generated. When a closed form of the CDF $F_X$ does not exist, as it is the case of the standard $Z⌢\mathrm{Gaussian}(0,1)$ RV, useful workarounds as the Box–Muller algorithm (Box and Muller [34]) transforming a random pair $(U_1,U_2)$ of independent standard uniform RVs in a random pair $(Z_1,Z_2)$ of independent standard Gaussian RVs are used to construct pseudo-random numbers sequences. Observe also that if $X⌢F_X$ and $Y⌢F_Y$, ${\displaystyle F_Y^{\leftarrow }\left[F_X\left(X\right)\right]\stackrel{\mathrm{d}}{=}Y}$. For thorough information on the Uniform RV, consult Johnson et al. [35], pp. 276–321.

In Section 2, we present some complementary information on uniform related RVs, and namely, an overview of the many instances of BetaBoop variables. Further, we show that the $\mathrm{BetaBoop}(p,n,1,1)$, $n\in \mathbb{N}$, and the $\mathrm{BetaBoop}(p,1,1,2)$ RVs are signed mixtures of positive power RVs, causing an interplay of expansion and contraction in their associated nonlinear maps, as indicated below.

Section 3 deals with $G_{1,\alpha ,1,1}$ in fractional calculus and in higher-order monotonicity and convexity (Hadjisavvas et al. [36]), and some consequences in Probability Theory, namely the definition of generalized unimodality and in refinements of Pólya’s class of characteristic functions. Further, $X_{1,\alpha ,1,1}^{1/\xi }$ is the Kumaraswamy RV with positive parameters $\xi$ and $\alpha$ (Kumaraswamy, [37]), whose CDF is $F\left(x\right)=\left\{\begin{array}{cc}0& x<0\\ 1-{(1-x^{\xi })}^{\alpha }& 0\le x<1\\ 1& x\ge 1\end{array}\right.$.

Section 4 is a preliminary presentation of two iconic population growth models, the Verhulst logistic model and the Gompertz model, and of their lattice counterparts, the logistic map $x_{n+1}=r{x}_n(1-x_n)=r{G}_{2,2,1,1}\left(x_n\right)$ and the Gompertz map $x_{n+1}=r{x}_n(-ln{x}_n)=r{G}_{2,1,1,2}\left(x_n\right)$. As $1-x$ is the first-order truncation of $-lnx$, we can view the logistic model as a first-order approximation of the Gompertz model. Hyper-logistic maps $x_{n+1}=r{G}_{p,q,1,1}\left(x_n\right)$ and hyper-Gompertz maps $x_{n+1}=r{G}_{p,1,1,\mathsf{\Psi }}\left(x_n\right)$ are natural extensions, and for appropriate values of the parameters, they are connected with Blumberg [16] and with Turner et al. [17,18] DE population models, respectively.

Section 5 is a brief sketch of extreme value theory (EVT). Aside from the classical independent, identically distributed (IID) framework, we include a brief description of Rachev and Resnick [38] GMS laws. The logistic model is a GMS law in the geometrically thinned framework, and the normalized Gompertz model is the Gumbel CDF of maxima of whole sequences, and the widespread use of those models seems to indicate that dominant species, or cells, have a tendency to grow as much as resources allow. This is further investigated in Section 6, in connection with modeling (linear or logarithmic) of the population growth rate. It is obvious that the Gumbel law dominates the logistic law—as expected when we recall the retroaction factor in the logistic and in the Gompertz maps. The Gompertz model is thus quite successful as a model of immoderate growth of cancerous tumor cells.

The goal of Section 6 is to discuss in depth population models, dealing both with the continuous differential definition and with their recurrence lattice discretization. In Section 6.1, starting from the Tsoularis and Wallace [39] “generalized logistic” encyclopaedic description of population models, we focus on the Blumberg [16] and Turner et al. [17,18] hyper-logistic and hyper-Gompertz models. The logistic paradigm in Section 6.2 serves as guideline for the presentation of those extensions, and the Schwarzian derivative (Singer [40]; Guckenheimer [41]) together with Sharkovskii [42] and Li and Yorke [43] results on period-3 orbits are the basis for the investigation of hyper-logistic and hyper-Gompertz maps. Further, in Section 6.3, we investigate population models tied to other extreme and geo-extreme laws. Finally, in Section 6.4, we discuss Dubois’ [28,44] incursive alternative $x_{n+1}=\frac{r{x}_n}{1+r{x}_n}$ to the classical logistic map, and his claim that May’s [27] “simple mathematical models with very complicated dynamics” result from instabilities of the Euler algorithm.

In the final discussion, Section 7, we indicate some open problems deserving research.

Appendix A establishes the existence of $X_{p,q,\mathsf{\Phi },\mathsf{\Psi }}⌢\mathrm{BetaBoop}(p,q,\mathsf{\Phi },\mathsf{\Psi })$ when $min\{p,q\}$$+max\left\{\mathsf{\Phi },\mathsf{\Psi }\right\}>1$. Appendix B sketches some details on the geo-EVT connected solutions for appropriate parameters of Blumberg’s differential population models (4), and Appendix C analyzes the max-EVT-connected solutions for appropriate parameters of Turner’s differential population models (5) and the min-EVT-connected solutions for appropriate parameters of modified Turner’s differential population models (6).

2. From Uniforms to BetaBoops—An Overview

Let $U⌢\mathrm{Uniform}(0,1)$ be a standard uniform RV and ${\left\{U_n\right\}}_{n\in \mathbb{N}}$ a sequence of IID standard uniform RVs. Obviously, ${\displaystyle 1-U\stackrel{\mathrm{d}}{=}U}$. The PDF of $U^{1/p}$ is $f_{U^{1/p}}\left(x\right)=p{x}^{p-1}{\mathbb{I}}_{(0,1)}\left(x\right)$, and the PDF of $1-U^{1/q}$ is $f_{1-U^{1/q}}\left(x\right)=q{(1-x)}^{q-1}{\mathbb{I}}_{(0,1)}\left(x\right)$. In other words, $U^{1/p}⌢\mathrm{Beta}(p,1)$ and $1-U^{1/q}⌢\mathrm{Beta}(1,q)$. Observe that if in the above statement $p=n$, $n\in \mathbb{N}$, then ${\displaystyle U^{1/n}\stackrel{\mathrm{d}}{=}{U}_{n:n}}$, the maximum of independent $U_1,\dots ,U_n$, and if $q=n$, then ${\displaystyle 1-U^{1/n}\stackrel{\mathrm{d}}{=}{U}_{1:n}}$, the minimum of independent $U_1,\dots ,U_n$. More generally, the k-th ascending OS, $U_{k:n}$, of the independent $U_1,\dots ,U_n$ is a $\mathrm{Beta}(k,n+1-k)$ RV (k and $n+1-k$ are the ascending and the descending ranks, respectively).

Mellin transforms, ${\mathcal{M}}_X\left(s\right)=\mathbb{E}\left[X^s\right]={\int }_0^{\infty }{x}^s\mathrm{d}{F}_X\left(x\right)$, $X\ge 0$, are very useful in the multiplicative algebra of non-negative independent RVs, since if $X\ge 0$ and $Y\ge 0$ are independent, then ${\mathcal{M}}_{XY}\left(s\right)={\mathcal{M}}_X\left(s\right){\mathcal{M}}_Y\left(s\right)$, for general information, see Zolotarev [45]. The Mellin transform of $U⌢\mathrm{Uniform}(0,1)$, ${\mathcal{M}}_U\left(s\right)=\mathbb{E}\left[U^s\right]$, is ${\int }_0^1{x}^s\mathrm{d}x=\frac{1}{1+s},\Re \left(s\right)>-1$. As ${\int }_0^1{x}^s\frac{{\left(-lnx\right)}^{n-1}}{\mathsf{\Gamma }\left(n\right)}\mathrm{d}x={\left(\frac{1}{1+s}\right)}^n$, $\Re \left(s\right)>-1$, we conclude that the PDF of the product of n independent standard uniform RVs, ${\prod }_{k=1}^n{U}_k$, is $f_{{\prod }_{k=1}^n{U}_k}\left(x\right)=\frac{{\left(-lnx\right)}^{n-1}}{\mathsf{\Gamma }\left(n\right)}{\mathbb{I}}_{(0,1)}\left(x\right)$. Observe also that ${\displaystyle {\prod }_{k=1}^n{U}_k\stackrel{\mathrm{d}}{=}{\mathcal{X}}_n}$, where ${\mathcal{X}}_0⌢\mathrm{Uniform}(0,1)$ and ${{\mathcal{X}}_{k+1}}_{|{\mathcal{X}}_k}⌢\mathrm{Uniform}(0,{\mathcal{X}}_k)$, $k=1,\dots ,n-1$, and that ${\displaystyle 1-{\prod }_{k=1}^n(1-U_k)\stackrel{\mathrm{d}}{=}{\mathcal{Y}}_n}$, where ${\mathcal{Y}}_0⌢\mathrm{Uniform}(0,1)$ and ${{\mathcal{Y}}_{k+1}}_{|{\mathcal{Y}}_k}⌢\mathrm{Uniform}({\mathcal{Y}}_k,1)$, $k=1,\dots ,n-1$.

As ${\mathcal{M}}_{U^{1/p}}\left(s\right)=\frac{p}{p+s},\Re \left(s\right)\ge -p$, and

$${\int }_0^1{x}^s\frac{p^n}{\mathsf{\Gamma }\left(n\right)}{x}^{p-1}{\left(-lnx\right)}^{n-1}\mathrm{d}x={\left(\frac{p}{p+s}\right)}^n,\Re \left(s\right)\ge -p,$$

the PDF of the product of n independent Beta $(p,1)$ RVs is

$$f_{{\prod }_{k=1}^n{U}_k^{1/p}}\left(x\right)=\frac{p^n}{\mathsf{\Gamma }\left(n\right)}{x}^{p-1}{\left(-lnx\right)}^{n-1}{\mathbb{I}}_{(0,1)}\left(x\right).$$

More generally, $f\left(x\right)=\frac{p^{\mathsf{\Psi }}}{\mathsf{\Gamma }\left(\mathsf{\Psi }\right)}{x}^{p-1}{(-lnx)}^{\mathsf{\Psi }-1}{\mathbb{I}}_{(0,1)}\left(x\right),p,\mathsf{\Psi }>0$, is a PDF.

Hence, a huge family of important random models may be defined starting from a sequence of n IID standard $\mathrm{Uniform}(0,1)$ RVs, $U_1,U_2,\dots ,U_n$, as exemplified below, where as usual, $U_{k:n}$ denotes the k-th ascending OS:

$$\underset{\frac{{\left(-lnx\right)}^{n-1}}{\mathsf{\Gamma }\left(n\right)}{\mathbb{I}}_{(0,1)}\left(x\right)}{\underbrace{\prod _{k=1}^n{U}_k}}\prec \underset{n{(1-x)}^{n-1}{\mathbb{I}}_{(0,1)}\left(x\right)}{\underbrace{U_{1:n}}}\prec \underset{\frac{x^{k-1}{(1-x)}^{n-k}}{B(k,n+1-k)}{\mathbb{I}}_{(0,1)}\left(x\right)}{\underbrace{U_{k:n}}}\prec \underset{n{x}^{n-1}{\mathbb{I}}_{(0,1)}\left(x\right)}{\underbrace{U_{n:n}}}\prec \underset{\frac{{\left[-ln(1-x)\right]}^{n-1}}{\mathsf{\Gamma }\left(n\right)}{\mathbb{I}}_{(0,1)}\left(x\right)}{\underbrace{1-\prod _{k=1}^n(1-U_k)}}.$$

Any of the PDFs in the underbrackets is readily extended considering non-integer positive exponents. Thus, whenever the BetaBoop function $G_{p,q,\mathsf{\Phi },\mathsf{\Psi }}\left(x\right)$, in (1), is integrable—a sufficient condition is $1<min\left\{\mathsf{\Phi },\mathsf{\Psi }\right\}+min\{p,q\}$, see Appendix A, $f_{p,q,\mathsf{\Phi },\mathsf{\Psi }}\left(x\right)=\frac{G_{p,q,\mathsf{\Phi },\mathsf{\Psi }}\left(x\right)}{{\int }_0^1{G}_{p,q,\mathsf{\Phi },\mathsf{\Psi }}\left(x\right)\mathrm{d}x}$, in (7), is a PDF of a $X_{p,q,\mathsf{\Phi },\mathsf{\Psi }}⌢\mathrm{BetaBoop}(p,q,\mathsf{\Phi },\mathsf{\Psi })$ RV.

So, the PDFs of products and of OSs of IID standard uniforms in the above inequalities underbrackets are the PDFs of, respectively, $\mathrm{BetaBoop}(1,1,1,n)$, $\mathrm{BetaBoop}(1,n,1,1)$, $\mathrm{BetaBoop}(k,n+1-k,1,1)$, $\mathrm{BetaBoop}(n,1,1,1)$, and $\mathrm{BetaBoop}(1,1,n,1)$, and the corresponding positive exponents extensions are:

1.

$\left\{f_{1,1,1,\mathsf{\Psi }}\left(x\right)=\frac{1}{\mathsf{\Gamma }\left(\mathsf{\Psi }\right)}{(-lnx)}^{\mathsf{\Psi }-1}{\mathbb{I}}_{(0,1)}\left(x\right),\mathsf{\Psi }>0\right\}$, generalized logarithmic PDFs, the $p=1$ subfamily of power-logarithmic PDFs,

$$\left\{f_{p,1,1,\mathsf{\Psi }}\left(x\right)=\frac{p^{\mathsf{\Psi }}}{\mathsf{\Gamma }\left(\mathsf{\Psi }\right)}{x}^{p-1}{\left(-lnx\right)}^{\mathsf{\Psi }-1}{\mathbb{I}}_{(0,1)}\left(x\right),p,\mathsf{\Psi }>0\right\};$$

2.

$\left\{f_{1,q,1,1}\left(x\right)=q{(1-x)}^{q-1}{\mathbb{I}}_{(0,1)}\left(x\right),q>0\right\}$, the subfamily of $\mathrm{Beta}(1,q)$ PDFs;

3.

$\left\{f_{p,q,1,1}\left(x\right)=\frac{1}{B(p,q)}{x}^{p-1}{(1-x)}^{q-1}{\mathbb{I}}_{(0,1)}\left(x\right),p,q>0\right\}$, the $\mathrm{Beta}(p,q)$ PDFs;

4.

$\left\{f_{p,1,1,1}\left(x\right)=p{x}^{p-1}{\mathbb{I}}_{(0,1)}\left(x\right),p>0\right\}$, the positive power PDFs, the subfamily of $\mathrm{Beta}(p,1)$ PDFs;

5.

$\left\{f_{1,1,\mathsf{\Phi },1}\left(x\right)=\frac{1}{\mathsf{\Gamma }\left(\mathsf{\Phi }\right)}{\left[-ln(1-x)\right]}^{\mathsf{\Phi }-1}{\mathbb{I}}_{(0,1)}\left(x\right),\mathsf{\Phi }>0\right\}$, reverted generalized logarithmic PDFs, the $q=1$ subfamily of

$$\left\{f_{1,q,\mathsf{\Phi },1}\left(x\right)=\frac{q^{\mathsf{\Phi }}}{\mathsf{\Gamma }\left(\mathsf{\Phi }\right)}{(1-x)}^{q-1}{\left[-ln(1-x)\right]}^{\mathsf{\Phi }-1}{\mathbb{I}}_{(0,1)}\left(x\right),q,\mathsf{\Phi }>0\right\}.$$

It is obvious that the factors $x^{p-1}(=G_{p,1,1,1}\left(x\right))$ and ${\left[-ln(1-x)\right]}^{\mathsf{\Phi }-1}(=G_{1,1,\mathsf{\Phi },1}\left(x\right))$ are left-skewed, with $x^{p-1}$ dominating ${\left[-ln(1-x)\right]}^{\mathsf{\Phi }-1}$, while the factors ${(1-x)}^{q-1}$ $(=G_{1,q,1,1}\left(x\right))$ and ${\left(-lnx\right)}^{\mathsf{\Psi }-1}(=\left(G_{1,1,1,\mathsf{\Psi }}\left(x\right)\right)$ are right-skewed, with ${(1-x)}^{q-1}$ dominating ${\left(-lnx\right)}^{\mathsf{\Psi }-1}$. In Figure 1, we plot some BetaBoop PDFs, where we observe symmetries with respect to $x=\frac{1}{2}$ due to $f_{p,q,\mathsf{\Phi },\mathsf{\Psi }}\left(x\right)=f_{q,p,\mathsf{\Psi },\mathsf{\Phi }}(1-x)$.

2.1. BetaBoop $(p,q,1,1)$, BetaBoop $(2,1,1,2)$ and Power Laws

The $\mathrm{BetaBoop}(2,2,1,1)$ PDF

$$f_{2,2,1,1,}\left(x\right)=6x(1-x){\mathbb{I}}_{(0,1)}\left(x\right)=3\times \underset{\mathrm{Beta}(2,1)}{\underbrace{2x{\mathbb{I}}_{(0,1)}\left(x\right)}}-2\times \underset{\mathrm{Beta}(3,1)}{\underbrace{3x^2{\mathbb{I}}_{(0,1)}\left(x\right)}}$$

is a signed mixture of the $\mathrm{Beta}(2,1)$, with weight 3, and the $\mathrm{Beta}(3,1)$, with weight $-2$. Similar results hold for $\mathrm{Beta}(p,n)$ RVs:

$$f_{p,n,1,1}\left(x\right)=\sum _{j=0}^{n-1}{(-1)}^j\left(\genfrac{}{}{0pt}{}{n-1}{j}\right)\frac{\mathsf{\Gamma }(p+n)}{(p+j)(n-1)!\mathsf{\Gamma }\left(p\right)}\underset{\mathrm{Beta}(p+j,1)}{\underbrace{(p+j)x^{p+j-1}{\mathbb{I}}_{(0,1)}\left(x\right)}}.$$

So, the $\mathrm{BetaBoop}(p,q,1,1)$ PDF is a signed mixture of power laws with natural exponents, each positive term contributing to the expansion being followed by a negative term moderating growth.

On the other hand, observe that the $\mathrm{BetaBoop}(2,1,1,2)$ PDF can be written as

$$\begin{array}{cc}\hfill f_{2,1,1,2}\left(x\right)& =4x(-lnx){\mathbb{I}}_{(0,1)}\left(x\right)=4x\sum _{k=1}^{\infty }\frac{{(1-x)}^k}{k}{\mathbb{I}}_{(0,1)}\left(x\right)=\hfill \\ \hfill & =\sum _{k=1}^{\infty }\underset{w_k}{\underbrace{\frac{4}{k(k+1)(k+2)}}}\underset{\mathrm{Beta}(2,k+1)}{\underbrace{\frac{1}{B(2,k+1)}x{(1-x)}^k{\mathbb{I}}_{(0,1)}\left(x\right)}},\hfill \end{array}$$

and, therefore, the $\mathrm{BetaBoop}(2,1,1,2)$ is a mixture of $\mathrm{BetaBoop}(2,k+1,1,1)$ with weights $w_k=\frac{4}{k(k+1)(k+2)}$. Therefore,

$$\begin{array}{cc}\hfill f_{2,1,1,2}\left(x\right)& =\sum _{k=1}^{\infty }\frac{4}{k(k+1)(k+2)}\frac{1}{B(2,k+1)}x\sum _{j=0}^k{(-1)}^j\left(\genfrac{}{}{0pt}{}{k}{j}\right)x^j{\mathbb{I}}_{(0,1)}\left(x\right)=\hfill \\ & =\sum _{k=1}^{\infty }\underset{w_k^{*}}{\underbrace{\sum _{j=0}^k{(-1)}^j\frac{4}{k(j+2)}\left(\genfrac{}{}{0pt}{}{k}{j}\right)}}\underset{\mathrm{Beta}(j+2,1)}{\underbrace{(j+2)x^{j+1}{\mathbb{I}}_{(0,1)}\left(x\right)}}.\hfill \end{array}$$

So, the $\mathrm{BetaBoop}(2,1,1,2)$ PDF is a mixture of power laws with natural exponents, positive weights contributing to expansion and negative weights moderating growth; similarly, the $\mathrm{BetaBoop}(p,1,1,2)$ PDF can be re-expressed as a mixture of power laws,

$$f_{p,1,1,2}\left(x\right)=\sum _{k=1}^{\infty }\sum _{j=0}^k{(-1)}^j\frac{p^2}{k(p+j)}\left(\genfrac{}{}{0pt}{}{k}{j}\right)\underset{\mathrm{Beta}(p+j,1)}{\underbrace{(p+j)x^{p+j-1}{\mathbb{I}}_{(0,1)}\left(x\right)}}.$$

This is indeed relevant for understanding the role of the expansion/contraction effects of the couple of factors in the BetaBoop-connected population growth models that we shall investigate in Section 6 as extensions of the Verhulst population growth model.

Power laws are at the root of scaling and self-similarity, and hence intervene in many areas of the theory of fractals and chaotic dynamics, as explained in detail by Schroeder [46], ch. 12. Karamata’s [47] regular variation theory, a useful tool in many branches of Mathematics (Bingham et al. [48]), essentially uses properties of power laws.

Feller [49] observed that Doeblin’s [50] and Gnedenko’s [51,52] theory of the domains of attraction of the additive stable laws could be simply rephrased in Karamata’s framework. Since then, slow variation and regular variation became the election tool for characterizing domains of attraction of stable limit laws, namely in de Haan [53] for a complete characterization of domains of attraction for stable EV laws, in Bingham [54] for a solution of a more general problem of domains of attraction of generalized convolution algebras, and in Kozubowski and Rachev [55] for a characterization of geometric-stable domains of attraction. As far as we are aware, no one pointed out that weak convergence of types problems in Probability Theory are connected to the renormalization group theory in Mathematical Physics, see Stanley [56], and that the appropriateness of the regular variation theory to deal with domains of attraction is a side-effect of the self-similarity resulting from power laws and scaling.

The fact that logistic and hyper-logistic maps and Gompertz and hyper-Gompertz maps are signed measures of power laws is obviously connected to the various forms of equilibrium appropriate to describe diverse population dynamics, from extinction, stagnation, sustainable growth, to the various forms of extreme growth that we shall discuss later.

3. BetaBoop(1,q,1,1), Fractional Calculus, Generalized Monotonicity and Convexity, and Applications in Probability Theory

3.1. Fractional Calculus and Some Applications in Probability Theory

The function $G_{1,\alpha ,1,1}$ plays an important role in fractional calculus, namely in the Riemann–Liouville fractional integral. Starting from the Cauchy’s iterated integral (Goursat [57]), using $G_{1,n,1,1}\left(t\right)={(1-t)}^{n-1}{\mathbb{I}}_{(0,1)}\left(t\right)$,

$$\begin{array}{cc}\hfill {\mathrm{I}}_{x_0}^n\left[f\left(x\right)\right]& ={\int }_{x_0}^x{\int }_{x_0}^{t_n}\cdots {\int }_{x_0}^{t_2}f\left(t_1\right)\mathrm{d}{t}_1\mathrm{d}{t}_2\dots \mathrm{d}{t}_n=\hfill \\ \hfill & {\displaystyle ={\int }_{x_0}^{x}f\left(t\right)\frac{{(x-t)}^{n-1}}{(n-1)!}\mathrm{d}t=\frac{x^n}{\mathsf{\Gamma }\left(n\right)}{\int }_{\frac{x_0}{x}}^1{(1-t)}^{n-1}f\left(tx\right)\mathrm{d}t.}\hfill \end{array}$$

Replacing the integer index n by $\alpha >0$, we obtain

$${\displaystyle {\mathrm{I}}_{x_0}^{\alpha }\left[f\left(x\right)\right]={\int }_{x_0}^x\frac{{(x-t)}^{\alpha -1}}{\mathsf{\Gamma }\left(\alpha \right)}f\left(t\right)\mathrm{d}t=\frac{x^{\alpha }}{\mathsf{\Gamma }\left(\alpha \right)}{\int }_{\frac{x_0}{x}}^1{(1-t)}^{\alpha -1}f\left(tx\right)\mathrm{d}t,}$$

a suitable generalization of the integral operator to non-integer order, since

1.

${\mathrm{I}}_{x_0}^{\alpha }\left[f\left(x\right)+g\left(x\right)\right]={\mathrm{I}}_{x_0}^{\alpha }\left[f\left(x\right)\right]+{\mathrm{I}}_{x_0}^{\alpha }\left[g\left(x\right)\right]$;

2.

${\mathrm{I}}_{x_0}^{\alpha +\beta }\left[f\left(x\right)\right]={\mathrm{I}}_{x_0}^{\alpha }\left[{\mathrm{I}}_{x_0}^{\beta }\left[f\left(x\right)\right]\right]$;

3.

$\underset{\alpha \downarrow 0}{lim}{\mathrm{I}}_{x_0}^{\alpha }\left[f\left(x\right)\right]=f\left(x\right)$.

From now on, we consider $x_0=0$, and we shall write simply

$${\displaystyle {\mathrm{I}}^{\alpha }\left[f\left(z\right)\right]=\frac{z^{\alpha }}{\mathsf{\Gamma }\left(\alpha \right)}{\int }_0^1{(1-t)}^{\alpha -1}f\left(tz\right)\mathrm{d}t=\frac{z^{\alpha }}{\mathsf{\Gamma }\left(\alpha \right)}{\int }_0^1{G}_{1,\alpha ,1,1}\left(t\right)f\left(tz\right)\mathrm{d}t,}$$

which is generally called the Riemann–Liouville fractional integral, although none of them used exactly this form for the fractional integral operator, see Ross [58] and Liouville [59], and interesting points of view on the history of fractional integrals in Dugowson [60].

The non-integer differential operator of order $\alpha$ is defined by

$$D^{\alpha }\left[f\left(x\right)\right]={\mathrm{I}}^{-\alpha }\left[f\left(x\right)\right]=\frac{{\mathrm{d}}^n}{\mathrm{d}{x}^n}{\mathrm{I}}^{n-\alpha }\left[f\left(x\right)\right],$$

where n is any integer greater than $\alpha$.

The definition of integrals and derivatives of non-integer order may be made in several ways, each of them having its own drawbacks, see, for instance, the use of the Prabhakar fractional integral and of the Liouville–Caputo fractional derivative in Area and Nieto [61]. For an excellent review of the existing definitions and of their limitations, see Lavoie et al. [62] and Luchko [63].

Fractional Calculus is a useful tool in several areas of Probability Theory (Pestana [64]; Gomes and Pestana [65]). The Riemann–Liouville integral played an important role in the theory of the addition of independent RVs. Lévy [66] used the Riemann–Liouville integral to establish the existence of non-normal additive stable laws, i.e., additive stable laws with characteristic exponent $\alpha \in (0,2)$. Feller [67] used the extension of M. Riesz’ potentials and the Riemann–Liouville integral to obtain an asymptotic series for the stable distribution functions, see also Wintner [68]. Wolfe [69] contains an overview of theory and applications of fractional moments.

Some useful references on Fractional Calculus: Oldham and Spanier [70], Samko et al. [71], Daftardar-Gejji [72], Katugampola [73], Herrmann [74], and Luchko [63].

3.2. Higher-Order Monotone Functions, Generalized Convexity, and Applications in Probability Theory

Denote $M_n$ the class of functions $f\ge 0$ and bounded with support $A\subseteq [0,\infty ]$; f is monotone of order $n,n=0,1,\dots \in [0,\infty ]$ if $f\in {\mathcal{C}}^n\left(A\right)$ and its derivative of order n, $f^{\left(n\right)}$ is monotone.

The boundedness of f implies that ${\displaystyle f^{\left(n\right)}\left(x\right)\underset{x\to \infty }{⟶}0}$, and hence, $f^{\left(n\right)}$ is either non-negative and non-increasing or non-positive and non-decreasing. $M_0$ is the class of monotone and bounded non-negative functions, and obviously, $M_n\subset M_{n-1}$. Functions f with support $A\in [-\infty ,0]$ are monotone of order n if $f(-x)$ is monotone of order n.

If $f\in M_n$, then $\left|f^{\left(n\right)}\right|$ is non-negative and non-increasing and ${(-1)}^k{f}^{\left(k\right)}$ is non-negative, non-increasing and convex for $k=0,1,\dots ,n-1$. Therefore, $M_{\infty }$ is the class of Bernstein [75] completely monotone functions, i.e., Laplace transforms of a bounded and right-continuous non-decreasing function (Feller [76], ch. XIII).

Denote $D_n,n\ge 1,$ the class of CDFs F such that $F^{\prime }\left(x\right){\mathbb{I}}_{[0,\infty ]}\left(x\right)$ and $F^{\prime }\left(x\right){\mathbb{I}}_{[-\infty ,0]}\left(x\right)$ are monotone of order $n-1$. The extreme points of the convex set $D_n$ are the degenerate distributions at $-\infty$, at 0 and at ∞, and the non-degenerate extreme points are $F_{n,a}$ defined by

$$F_{n,a}^{\prime }\left(x\right)=\frac{n}{{\left|a\right|}^n}{\left(\right|a|-x\mathrm{sign}\left(x\right))}^{n-1}{\mathbb{I}}_{(min(0,a),max(0,a\left)\right)}\left(x\right),a\in \mathbb{R}\setminus \left\{0\right\},$$

for details, see Pestana and Mendonça [77].

From the Choquet extreme point integral representation theory (Choquet [78]; Phelps [79]),

$$f\in M_n⟺f\left(x\right)={\int }_0^{1/x}{(1-xt)}^n\mathrm{d}G\left(t\right),x>0,$$

(9)

with G bounded below and non-decreasing, and consequently,

$$F_Y\in D_n⟺Y\stackrel{d}{=}(1-U^{1/n})Y_n,$$

with $U⌢\mathrm{Uniform}(0,1)$ and $Y_n$ and $1-U^{1/n}$ independent RVs.

For simplicity, assume that $G\left(0\right)=0$. From (9),

$$x^{n}f\left(\frac{1}{x}\right)={\int }_0^x{(x-t)}^n\mathrm{d}G\left(t\right)⟹\frac{{\mathrm{d}}^n}{\mathrm{d}{x}^n}\left[x^{n}f\left(\frac{1}{x}\right)\right]=n!G\left(x\right).$$

On the other hand, if we assume that $\frac{{\mathrm{d}}^n}{\mathrm{d}{x}^n}\left[x^{n}f\left(\frac{1}{x}\right)\right]$ is non-negative and non-decreasing, defining $G\left(x\right)=\frac{1}{n!}\frac{{\mathrm{d}}^n}{\mathrm{d}{x}^n}\left[x^{n}f\left(\frac{1}{x}\right)\right]{\mathbb{I}}_{(0,\infty )}\left(x\right)$, then $h\left(x\right)={\int }_0^{1/x}{(1-xt)}^n\mathrm{d}G\left(t\right)\in M_n$ and $\frac{1}{n!}\frac{{\mathrm{d}}^n}{\mathrm{d}{x}^n}\left[x^{n}h\left(\frac{1}{x}\right)\right]=\frac{1}{n!}\frac{{\mathrm{d}}^n}{\mathrm{d}{x}^n}\left[x^{n}f\left(\frac{1}{x}\right)\right]$ implying that $x^{n}h\left(\frac{1}{x}\right)=x^{n}f\left(\frac{1}{x}\right)+{\sum }_{k=0}^n{\alpha }_k{x}^k$.

Thus, assuming that ${\displaystyle \underset{x\to \infty }{lim}f\left(x\right)<\infty }$ and ${\displaystyle \underset{x\to \infty }{lim}h\left(x\right)<\infty }$, we conclude that $h\equiv f$.

So, f with support $A\subseteq (0,\infty )$ is monotone of order n iff ${\displaystyle 0\le \underset{x\to \infty }{lim}f\left(x\right)<\infty }$ and $\frac{{\mathrm{d}}^n}{\mathrm{d}{x}^n}\left[x^{n}f\left(\frac{1}{x}\right)\right]$ is bounded below and non-decreasing in $(0,\infty )$.

All this can almost trivially be extended using fractional $\nu$ instead of n: A function f is said to be monotone of order $\nu ,\nu >0$, iff ${\displaystyle 0\le \underset{x\to \infty }{lim}f\left(x\right)<\infty }$ and the fractional derivative ${\mathrm{I}}^{-\nu }\left[x^{\nu }f\left(\frac{1}{x}\right)\right]$ is bounded below and non-decreasing.

Equivalently, $f\in M_{\nu }$ iff $f\left(x\right)={\int }_0^{1/x}{(1-xt)}^{\nu }\mathrm{d}G\left(t\right)$ with G bounded below and non-decreasing in $(0,\infty )$.

Defining the class of fractional unimodal distributions $D_{\nu }$ in analogy with $D_n$,

$$F_Y\in D_{\nu },\nu =1,2,\dots ⟺Y\stackrel{\mathrm{d}}{=}(1-U^{1/\nu })Y_{\nu },$$

where $U⌢\mathrm{Uniform}(0,1)$ and $Y_{\nu }$ and $1-U^{1/\nu }$ are independent.

Additional details can be found in Pestana and Mendonça [77]. Further, observe that if $0<\nu <\mu$, then $M_{\mu }\subset M_{\nu }$ and $D_{\mu }\subset D_{\nu }$, since simple algebra with Mellin transforms

$$\frac{\mathsf{\Gamma }(1+s)\mathsf{\Gamma }(1+\mu )}{\mathsf{\Gamma }\left(1\right)\mathsf{\Gamma }(1+\mu +s)}=\frac{\mathsf{\Gamma }(1+s)\mathsf{\Gamma }(1+\nu )}{\mathsf{\Gamma }\left(1\right)\mathsf{\Gamma }(1+\nu +s)}\frac{\mathsf{\Gamma }(1+\nu +s)\mathsf{\Gamma }(1+\mu )}{\mathsf{\Gamma }(1+\nu )\mathsf{\Gamma }(1+\mu +s)}$$

shows that $1-U^{1/\mu }=(1-U^{1/\nu })X_{1+\nu ,\mu -\nu ,1,1}$.

The classical result on unimodality—the RV X is unimodal with mode at 0 iff $X=UY$ with U and Y independent (Khinchine [80])—is the special case $\nu =1$ of the above generalized unimodality $Y=(1-U^{1/\nu })Y_{\nu }$, using $1-U^{1/\nu }⌢\mathrm{BetaBoop}(1,\nu ,1,1)$ (Pestana [64,81]). As $X_{1,\nu ,1,1}$ has a mode at 0, generalized unimodality is a natural extension of unimodality. Compare with the Olshen and Savage [82] star-shaped unimodality $X=U^{1/\alpha }Y$ using $X_{\alpha ,1,1,1}$ with a mode at 1.

The function $G_{1,q,1,1}$ has been used by Pestana [64] and Pestana and Mendonça [77] to investigate higher-order monotone functions in Probability Theory, and also on refinements of Pólya’s characteristic functions (Pestana [83]), namely Sakovič [84] concave characteristic functions.

Further information on convexity: Roberts and Varberg [85] is a good read, and is the source for the following quotation of Jensen: “It seems to me that the notion of convex function is just as fundamental as positive function or increasing function. If I am not mistaken in this, the notion ought to find its place in elementary expositions of the theory of real functions.” Hadjisavvas et al. [36] cover a broad diversity of topics, and Vivas and Hernández’s [86] preprint presents recent developments. For advanced recent work, consult Sitthiwirattham et al. [87], Sahoo et al. [88], and references therein.

3.3. Fractional Powers of BetaBoop Random Variables

Observe also in what concerns fractional powers of $X_{1,\alpha ,1,1}⌢\mathrm{BetaBoop}(1,\alpha ,1,1)$, that the PDF of ${\mathcal{K}}_{\xi ,\alpha }=X_{1,\alpha ,1,1}^{1/\xi },\alpha ,\xi >0$, is $f_{{\mathcal{K}}_{\xi ,\alpha }}\left(x\right)=\xi \alpha x^{\xi -1}{\left(1-x^{\xi }\right)}^{\alpha -1}{\mathbb{I}}_{(0,1)}\left(x\right)$. So, $X_{1,\alpha ,1,1}^{1/\xi }$ is the Kumaraswamy RV with parameters $\xi$ and $\alpha$ (Kumaraswamy, [37]; Jones, [89]). Obviously, the PDFs of Kumaraswamy’s generalized RVs $X_{p,q,1,1}^{1/\xi }$ or even $X_{p,q,\mathsf{\Phi },\mathsf{\Psi }}^{1/\xi }$ are readily obtained:

Let $Y=X_{p,q,\mathsf{\Phi },\mathsf{\Psi }}^{1/\xi }$, with $p,q,\mathsf{\Phi },\mathsf{\Psi },\xi >0$. The PDF of Y is given by

$$f_Y\left(y\right)=f_{p,q,\mathsf{\Phi },\mathsf{\Psi }}\left(y^{\xi }\right)\xi y^{\xi -1}{\mathbb{I}}_{(0,1)}\left(y\right).$$

Therefore,

$$f_Y\left(y\right)=C{\xi }^{\mathsf{\Psi }}{y}^{\xi p-1}{\left(1-y^{\xi }\right)}^{q-1}{\left[-ln\left(1-y^{\xi }\right)\right]}^{\mathsf{\Phi }-1}{\left(-lny\right)}^{\mathsf{\Psi }-1}{\mathbb{I}}_{(0,1)}\left(y\right),$$

where $C=\frac{1}{{\int }_0^1{G}_{p,q,\mathsf{\Phi },\mathsf{\Psi }}\left(x\right)\mathrm{d}x}$.

It is worth noticing that ${\displaystyle X_{p,1,1,\mathsf{\Psi }}^{1/\xi }\stackrel{\mathrm{d}}{=}{X}_{p\xi ,1,1,\mathsf{\Psi }}}$, and in particular, ${\displaystyle X_{p,1,1,1}^{1/\xi }\stackrel{\mathrm{d}}{=}{X}_{p\xi ,1,1,1}}$ and ${\displaystyle X_{1,1,1,\mathsf{\Psi }}^{1/\xi }\stackrel{\mathrm{d}}{=}{X}_{\xi ,1,1,\mathsf{\Psi }}}$.

4. Logistic Growth, Gompertz Growth, and Extensions of the Logistic Map

As implied in the observation on $f_{2,2,1,1}=f_{U_{2:3}}$ and the logistic map, there is a strong connection tying $r{G}_{2,2,1,1}:(0,1]\to (0,1],r\in (0,4]$, to the Verhulst population growth model defined in Equation (2)

$$\frac{\mathrm{d}}{\mathrm{d}\mathrm{t}}N\left(t\right)=rN\left(t\right)\left[1-\frac{N\left(t\right)}{K}\right],$$

with logistic solution $N\left(t\right)=\frac{K}{1+\frac{K-N\left(0\right)}{N\left(0\right)}{\mathrm{e}}^{-rt}}$ (K is the saturation level or carrying capacity, r is the intrinsic growth rate, and in the Verhulst model, the relative growth $\frac{\frac{\mathrm{d}}{\mathrm{d}t}N\left(t\right)}{N\left(t\right)}$ decreases linearly with the population size, which will ultimately reach the carrying capacity K). Tjørve and Tjørve [90] provide a thorough discussion of the logistic model as a member of the Richards [30] family of growth models, with a useful overview of parameter interpretation.

The quadratic logistic map

$$x_{n+1}=r{x}_n(1-x_n),r\in (0,4],$$

(10)

is the classical discrete version obtained from the Verhulst differential population growth model $\frac{\mathrm{d}}{\mathrm{d}t}\alpha \left(t\right)=r\alpha \left(t\right)\left[1-\alpha \left(t\right)\right],r\in (0,4]$, $\alpha \left(t\right)=\frac{N\left(t\right)}{K}\in [0,1]$, by using Euler’s algorithm (for simplicity, we use the same r, but this is a shortcut for: “$\alpha (t+\Delta t)=\alpha \left(t\right)+\Delta t[r\alpha \left(t\right)-r{\alpha }^2\left(t\right)]$, where the time step $\Delta t$ may be taken equal to 1 in rescaling the growth rate r, see Equations (21) and (22) in Section 6.2”). It plays a prominent role in what concerns Chaos Theory since the publication of the pathbreaking analysis of May [27], interesting overviews are given in Schroeder [46], ch. XII, or in Peitgen et al. [91], ch. 10–11.

Observe that for $x\in (0,1)$, we have $\frac{\mathrm{d}}{\mathrm{d}x}{\sum }_{k=1}^{\infty }\frac{x^k}{k}=\frac{1}{1-x}$ and $\frac{\mathrm{d}}{\mathrm{d}x}{\sum }_{k=1}^{\infty }\frac{{(1-x)}^k}{k}=-\frac{1}{x}$, so $-ln(1-x)={\sum }_{k=1}^{\infty }\frac{x^k}{k}$ and $-lnx={\sum }_{k=1}^{\infty }\frac{{(1-x)}^k}{k}$. So, we may regard the logistic map $x_{n+1}=r{x}_n(1-x_n)$ as a first-order approximation of the “Gompertz map” $x_{n+1}=r{x}_n(-ln{x}_n)$ (or of the map $x_{n+1}=r(1-x_n)(-ln(1-x_n))$, or even of $x_{n+1}=r(-ln(1-x_n))(-ln{x}_n)$).

The Gompertz map $r{G}_{2,1,1,2}:(0,1]\to (0,1],r\in (0,\mathrm{e}]$ is the lattice version $x_{n+1}=r{x}_n(-ln{x}_n)$ of the DE $\frac{\mathrm{d}}{\mathrm{d}\mathrm{t}}\alpha \left(t\right)=r\alpha \left(t\right)\left[-ln\left(\alpha \left(t\right)\right)\right]$ (Laird [5]; Laird et al. [6]; Waliszewski and Konarski [12,13]), obtained using the saturation level $\alpha \left(t\right)=\frac{N\left(t\right)}{K}$ in the Gompertz population model defined in Equation (3)

$$\frac{\mathrm{d}}{\mathrm{d}\mathrm{t}}N\left(t\right)=rN\left(t\right)\left[-ln\left(\frac{N\left(t\right)}{K}\right)\right]⟹N\left(t\right)=Kexp(-{\mathrm{e}}^{-rt+C}),$$

where ${\mathrm{e}}^C=-ln\frac{N\left(0\right)}{K}$, proportional to the Gumbel CDF with location and scale parameters $-C/r$ and $1/r$, respectively, see details in Appendix B. Tjørve and Tjørve [15] have a coordinated presentation of useful reparameterizations of the Gompertz model in growth modeling.

Blumberg [16] introduced a hyper-logistic growth function, solution of the modified Verhulst DE, $\frac{\mathrm{d}}{\mathrm{d}\mathrm{t}}N\left(t\right)=r{\left[N\left(t\right)\right]}^{\alpha }{\left[1-\frac{N\left(t\right)}{K}\right]}^{\gamma }$, which we reparameterize as Equation (4)

$$\frac{\mathrm{d}}{\mathrm{d}\mathrm{t}}N\left(t\right)=r{\left[N\left(t\right)\right]}^{p-1}{\left[1-\frac{N\left(t\right)}{K}\right]}^{q-1},$$

for direct association with the BetaBoop functions.

The corresponding $r{G}_{p,q,1,1}:(0,1]\to (0,1],r\in \left(0,{\left(\frac{p-1}{p+q-2}\right)}^{1-p}{\left(\frac{q-1}{p+q-2}\right)}^{1-q}\right)$, hyper-logistic maps or $\mathrm{Beta}(p,q)$ maps $x_{n+1}=r{x}_n^{p-1}{(1-x_n)}^{q-1}$ have been investigated by Aleixo et al. [92,93,94], Brilhante et al. [19,20,21], Rocha and Aleixo [95], and Gomes et al. [22].

Further, the hyper-Gompertz map $r{G}_{p,1,1,\mathsf{\Psi }}:(0,1]\to (0,1]$, $r\in \left(0,{\left(\frac{p-1}{\mathsf{\Psi }-1}\mathrm{e}\right)}^{\mathsf{\Psi }-1}\right)$, $p,\mathsf{\Psi }>1$, $x_{n+1}=r{x}_n^{p-1}{\left(-ln{x}_n\right)}^{\mathsf{\Psi }-1}$ (Waliszewski and Konarski [12,13]; Pestana et al. [96]; Rocha and Aleixo [97]) is connected to the Turner et al. [17,18] hyper-Gompertz generalized ecological growth model $\frac{\mathrm{d}}{\mathrm{d}t}N\left(t\right)=rN\left(t\right){\left[-ln\frac{N\left(t\right)}{K}\right]}^{\gamma },r>0,\gamma \ne 1$, which we reparameterize as Equation (5),

$$\frac{\mathrm{d}}{\mathrm{d}t}N\left(t\right)=r{\left[N\left(t\right)\right]}^{p-1}{\left[-ln\left(\frac{N\left(t\right)}{K}\right)\right]}^{\mathsf{\Psi }-1}.$$

The discussion of the Blumberg [16], Turner et al. [17,18] and other population growth models and their associated maps is postponed to Section 6.

5. Stable and Geo-stable EV Distributions

5.1. Extremes of IID Sequences

Let ${\left\{X_n\right\}}_{n=1}^{\infty }$ be a sequence of IID RVs. For any $n\in \mathbb{N},$ $F_{X_{n:n}}\left(x\right)=F_X^n\left(x\right)$, with degenerate limit when $n\to \infty$. Non-degenerate limit EV laws are obtained in the context of Khinchine-type limits of classes, when there exist attraction coefficients $a_n>0,b_n\in \mathbb{R}$ such that

$$F_{\frac{X_{n:n}-b_n}{a_n}}\left(x\right)=F_X^n(a_{n}x+b_n)\underset{n\to \infty }{⟶}G\left(x\right)\mathrm{non}\mathrm{degenerate},$$

and we say that $F_X$ is in the domain of attraction of the MS CDF G.

The possible non-degenerate limit distributions G satisfy the max-stability equation

$$G(a_{\theta }x+b_{\theta })=G^{\frac{1}{\theta }}\left(x\right).$$

(11)

When $b_{\theta }$, the shift parameter, is zero, then $G\left(a_{\theta }x\right)=G^{\frac{1}{\theta }}\left(x\right)$ (strict max-stability equation) and we say that the CDF G is a strict MS distribution.

Fréchet [98] obtained a first solution of the aforementioned strict max-stability equation, given by $G_{\xi }\left(x\right)=exp\left(-x^{-1/\xi }\right){\mathbb{I}}_{[0,\infty )}\left(x\right),\xi >0$, rightly called Fréchet CDF. Fisher and Tippet [99] have shown that CDFs G of the type

$$G\left(x\right)\equiv G_{\xi }\left(x\right)=exp\left[-{(1+\xi x)}^{-1/\xi }\right]{\mathbb{I}}_{\{x:1+\xi x>0\}}\left(x\right),\xi \in \mathbb{R},$$

(12)

which, when $\xi \to 0$, converges to the Gumbel distribution, $G_0\left(x\right)=exp\left(-{\mathrm{e}}^{-x}\right)$ for real x, satisfy the functional max-stability Equation (11), and Gnedenko [52] has shown that the general extreme value (GEV) distributions in (12) are the unique MS laws.

Expression (12) is known as von Mises–Jenkinson GEV family of distributions. It is sometimes more convenient to split it into three separate branches:

Fréchet-$\alpha$ distributions: for $\xi >0$, and using $\alpha =\frac{1}{\xi }$,

$${\mathsf{\Phi }}_{\alpha }\left(x\right)=exp\left(-x^{-\alpha }\right){\mathbb{I}}_{[0,\infty )}\left(x\right);$$

(13)

Gumbel distribution: for $\xi =0$,

$$\mathsf{\Lambda }\left(x\right)=exp\left(-{\mathrm{e}}^{-x}\right){\mathbb{I}}_{\mathbb{R}}\left(x\right);$$

(14)

max-Weibull-$\alpha$ distributions: for $\xi <0$ and using $\alpha =-\frac{1}{\xi }$,

$${\mathsf{\Psi }}_{\alpha }\left(x\right)=\left\{\begin{array}{cc}exp\left[-{(-x)}^{\alpha }\right]& x<0\\ 1& x\ge 0\end{array}\right..$$

(15)

As $min\{X_1,\dots ,X_n\}=-max\{-X_1,\dots ,-X_n\}$, we obviously have $F_{min\{X_1,\dots ,X_n\}}\left(x\right)=1-F_{max\{X_1,\dots ,X_n\}}(-x)$, and so, in what concerns min-stability, we get:

• min-Fréchet-$\alpha$ distributions: ${\mathsf{\Phi }}_{\alpha }^{*}\left(x\right)=\left\{\begin{array}{cc}1-exp\left[-{(-x)}^{-\alpha }\right]& x<0\\ 1& x\ge 0\end{array}\right.$;
• min-Gumbel distribution: ${\mathsf{\Lambda }}^{*}\left(x\right)=1-exp\left(-{\mathrm{e}}^x\right){\mathbb{I}}_{\mathbb{R}}\left(x\right)$;
• Weibull-$\alpha$ distributions: ${\mathsf{\Psi }}_{\alpha }^{*}\left(x\right)=\left[1-exp\left(-x^{\alpha }\right)\right]{\mathbb{I}}_{[0,\infty )}\left(x\right)$.

Classical EVT deals essentially with the identification of the extreme stable laws and the characterization of the domains of attraction of the Fréchet and Weibull-types achieved by Gnedenko [52], and of the Gumbel-type due to de Haan [53]. For more on EV distributions, see Johnson et al. [35], pp. 1–112. See also the recent overviews by Davison and Huser [100] and by Gomes and Guillou [101].

5.2. Extremes of Geometrically Thinned Sequences

If $X_1,X_2,\dots ,X_n,\dots$ are independent replicas of X, with CDF $F_X$, and $N⌢\mathrm{Geometric}\left(\theta \right)$ independent of the $X_k$’s, the CDF of $Y=max\{X_1,\dots ,X_N\}$ is

$$\sum _{k=1}^{\infty }{F}_X^k\left(x\right)\theta {(1-\theta )}^{k-1}=\frac{\theta F_X\left(x\right)}{1-(1-\theta )F_X\left(x\right)}.$$

(16)

From (16), we say that Y is a GMS RV (or that $F_Y$ is a GMS CDF) if for all $\theta \in (0,1)$ there exist $a_{\theta }>0$ and $b_{\theta }\in \mathbb{R}$ such that

$$F_Y(a_{\theta }x+b_{\theta })=\frac{\theta F_Y\left(x\right)}{1-(1-\theta )F_Y\left(x\right)}.$$

(17)

In other words, an RV X is GMS iff it shares the same Khinchine’s type with the geometric randomly stopped maxima of independent replicas of itself.

The core theory of max-geo-stability has been developed by Rachev and Resnick [38], and results on thinned sequences with count random subordinator are due to Rényi [102], to Kovalenko [103] and, in all generality, to Kozubowski [104], see also Gnedenko and Korolev [105] for a thorough discussion.

Let us define $G\left(x\right)={\mathrm{e}}^{1-1/F_Y\left(x\right)},x>{\alpha }_{F_Y}$, where ${\alpha }_F=inf\{x:F\left(x\right)>0\}$ is the left-endpoint of F, and rewrite (17) as

$$\frac{1}{1-ln\left(G(a_{\theta }x+b_{\theta })\right)}=\frac{\frac{\theta }{1-ln\left(G\left(x\right)\right)}}{1-(1-\theta )\frac{1}{1-ln\left(G\left(x\right)\right)}}⟺\frac{1}{1-ln\left(G(a_{\theta }x+b_{\theta })\right)}=\frac{1}{1-\frac{ln\left(G\left(x\right)\right)}{\theta }}.$$

(18)

Then, (17) is equivalent to the max-stability Equation (11), $G(a_{\theta }x+b_{\theta })=G^{\frac{1}{\theta }}\left(x\right)$, i.e., G is an MS distribution.

Hence, the GMS CDFs, ${}^g{G}_{\xi }\left(x\right)=\frac{1}{1-ln{G}_{\xi }\left(x\right)}$, $x>{\alpha }_F$, are given by

$${}^g{G}_{\xi }\left(x\right)\equiv F_{\xi }\left(x\right)=\frac{1}{1-ln{G}_{\xi }\left(x\right)}=\frac{1}{1+{(1+\xi x)}^{-1/\xi }},1+\xi x>0.$$

(19)

The GMS CDFs in (19), using $\alpha =\frac{1}{\left|\xi \right|}$, can thus be written as one of the following types:

Log-logistic distributions, whose natural logarithm follows the logistic distribution, from the classical MS Fréchet-$\alpha$ distribution,

$${}^g{\mathsf{\Phi }}_{\alpha }\left(x\right)=\frac{1}{1+x^{-\alpha }}{\mathrm{I}}_{[0,\infty )}\left(x\right),\alpha >0;$$

Logistic distribution, from the classical MS Gumbel distribution,

$${}^g\mathsf{\Lambda }\left(x\right)=\frac{1}{1+{\mathrm{e}}^{-x}}{\mathrm{I}}_{\mathbb{R}}\left(x\right);$$

Backward log-logistic distributions, from the MS max-Weibull-$\alpha$ distribution,

$${}^g{\mathsf{\Psi }}_{\alpha }\left(x\right)=\left\{\begin{array}{cc}\frac{1}{1+{(-x)}^{\alpha }}& x<0\\ 1& x\ge 0\end{array}\right.,\alpha >0.$$

The characterization of the domains of attraction of GMS laws are similar to the characterization of the domains of attraction of the classical maxima EV laws, a consequence of tail equivalence results from Resnick [106,107], see also Cline [108].

Obviously, the maximum of a geometrically thinned sequence is necessarily stochastically smaller than the maximum of the full sequence, for the same parent population. Therefore, GMS laws are stochastically smaller than the corresponding classical EV laws, as can be seen in Figure 2.

Noticing again that $min\{X_1,\dots ,X_n\}=-max\{-X_1,\dots ,-X_n\}$, similar results exist for min-geo-stability. The min-geo-stable laws are ${}^g{G}_{\xi }^{*}\left(x\right)=1-{}^g{G}_{\xi }(-x)$, with ${}^g{G}_{\xi }$ given in (19). We thus get ${}^g{\mathsf{\Phi }}_{\alpha }^{*}={}^g{\mathsf{\Psi }}_{\alpha }$, ${}^g{\mathsf{\Lambda }}^{*}={}^g\mathsf{\Lambda }$ and ${}^g{\mathsf{\Psi }}_{\alpha }^{*}={}^g{\mathsf{\Phi }}_{\alpha }$; and obviously, min-geo-stable laws are stochastically greater than the corresponding classical minimum EV laws, in what regards stochastic ordering.

5.3. Verhulst Growth and Thinned Maxima, Gompertz Growth, and Maxima of IID Sequences

In particular, in what concerns population growth models (2) and (3), we can now state that the logistic and the Gumbel solutions are both EV models, respectively, in the thinned IID framework and in the full IID framework. In $f_{2,1,1,2}$, the expansion factor x dominates the retroaction $-lnx$ and the PDF is skewed to the right, while the symmetry of $f_{2,2,1,1}$ implies perfect equilibrium of the expansion factor x and the retroaction factor $(1-x)$. Therefore, $\mathbb{P}[X_{2,2,1,1}>x]>\mathbb{P}[X_{2,1,1,2}>x]$ for all $x\in (0,1)$, and this implies that the Gompertz model stochastically dominates the Verhulst model—as stated before, ${}^g\mathsf{\Lambda }⪯\mathsf{\Lambda }$. Looking at Figure 3, it is obvious that the retroaction factor $-lnx$ in $f_{2,1,1,2}$ exerts a weaker control than the retroaction factor $1-x$ in $f_{2,2,1,1,}$ on the growth of the population. This explains the success of the Gompertz (Gumbel) model in cancerous tumor growth (Laird [5]; Laird et al. [6]; Bajzer et al. [109]; Waliszewski et al. [9,10]; Waliszewski and Konarski [12]).

Further, Mejzler [110,111,112] characterized the class $\mathcal{M}$ of limit laws of maxima of non-identically distributed independent RVs, similar to Khinchine’s self-decomposable $\mathcal{L}$-class in the asymptotic additive theory, in terms of log-concavity. Graça Martins and Pestana [113] investigated EVT self-decomposable distributions of Mejzler’s class using an idea of Urbanik [114] on refinements of Khinchine’s $\mathcal{L}$ class of additive self-decomposable characteristic functions. In this framework, ${\mathcal{M}}_r$ is the class of CDFs $F\left(x\right)=exp(-K(x\left)\right)$ with K monotone of order r—thus, $M_{\infty }$ is the class of CDFs $F\left(x\right)=exp(-K(x\left)\right)$ with K completely monotone, and so a Laplace transform of a non-decreasing function. Observe that ${\mathrm{e}}^{-x}$, $x^{-\alpha }$, and ${(-x)}^{\alpha }$ are completely monotone, and therefore, ${\mathcal{M}}_{\infty }$ is a natural but non-trivial extension of EV models.

6. Population Growth Models and New Routes to Chaos

The logistic (2) and the Gompertz (3) population growth models, and their Blumberg hyper-logistic and hyper-Gompertz–Turner extensions are special cases of the Tsoularis and Wallace [39] generalized logistic equation (Section 6.1). The basics on the association of the logistic map to the Verhulst DE in Section 6.2 suggest the association of generalized logistic maps to investigate the dynamics of generalized logistic population growth; the Schwarzian derivative is used to investigate the dynamics of the hyper-logistic and of the hyper-Gompertz maps. In Section 6.3, we investigate population models’ connection to EV laws. In Section 6.4, we consider incursive maps and Dubois’ [28,44] denial that chaotic complicated dynamics arises from the simple discrete version of the Verhulst equation.

6.1. Generalized Logistic Population Models

Tsoularis and Wallace [39], p. 22, defined a generalized logistic population model $N\left(t\right)$

$$\frac{\mathrm{d}}{\mathrm{d}t}N\left(t\right)=r{\left[N\left(t\right)\right]}^{\alpha }{\left[1-{\left(\frac{N\left(t\right)}{K}\right)}^{\beta }\right]}^{\gamma },\alpha ,\beta ,\gamma >0.$$

(20)

Special cases relevant in our framework are:

1.

$\gamma =1$, the model $\frac{\mathrm{d}}{\mathrm{d}t}N\left(t\right)=r{\left[N\left(t\right)\right]}^{\alpha }\left[1-{\left(\frac{N\left(t\right)}{K}\right)}^{\beta }\right]$, investigated by Richards [30], with solution $N\left(t\right)=K{\left\{1-\left[1-{\left(\frac{K}{N_0}\right)}^{\beta }\right]{\mathrm{e}}^{-\beta rt}\right\}}^{-1/\beta }$, where $N_0=N\left(0\right)$, i.e., the initial population size. The special case $\alpha =\beta =\gamma =1$ is the Verhulst logistic model, and the special case $\alpha =1$ and $\gamma =0$ is the Malthus exponential growth model.

2.

$\beta =1$, $\frac{\mathrm{d}}{\mathrm{d}t}N\left(t\right)=r{\left[N\left(t\right)\right]}^{\alpha }{\left[1-\frac{N\left(t\right)}{K}\right]}^{\gamma }$, Blumberg [16] hyper-logistic DE, in (4), whose dynamics has been investigated by Aleixo et al. [92,93,94] and Rocha and Aleixo [95].

3.

$\alpha =1+\beta (1-\gamma )$, the model $\frac{\mathrm{d}}{\mathrm{d}t}N\left(t\right)=r{\left[N\left(t\right)\right]}^{1+\beta (1-\gamma )}{\left[1-{\left(\frac{N\left(t\right)}{K}\right)}^{\beta }\right]}^{\gamma }$, $\alpha ,\beta ,\gamma >0$, extensively studied by Turner et al. [17,18] under the name generic growth function, with solution

$$N\left(t\right)=K{\left\{1+{\left[\beta (\gamma -1)K^{\beta (1-\gamma )}rt+{\left[{\left(\frac{K}{N_0}\right)}^{\beta }-1\right]}^{1-\gamma }\right]}^{1/(1-\gamma )}\right\}}^{-1/\beta }.$$

The case $\frac{\mathrm{d}}{\mathrm{d}t}N\left(t\right)=r{\left[N\left(t\right)\right]}^{\beta -1}{\left[1-{\left(\frac{N\left(t\right)}{K}\right)}^{\beta }\right]}^{\frac{2}{\beta }}$ (i.e., $1+\beta (1-\gamma )=\beta -1,\beta >1$) is tied to the PDF of the ${\mathcal{K}}_{\beta ,1+\frac{2}{\beta }}$ Kumaraswamy RV.

4.

If $\beta \downarrow 0$, the model $\frac{\mathrm{d}}{\mathrm{d}t}N\left(t\right)=r{\left[N\left(t\right)\right]}^{\alpha }{\left[-ln\left(\frac{N\left(t\right)}{K}\right)\right]}^{\gamma }$ (for $\alpha =\gamma =1$, this is the Gompertz model). If $\alpha =1$, this hyper-Gompertz DE has the solution

$$N\left(t\right)=Kexp\left\{-{\left[(\gamma -1)rt+{\left[-ln\left(\frac{N_0}{K}\right)\right]}^{1-\gamma }\right]}^{\frac{1}{1-\gamma }}\right\}.$$

For $\alpha =\gamma =1$, this is the Gompertz model.

6.2. The Logistic Paradigm and Extensions

May [27] pathbreaking investigation of the surprising very complicated dynamics of the rather simple logistic map connected to the Verhulst population growth model originated a renewal of interest on population growth equilibrium, see the outstanding overview of growth models in Tsoularis and Wallace [39] and references therein.

6.2.1. Logistic Growth

Consider a population of size $N\left(t\right)$, and let $R\left(t\right)$ be the amount of resources at time t. It is sensible to consider that the amount of resources available decreases with the population growth, the simplest model being that $R\left(t\right)$ decreases at constant rate with respect to the population size derivative, $\frac{\mathrm{d}}{\mathrm{d}t}R\left(t\right)=-\xi \frac{\mathrm{d}}{\mathrm{d}t}N\left(t\right)$, where $\xi$ is the amount of resources consumed to produce a new population unit. Solving $R\left(t\right)$ in terms of $N\left(t\right)$ yields $R\left(t\right)=\xi [K-N(t\left)\right]=R\left(0\right)-\xi N\left(t\right)$, where $K=\frac{R\left(0\right)}{\xi }$ is the carrying capacity, the overpopulation size that consumes all the resources.

On the other hand, it is natural to assume that the population growth rate $\frac{\frac{\mathrm{d}}{\mathrm{d}t}N\left(t\right)}{N\left(t\right)}$ is proportional to the availability of nutrients, i.e., $\frac{\frac{\mathrm{d}}{\mathrm{d}t}N\left(t\right)}{N\left(t\right)}=KR\left(t\right)$, and hence, $\frac{\mathrm{d}}{\mathrm{d}t}N\left(t\right)=KR\left(0\right)N\left(t\right)\left(1-\frac{N\left(t\right)}{R\left(0\right)/\xi }\right)$, which we shall rewrite in the traditional Verhulst form in (2):

$${\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}N\left(t\right)=rN\left(t\right)\left(1-\frac{N\left(t\right)}{K}\right)⟹N\left(t\right)=\frac{KN\left(0\right)}{N\left(0\right)+(K-N\left(0\right)){\mathrm{e}}^{-rt}}.}$$

It is sometimes more convenient to rewrite (2) in terms of the rate $\alpha \left(t\right)=\frac{N\left(t\right)}{K}\in [0,1]$,

$$\frac{\mathrm{d}}{\mathrm{d}t}\alpha \left(t\right)=r\alpha \left(t\right)\left(1-\alpha \left(t\right)\right),$$

(21)

and for lattice observations of the populations size, using $\alpha \left(n\right)=\frac{1+\rho }{\rho }{x}_n,r=1+\rho$, to approximate (21) with the difference equation

$$x_{n+1}=r{x}_n(1-x_n),n=0,1,\dots ,$$

(22)

known as the logistic map. As pointed out before, the logistic parabola $G_{2,2,1,1}\left(x\right)=x(1-x){\mathbb{I}}_{(0,1)}\left(x\right)$ is proportional to the PDF $f_{U_{2:3}}\left(x\right)=6x(1-x){\mathbb{I}}_{(0,1)}\left(x\right)$ of the middle OS $U_{2:3}$ of $n=3$ IID standard uniform RVs.

It is obvious that the steady state $x_{n+1}=x_n$ is attained either at the root $x_n=0$ or at the root $x_n=1-\frac{1}{r}$. It seems far-fetched to use the fixed point method and the iterative expression (22) to solve $x=f\left(x\right)=rx(1-x)$, namely because whenever $r\notin (1,3)$, we have that $\left|f^{\prime }\left(1-\frac{1}{r}\right)\right|\ge 1$, and hence, numerical instabilities are observed outside the range $(1,3)$. As the condition $x_n\in [0,1],n=1,2,\dots$, is fulfilled for $r\in (0,4]$, for $r\in [3,4]$, we observe that as r increases, there are first periodic oscillations in the sequence of iterates (Feigenbaum bifurcations), and then, for $r>1+\sqrt{5}$, chaos. For details, consult Schroeder [46], ch. XII, or Peitgen et al. [91], ch. 10–11.

6.2.2. Gompertz and Hyper-Logistic Growth

The assumption on the decrease of resources can be modified, for instance, assuming a logarithmic or polynomial decrease rate of resources.

If instead of $\frac{\mathrm{d}}{\mathrm{d}t}R\left(t\right)=-\xi \frac{\mathrm{d}}{\mathrm{d}t}N\left(t\right)$, we assume now that $R\left(t\right)$ decreases at constant rate with the intrinsic growth rate of the population, $\frac{\mathrm{d}}{\mathrm{d}t}R\left(t\right)=-\xi \frac{\frac{\mathrm{d}}{\mathrm{d}t}N\left(t\right)}{N\left(t\right)}$, whose solution is $R\left(t\right)=-\xi ln\left(N\left(t\right)\right)+c=-\xi ln\left[\frac{N\left(t\right)}{M}\right]$ since at overpopulation M, all resources are consumed. Hence, in terms of the rate $\alpha \left(t\right)=\frac{N\left(t\right)}{K}\in [0,1]$,

$${\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}\alpha \left(t\right)=\xi \alpha \left(t\right)\left[-ln\left(\alpha \right(t\left)\right)\right]⟹\alpha \left(t\right)={\mathrm{e}}^{{-\mathrm{e}}^{-\xi t}},}$$

the Gompertz (or Gumbel) CDF, see Section 4, and below a discussion of the Gompertz map.

With polynomial decreasing resources

$$\frac{\mathrm{d}}{\mathrm{d}t}R\left(t\right)=-\xi {[K-N\left(t\right)]}^{\gamma -1}\frac{\mathrm{d}}{\mathrm{d}t}N\left(t\right)⟹R\left(t\right)=\frac{\xi }{\gamma }{[K-N\left(t\right)]}^{\gamma }$$

(obviously, the case $K=1$ is what we assumed for the logistic model) and in terms of the rate $\alpha \left(t\right)=\frac{N\left(t\right)}{K}\in [0,1]$, we obtain the hyper-logistic population dynamics model

$$\frac{\mathrm{d}}{\mathrm{d}t}\alpha \left(t\right)=K\frac{\xi K^{\gamma }}{\gamma }\alpha \left(t\right){\left[1-\alpha \left(t\right)\right]}^{\gamma },$$

a Blumberg equation, see case 2 associated with (20).

6.2.3. The Logistic Map with Random Reproducing Rate

Rewriting (22) as

$$x_{n+1}=4r^{*}{x}_n(1-x_n),n=0,1,\dots ,r^{*}\in (0,1],$$

the parameter $r^{*}$ reflecting the population reproduction rate; in many situations, it makes sense to consider that, instead of a parameter, the reproduction rate is a random variable R with support $(0,1)$. For instance, contrarily to intuition, reproduction rate can increase with population density (Allee effect, see Stephens et al. [115], Berec et al. [116], Kramer et al. [117] and Aleixo et al. [92,94]).

The family of PDFs of Beta RVs, $X⌢\mathrm{Beta}(p,q)$, in (8), provides a huge diversity of shapes, and has the advantage of staying within the family of Blumberg hyper-logistic population models. We shall therefore consider the system

$$x_{n+1}=4R{x}_n(1-x_n),n=0,1,\dots ,R⌢Beta(p,q).$$

A judicious choice of parameters so that the mode $\frac{p+1}{p+q-2}$ matches the empirical mode can be a useful model. Alternatively, previous empiric information on the central location and on the dispersion of R can serve as a guideline to choose the range of the parameters p and q.

6.2.4. Schwarzian Derivative of Hyper-Logistic and Hyper-Gompertz Maps

The orbit of any point $x\in \mathbb{R}$, i.e., its behavior under iteration of f, is given by the sequence $x,f\left(x\right),f\left(f\right(x\left)\right),f\left(f\right(f\left(x\right)\left)\right),\dots$. We denote the points in the orbit $x,f\left(x\right),f^{\circ 2}\left(x\right),$$f^{\circ 3}\left(x\right),\dots$. If $f\left(p\right)=p$, implying that $f^{\circ n}\left(p\right)=p$, we say that p is a fixed point for f, and q is a periodic point of period n for f if $f^{\circ n}\left(q\right)=q$. This is the prime period of q if it is the least positive integer n for which $f^{\circ n}\left(q\right)=q$.

The Schwarzian derivative of a smooth interval map $f:(0,1)\to [0,1]$, $f\in {\mathcal{C}}^3$ is defined as

$$S\left(f\left(x\right)\right)=\frac{f^{\prime \prime \prime }\left(x\right)}{f^{\prime }\left(x\right)}-\frac{3}{2}{\left(\frac{f^{\prime \prime }\left(x\right)}{f^{\prime }\left(x\right)}\right)}^2.$$

Singer [40] established its relevance in one-dimensional dynamics showing that the negative Schwarzian derivative of the map $f:I\to I$ implies that for every stable periodic point x of period n, there exists an $i<n$ and a critical point c or endpoint of I such that $y\in [c,f^{\circ i}\left(x\right)]$ implies $f^{\circ kn}\left(y\right)\to f^{\circ i}\left(x\right)$ as $k\to \infty$.

Further, consider the Sharkovskii total ordering (but not well-ordering) of the positive integers

$$\begin{array}{ccccccccc}3& \prec & 5& \prec & 7& \prec & 9& \prec \dots \prec & (2n+1)·2^0\prec \dots \\ 3·2& \prec & 5·2& \prec & 7·2& \prec & 9·2& \prec \dots \prec & (2n+1)·2^1\prec \dots \\ 3·2^2& \prec & 5·2^2& \prec & 7·2^2& \prec & 9·2^2& \prec \dots \prec & (2n+1)·2^2\prec \dots \\ 3·2^3& \prec & 5·2^3& \prec & 7·2^3& \prec & 9·2^3& \prec \dots \prec & (2n+1)·2^3\prec \dots \\ & ⋮\\ \dots & \prec & 2^n& \prec \dots \prec 2^4& \prec & 2^3& \prec & 2^2\prec 2& \prec 1\end{array}$$

used by Sharkovskii [42] to prove that the existence of period $p=3$ implies the existence of any other period length. However, most of those orbits are unstable, as Schroeder [46], pp. 285–286, poetically states, “they are the ‘ghosts’ of orbits that were stable for smaller r values”. Further, Li and Yorke [43] proved that period $p=3$ implies chaos: the existence of a period-3 cycle implies the existence of an uncountable infinitude of chaotic points that never map to any cycle. For further information on wandering intervals and sensitivity on initial conditions, consult Harrison [118], and for general information on Schwarzian derivatives, see Cooper [119] and Devaney [120] (ch. 11).

Below, we establish that the hyper-logistic and hyper-Gompertz maps verify Singer [40] and Guckenheimer [41] assumptions, implying thus the nonexistence of wandering intervals. From Sharkovskii [42] and Li and Yorke [43] results, they share with the logistic map instabilities, bifurcations, and chaos when, for the fixed point $x^{*}$, r increases beyond the values for which $\left|f^{\prime }\left(x^{*}\right)\right|<1$. For thorough information on the dynamics of one-dimensional maps, consult Sharkovsky et al. [121].

6.2.5. Hyper-Logistic Maps

In what concerns the family of unimodal maps

$$r{G}_{p,q,1,1}\left(x\right)=r{x}^{p-1}{\left(1-x\right)}^{q-1}:(0,1)\to (0,1),$$

with $r\in \left(0,{\left(\frac{p-1}{p+q-2}\right)}^{1-p}{\left(\frac{q-1}{p+q-2}\right)}^{1-q}\right)$ and $p,q>1$,

$$G{\prime }_{p,q,1,1}\left(x\right)=0\mathrm{for}x=c=\frac{p-1}{p+q-2},G_{p,q,1,1}^{\prime \prime }\left(c\right)<0$$

and

$${\left(\frac{p-1}{p+q-2}\right)}^{1-p}{\left(\frac{q-1}{p+q-2}\right)}^{1-q}{G}_{p,q,1,1}\left(c\right)=1.$$

The Schwarzian derivatives are

$$S\left(G_{p,q,1,1}\left(x\right)\right)=\frac{A}{2x^2{(x-1)}^2{[1+p(x-1)+(q-2)x]}^2},$$

where

$$\begin{array}{c}A=-p^4{(x-1)}^4-4p^3{(x-1)}^3[1+(q-2)x]-\hfill \\ -(q^2-3q+2)x^2[6+4(q-3)x+(q^2-5q+6)x^2]-\\ -p^2{(x-1)}^2[5+(12q-22)x+(6q^2-24q+23)x^2]-\\ \hfill -2p(x-1)[1+(4q-7)x+(6q^2-23q+20)x^2+(2q^3-12q^2+23q-14)x^3].\end{array}$$

See Figure 4 and Figure 5, and a sample of orbits and of bifurcations in Figure 6 and Figure 7.

6.2.6. Hyper−Gompertz Maps

In what concerns the family of unimodal maps

$$r{G}_{p,1,1,\mathsf{\Psi }}\left(x\right)=r{x}^{p-1}{\left(-lnx\right)}^{\mathsf{\Psi }-1}:(0,1)\to (0,1),$$

with $r\in \left(0,{\left(\frac{p-1}{\mathsf{\Psi }-1}\mathrm{e}\right)}^{\mathsf{\Psi }-1}\right)$ and $p,q>1$,

$$G_{p,1,1,\mathsf{\Psi }}^{\prime }\left(x\right)=0\mathrm{for}x=c={\mathrm{e}}^{-\frac{\mathsf{\Psi }-1}{p-1}},G_{p,1,1,\mathsf{\Psi }}^{\prime \prime }\left(c\right)<0$$

and

$${\left(\frac{p-1}{\mathsf{\Psi }-1}\mathrm{e}\right)}^{\mathsf{\Psi }-1}{G}_{p,1,1,\mathsf{\Psi }}\left(c\right)=1.$$

The Schwarzian derivatives are

$$S\left(G_{p,1,1,\mathsf{\Psi }}\left(x\right)\right)=\frac{B}{2x^2{(lnx)}^2{\left[\mathsf{\Psi }-1+(p-1)lnx\right]}^2},$$

where

$$\begin{array}{c}B=-[\mathsf{\Psi }{(\mathsf{\Psi }-1)}^2(\mathsf{\Psi }-2)+4(p-1)\mathsf{\Psi }({\mathsf{\Psi }}^2-3\mathsf{\Psi }+2)lnx+\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}+(6p^2-12p+5){(\mathsf{\Psi }-1)}^2{(lnx)}^2+(4p^3-12p^2+10p-2)(\mathsf{\Psi }-1){(lnx)}^3+\\ \hfill +p{(p-1)}^2(p-2){(lnx)}^4].\end{array}$$

See Figure 8 and Figure 9, and some orbits and bifurcations in Figure 10 and Figure 11.

6.3. Population Growth Models and EV Distributions

Rewriting $N\left(t\right)=\frac{K}{1+{\mathrm{e}}^{-rt-C}},\mathrm{where}{\mathrm{e}}^{-C}=\frac{K-N\left(0\right)}{N\left(0\right)}$, the solution of the basic Verhulst population growth DE in (2) is proportional to the GMS logistic CDF with location $-C/r$ and scale $1/r$.

The solution $N\left(t\right)=Kexp(-{\mathrm{e}}^{-rt+C})$ of the Gompertz DE in (3) is proportional to the Gumbel CDF in the classical IID framework, with location and scale parameters $-C/r$ and $1/r$, respectively.

Further, the solution $N\left(t\right)=\frac{K}{1+t^{-rr}}{\mathbb{I}}_{[0,\infty )}\left(x\right)$ of $\frac{\mathrm{d}}{\mathrm{d}\mathrm{t}}N\left(t\right)=\frac{r}{t}N\left(t\right)\left[1-\frac{N\left(t\right)}{K}\right]$ is proportional to the log-logistic CDF, also a stable distribution of maxima of geometrically thinned IID RVs.

On the other hand, $\frac{\mathrm{d}}{\mathrm{d}t}{N}_{\nu }\left(t\right)=r{x}^{\nu }{N}_{\nu }\left(t\right)\left[-ln\frac{N_{\nu }\left(t\right)}{K}\right],r>0,$ has solution $N_{\nu }\left(t\right)=exp\left({\mathrm{e}}^{-\frac{r{x}^{\nu +1}}{\nu +1}+C}\right)$, C constant. Hence, letting $\nu \to -1$ with $C=\frac{r}{\nu +1}$ and using L’Hôpital’s rule, we get ${\displaystyle \underset{\nu \to -1}{lim}-r\frac{x^{\nu +1}-1}{\nu +1}=-rlnx}$ and we obtain the solution $N_{-1}\left(t\right)=exp(-x^{-1}){\mathbb{I}}_{[0,\infty )}\left(x\right)$, the Fréchet-1 distribution of maxima in the IID classical Fréchet [98] and Fisher and Tippett [99] setting.

Three equations modeling three patterns of extremes: sustainable logistic growth with the Verhulst model (2), moderate EVT Gumbel growth with the Gompertz model (3), and heavy tailed EVT Fréchet-1 growth with the above ${\displaystyle \underset{\nu \to -1}{lim}\frac{\mathrm{d}}{\mathrm{d}t}{N}_{\nu }\left(t\right)=r{x}^{\nu }{N}_{\nu }\left(t\right)\left[-ln\frac{N_{\nu }\left(t\right)}{K}\right],}$$r>0,$ model. The EVT Weibull model has also been used to fit population growth, namely in forests (Yang et al. [3]; Payandeh and Wang [4]).

We can speculate that evolution and adaptation, together with the need to preserve ecological equilibrium frameworks, are at the source of sustainable logistic growth in most natural populations, while viral spread or cancer tumor cell growth tends to overpass the availability of resources.

On the other hand, in adverse situations such as extreme climate, spread of epidemics or catastrophic circumstances such as the periodic “el Niño”, populations can decrease to minima compatible with collective species survival. Observe that in subfigures (c) of Figure 6 and Figure 10, there is a fixed point to the left of the mode of that unimodal map. Observe also that the solution of

$$\frac{\mathrm{d}}{\mathrm{d}\mathrm{t}}N\left(t\right)=r\left[1-\frac{N\left(t\right)}{K}\right]\left[-ln\left(1-\frac{N\left(t\right)}{K}\right)\right]$$

is proportional to a min-Gumbel CDF. Figure 12 exhibits the pattern of bifurcations and ultimate chaos of the non-linear maps $x_{n+1}=r{x}_n(-ln{x}_n)$ and $x_{n+1}=r(1-x_n)[-ln(1-x_n)]$.

In what follows, we examine geo-extreme solutions of the Blumberg [16] hyper-logistic model, and EVT solutions of the Turner et al. [17,18] hyper-Gompertz model.

It is curious to observe that the special case $p=q=2$ of the Blumberg model corresponds to the Verhulst Equation (2) with logistic solution, and the special case $p=\mathsf{\Psi }=2$ in the Turner equation is the Gompertz Equation (3) with the solution proportional to the Gumbel CDF —in other words, the limiting case $\xi \to 0$ in the general GMS distribution (19) and in the GEV (12), respectively.

Below, we indicate that in order to obtain the other GMS solutions via the Blumberg equation, we must consider fractional exponents $p-1=1\pm \frac{1}{\alpha }$ and $q-1=1\mp \frac{1}{\alpha }$, and that the other EVT-connected solutions via the Turner equation are obtained with $p=1$ and $2\ne \mathsf{\Psi }\in \mathbb{R}$.

6.3.1. Blumberg Equation and Geo-Extreme Models

The solutions of the Blumberg Equation (4), $\frac{\mathrm{d}}{\mathrm{d}\mathrm{t}}N\left(t\right)=r{\left[N\left(t\right)\right]}^{p-1}{\left[1-\frac{N\left(t\right)}{K}\right]}^{q-1}$, for any choice of parameters $p,q>0$ such that $p+q=4$ is a geo-stable EV distribution. More precisely, using fractional exponents $p-1=1-\frac{1}{\alpha }$ and $q-1=1+\frac{1}{\alpha }$,

$$\frac{\mathrm{d}N\left(t\right)}{\mathrm{d}t}=r{\left[N\left(t\right)\right]}^{1-\frac{1}{\alpha }}{\left[1-\frac{N\left(t\right)}{K}\right]}^{1+\frac{1}{\alpha }}:$$

• When $\alpha >0$, we get a solution that is proportional to the log-logistic CDF with shape parameter $\alpha$, location parameter $-\frac{C}{r{K}^{-1/\alpha }}$, and scale parameter $\frac{\alpha }{r{K}^{-1/\alpha }}$;
• When $\alpha <0$, we get a solution that is proportional to the backward log-logistic CDF with shape parameter $-\alpha$, location parameter $-\frac{C}{r{K}^{-1/\alpha }}$, and scale parameter $-\frac{\alpha }{r{K}^{-1/\alpha }}$;
• When $\frac{1}{\alpha }\to 0$, the solution of the Verhulst equation is proportional to the logistic CDF.

Figure 13 exhibits the pattern of bifurcations and ultimate chaos of the non-linear map $x_{n+1}=r{x}_n^{1-\frac{1}{2}}{(1-x_n)}^{1+\frac{1}{2}}$, and the animation in Supplementary materials shows the corresponding evolution of different initial conditions as a function of r.

For mathematical details, see Appendix B.

6.3.2. Turner Equation and EV Models

• For $p=2$ and $\mathsf{\Psi }=2+\frac{1}{\alpha }$, in (5), the solution of

  $$\frac{\mathrm{d}}{\mathrm{d}\mathrm{t}}N\left(t\right)=rN\left(t\right){\left[-ln\left(\frac{N\left(t\right)}{K}\right)\right]}^{1+\frac{1}{\alpha }}$$

  is proportional to the EV Fréchet CDF, in (13), or to the max-Weibull CDF, in (15), when $\alpha >0$ or $\alpha <0$, respectively;
• When $\frac{1}{\alpha }\to 0$—in other words, when we have the Gompertz Equation (3)—the solution is proportional to the Gumbel CDF, in (14).

On the other hand, the modified Turner equation

$$\frac{\mathrm{d}N\left(t\right)}{\mathrm{d}t}=r\left[1-\frac{N\left(t\right)}{K}\right]{\left[-ln\left(1-\frac{N\left(t\right)}{K}\right)\right]}^{1+\frac{1}{\alpha }},$$

i.e., $q=2$ and $\mathsf{\Phi }=2+\frac{1}{\alpha }$, has the solution

$$N\left(t\right)=K\left\{1-exp\left[-{\left(-\frac{t-KC/r}{K\alpha /r}\right)}^{-\alpha }\right]\right\}.$$

Therefore,

• When $\alpha >0$, the solution is proportional to the min-Fréchet CDF with shape parameter $\alpha$, location parameter $KC/r$, and scale parameter $K\alpha /r$;
• When $\alpha <0$, the solution is proportional to the Weibull CDF with shape parameter $-\alpha$, location parameter $KC/r$, and scale parameter $-K\alpha /r$.
  Figure 14 exhibits the pattern of bifurcations and ultimate chaos of the Turner map $x_{n+1}=r{x}_n{(-ln{x}_n)}^{1+\frac{1}{2}}$ (see also the animation in Supplementary materials showing the corresponding evolution of different initial conditions as a function of r) and of the modified Turner map $x_{n+1}=r(1-x_n){[-ln(1-x_n)]}^{1+\frac{1}{2}}$.

For mathematical details, see Appendix C.

6.4. Chaos, Indeed?

Dubois [28,44] claimed that the logistic map (10) is not the correct discrete equivalent to the Verhulst logistic DE in (2). His incursive equation

$$x_{n+1}=r{x}_n\left(1-x_{n+1}\right)$$

can be transformed into

$$x_{n+1}=r{x}_n\left(1-\frac{r{x}_n}{1+r{x}_n}\right)=\frac{r{x}_n}{1+r{x}_n},$$

which has a stable solution

$$x_n=\frac{r^n{x}_o}{1+\frac{(r^n-1)}{1-\frac{1}{r}}}$$

(23)

for any value of r, since $\left|\frac{r}{{(1+r{x}_n)}^2}\right|<1$.

Thus, (23) is the discrete equivalent to the continuous logistic Verhulst population model, and Dubois [28] concluded that “the chaos emerging from the so-called chaos map is due to instabilities of the Euler algorithm and not from fundamental biological properties of the logistic differential equation”.

An investigation of a similar incursive discrete approximation to the Gompertz DE, whose solution is the Gompertz curve, i.e., the Gumbel distribution, when properly normalized, is less satisfactory, since the incursive solution we get approximates a Fréchet-1 distribution instead of the Gumbel distribution.

In fact, consider that the retroaction acts at $n+1$, obtaining the difference equation $x_{n+1}=r{x}_n(-ln{x}_{n+1})$, with the same stationary solutions of the Gompertz map $x_{n+1}=r{x}_n(-ln{x}_{n+1})$—but without bifurcation and ultimate chaos. From $x_{n+1}=r{x}_n(-ln{x}_{n+1})$, we obtain a solution $f_r\left(x\right)=rxW\left(\frac{1}{rx}\right)$, where W is Lambert’s Logarithmic Product function, taking real values for $x>-0.5$.

Figure 15a shows that $cxW\left(\frac{1}{cx}\right){\mathbb{I}}_{[0,\infty )}\left(x\right)$ is a CDF, which is log-concave, and therefore, it is a Mejzler self-decomposable non-stable EV law.

Observe, however, that it is not a good approximation to the Gumbel distribution in $[0,\infty ]$. In the graph, $c=2$ and the Gumbel distribution is such that the lines meet at the 0.9 quantile, which means that the scale parameter is 0.52688. Although the solution of the Gompertz DE is proportional to the Gumbel distribution, this non-chaotic solution of its discrete counterpart incursive difference equation is patently a rather poor approximation, even for quite large values.

On the other hand, this Logarithmic Product law is in the domain of attraction of the Fréchet distribution with shape parameter 1, whatever the value $c>0$, since

$$\underset{x\to \infty }{\overset{}{lim}}\frac{1-ctxW\left(1/\left(ctx\right)\right)}{1-cxW\left(1/\left(cx\right)\right)}=t^{-1}.$$

In Figure 15b, the comparison of the CDF $2xW\left(\frac{1}{2x}\right){\mathbb{I}}_{[0,\infty )}\left(x\right)$ and the Fréchet distribution ${\mathrm{e}}^{-\frac{0.44995}{x}}{\mathbb{I}}_{(0,\infty )}\left(x\right)$ (with the scale parameter chosen so that the lines cross at the common 0.9 quantile) shows that this approximation is quite good.

7. Conclusions and Open Problems

Hamlet, who suffered from May–Dubois split personality, annoyed everyone with his “Chaos or no chaos, that’s the question!”, causing poor Ophelia to die of acute boredom cursing Hamlet and the heralds of Doomsday (Foerster et al. [23]). Although we recognize that Dubois [28] presents strong arguments—namely that his incursive solution agrees with the logistic solution of Verhulst’s DE in (2)—the fact that following his methodology with the Gompertz-map the incursive recurrence solution is Fréchet-1, while the solution of the DE in (3) is Gumbel, is indeed a weak point and does not contribute to substantiate Dubois’ claims.

On the other hand, in biology, there are so many instances of the tendency of immoderate growth, from virus and tumor cells to dominant species or human overpopulation, that instabilities and ultimate chaos seem to be the ubiquitous biological result of such selfish framework (in here, the term “selfish” is used in a loose sense, as in Dawkin’s [122], The Selfish Gene, voted as the most influential science book of all times). Therefore, although we consider Dubois [28] a challenging paper deserving high credit, May’s [27] pathbreaking absurd idea of using the fixed point method to solve a trivial second-order equation has a scientific descent that fully illustrates Einstein’s quote “If at first the idea is not absurd, then there is no hope for it”.

Nonetheless, an interesting open problem is to investigate further incursive recurrence solutions of other nonlinear maps, namely whether they are the right way to approximate the DEs with difference equations.

In addition, a substantial part of this review is to explore EV laws as possible patterns of population growth. It is worth noticing that the widely used logistic and Gompertz (Gumbel) models correspond to the Blumberg [16] $p=q=2$ equation and Turner et al. [17,18] $p=\mathsf{\Psi }=2$ equation, respectively.

Further, fractional powers $p=2\mp \frac{1}{\alpha }$ and $q=2\pm \frac{1}{\alpha }$ in the Blumberg equation lead to log-logistic and backward log-logistic extreme geo-stable models multiplied by the carrying capacity K. In what regards the Turner et al. [17,18] equation, $p=2$ and $\mathsf{\Psi }=2\pm \frac{1}{\alpha }$, eventually fractional, lead to the remaining Fréchet-$\alpha$ and min-Weibull-$\alpha$ EV distributions multiplied by the carrying capacity.

The delimitation of several regions (sudden extinction; stability with symbolic sequences $C{L}^{\infty }$; stability with symbolic sequences $C{R}^{\infty }$; period doubling; chaotic; Allee effect caused extinction; differed extinction) exhibited in Figure 1.3 of Aleixo et al. [94] in what concerns the hyper-logistic maps $x_{n+1}=r{x}^p(1-x)$, is still a challenging open problem for the hyper-Gompertz map, namely for the EVT and geo-thinned EVT maps.

In what regards the observation that $1-x$ is the first-order approximation of $-lnx$, implying that the logistic model may be considered a first-order approximation of the Gompertz model, we might as well say that ${\mathrm{e}}^{-x}=1-x+\cdots$. Therefore, differential equation models connected to gamma instead of beta or BetaBoop functions present a new field of research.

Recent research work on the logistic equation, namely Area and Nieto [61], using Prabhakar fractional calculus, deserve intense follow-up, namely applying their ideas to the hyper-logistic and to the hyper-Gompertz maps addressed in the present research work. For more information on the area and working methodology, consult El-Sayed et al. [123] and Giusti et al. [124].

Golmankhaneh and Cattani [125] used an analogue of the classical Euler method in fractal calculus useful for solving fractal DEs or for finding approximate analytical solutions, studying a fractal logistic equation using fractal derivative instead of the usual derivative in the Verhulst equation. Golmankhaneh and Fernandez [126] investigated stable distributions in fractal Cantor sets. It would be interesting to develop EVT and geo-EVT stability in fractal sets, namely middle-$\kappa$ Cantor sets, and to characterize domains of attraction. In what regards extensions of EVT, an explicit description of the ${\mathcal{M}}_{\infty }$ Mejzler class, namely identification of its extreme points and obtaining a Choquet representation, would be an interesting development.

So far, we could not find closed results connected to the $G_{1,1,\mathsf{\Phi },\mathsf{\Psi }}$ function, even for simple choices of $\mathsf{\Phi },\mathsf{\Psi }\ne 1$. On the other hand, it would be interesting to investigate further the generalized Pareto growth model investigated in Brilhante et al. [21], since the generalized Pareto distribution is closely tied to the GEV distribution.