1 Introduction
In 1925, Udny Yule published the paper [45] in which he described a possible model for macroevolution. Genera (species grouped by similar characteristics) appear in the system at random times according to a linear birth process. Each genus is initially composed by a single species. As soon as the genus appears, an independent linear birth process modelling the evolution of the species belonging to it, starts. The initial species thus generates offsprings (representing related species) with constant individual intensities. The model is now classical and it is usually named after him. The Yule model therefore admits the existance of two independent and superimposed mechanisms of evolution, one for genera and one for species, both at possibly different constant intensities.
One of the most interesting characteristics the model exhibits is its intrinsic “preferential attachment” mechanism. And this is exactly implied by the presence of the two distinct and superimposed random growths for the appearance of genera and of species.
The preferential attachment mechanism is a fundamental ingredient of many modern models of random graphs growth. The literature is vast and covers many fields. We recall here only some relevant papers and books [1, 2, 5, 6, 8, 12, 13, 16, 18, 23–25, 34, 37, 44] leaving the reader the possibility of widening the number of sources by looking at the references cited therein.
The Barabási–Albert model of random graph growth is by far the most famous example of a stochastic process based on preferential attachment [2, 7]. In [35] we discuss the relationships between the Barabási–Albert graph, the Yule model and a third model based on preferential attachment introduced by Herbert Simon in 1955 [26, 42, 43]. In [36] we have further analyzed and described the exact relation between the Barabási–Albert model and the Yule model. Briefly, the finite-dimensional distributions of the degree of a vertex in the Barabási–Albert model converges to the finite-dimensional distributions of the number of individuals in a Yule process with initial population size equal to the number m of attached edges in each time step. This further entails that the asymptotic degree distribution of a vertex chosen uniformly at random in the Barabási–Albert model coincides with the asymptotic distribution of the number of species belonging to a genus chosen uniformly at random from the Yule model with an initial number m of species. This result suggests that asymptotic models similar to the Yule model can be linked to different preferential attachment random graph processes in discrete time.
Therefore, a direct analysis of the Yule model, a model in continuous time and hence possessing a greater mathematical treatability, has the potential to uncover important aspects and characteristics of the discrete-time model to which it is related.
In [27] and [28] we have started a study of macroevolutionary models similar to the classical Yule model where the process governing the appearance of genera is left unchanged, while those describing the growth of species account for more realistic features. Specifically, in [27] we have generalized the latter allowing the possibility of extinction of species while in [28] we have studied the effect of a slowly-decaying memory by considering a fractional nonlinear birth process. Notice that, by suitably specializing the nonlinear rates, other peculiar behaviours such as saturation or logistic growth may be observed.
In this paper we proceed with the analysis by looking at modifications of the classical Yule model in which the appearance of genera follows a different dynamics. We will show that a change in the dynamics of genera will lead to radical changes in the model. This is a preliminary step before aiming at deriving models of random graphs with different features than those graphs connected to the Yule model. In the following we consider the class of mixed Poisson processes time-changed by means of a deterministic function. The rationale which justifies this choice will be clear in Section 3 where we state the main results. Section 2 contains the mathematical background necessary to develop the results presented later. In particular we will connect a member of the class of the suitably time-changed mixed Poisson processes with the time-fractional Poisson process, a non-Markov renewal process governed by a time-fractional difference-differential equations involving the Caputo–Džrbašjan derivative.
2 Preliminaries
We consider two different classes of point processes, namely the time-fractional Poisson processes (shortly tfPp) and an extension of the mixed Poisson processes, i.e. the mixed Poisson processes up to a (deterministic) time transformation (mPp-utt in the following). We will analyze briefly their properties and state a result linking the two classes. Besides, we will introduce the definitions and mathematical tools needed to understand Section 3.
2.1 Time-fractional Poisson processes (tfPp)
The tfPp has been introduced in the literature in [41, 40] (see also Laskin’s paper [29]). We show here the construction by means of random time-change with an inverse stable subordinator [31]. Alternatively, the tfPp can be defined as a specific renewal process (see e.g. [3, 4, 30], and see [31] for the proof of the equivalence between the two constructions).
Let us consider a homogeneous Poisson process ${({\widetilde{N}_{t}})}_{t\ge 0}$ of parameter $\lambda >0$ and an independent inverse stable subordinator, that is a one-dimensional time-continuous stochastic and non-Markov process defined as follows. Consider the subordinator ${({D_{t}})}_{t\ge 0}$ with Lévy measure $\nu (\text{d}x)=[\alpha /\varGamma (1-\alpha )]{x^{-1-\alpha }}$, $\alpha \in (0,1)$, and define its stochastic inverse ${({E_{t}})}_{t\ge 0}$ as the first time at which it exceeds a given threshold, i.e.
Now, consider the time-changed point process $\mathcal{N}={({\mathcal{N}_{t}})}_{t\ge 0}={({\widetilde{N}_{{E_{t}}}})}_{t\ge 0}$. The process $\mathcal{N}$ is called tfPp of parameters λ and α.
Many properties are known for the tfPp. Let us review some of them. The marginal probability distribution of the process generalizes the Poisson distribution (for $\alpha \to 1$) and can be written as [3, 4]
where
is the Prabhakar function [38] for complex parameters $\nu ,\beta ,\gamma $, with $\mathrm{\Re }(\nu )>0$. It is interesting to note that the mean value of the process grows in time less than linearly for each allowed value of α,
and that the variance can be written as
thus highlighting the overdispersion of the process. Furthermore, the probability generating function reads
where ${E_{\alpha }}(z)={E_{\alpha ,1}^{1}}(z)$ is the classical Mittag-Leffler function.
(2)
\[ \mathbb{P}({\mathcal{N}_{t}}=k)={\big(\lambda {t^{\alpha }}\big)^{k}}{E_{\alpha ,\alpha k+1}^{k+1}}\big(-\lambda {t^{\alpha }}\big),\hspace{1em}k\ge 0,\hspace{2.5pt}t\ge 0,\](3)
\[ {E_{\nu ,\beta }^{\gamma }}(z)={\sum \limits_{r=0}^{\infty }}\frac{{(\gamma )_{r}}{z^{r}}}{r!\varGamma (\nu r+\beta )},\hspace{2em}{(\gamma )_{r}}=\frac{\varGamma (\gamma +r)}{\varGamma (\gamma )},\hspace{1em}z\in \mathbb{C},\](4)
\[ \mathbb{E}({\mathcal{N}_{t}})=\lambda {t^{\alpha }}/\varGamma (\alpha +1),\hspace{1em}t\ge 0,\](5)
\[ \mathbb{V}\text{ar}({\mathcal{N}_{t}})=\mathbb{E}({\mathcal{N}_{t}})+\frac{{(\lambda {t^{\alpha }})^{2}}}{\alpha }\bigg(\frac{1}{\varGamma (2\alpha )}-\frac{1}{\alpha \varGamma {(\alpha )^{2}}}\bigg),\]Considering the renewal nature of $\mathcal{N}$, and calling ${U_{j}}$, $j\ge 1$, the random inter-arrival times between the $(j-1)$-th and the j-th event, it is possible to give an explicit expression for the common probability density function. Indeed,
where ${E_{\alpha ,\alpha }}(z)={E_{\alpha ,\alpha }^{1}}(z)$ is the generalized Mittag-Leffler function. By analyzing the density, its slowly decaying right tail (it is actually an ultimately monotone regularly varying function of order $-\alpha -1$) and the asymptote in zero, the clusterization of events in time appears evident.
(7)
\[ {f_{{U_{j}}}}(t)=\lambda {t^{\alpha -1}}{E_{\alpha ,\alpha }}\big(-\lambda {t^{\alpha }}\big){\mathbb{I}_{{\mathbb{R}_{+}}}}(t),\]Lastly, let us recall the direct relation linking the tfPp with fractional calculus: the probability distribution of the tfPp solves a specific difference-differential equation in which the time derivative appearing in the difference-differential equations related to the homogeneous Poisson process is replaced by a fractional derivative of the Caputo–Džrbašjan type. Regarding this, for $n\in \mathbb{N}$, denote by ${\text{AC}^{m}}[a,b]$ the space of real-valued functions with continuous derivatives up to order $m-1$ such that the $(m-1)$-th derivative belongs to the space of absolutely continuous functions $\text{AC}[a,b]$. In other words,
Then, for $\alpha >0$, $m=\lceil \alpha \rceil $, and $f\in A{C^{m}}[a,b]$, the Caputo–Džrbašjan derivative of order $\alpha >0$ is defined as
Let us now denote the state probabilities $\mathbb{P}({\mathcal{N}_{t}}=k)$ of the tfPp by ${p_{k}}(t)$, $k\ge 0$, $t\ge 0$. Then the probabilities ${p_{k}}(t)$ satisfy the equations
where we consider ${p_{-1}}(t)$ being equal to zero. In Example 2.2, the tfPp will be compared with a member of the class of the mPp-utt.
(8)
\[ {\text{AC}^{m}}[a,b]=\bigg\{f:[a,b]\mapsto \mathbb{R}:\frac{{\text{d}^{m-1}}}{\text{d}{x^{m-1}}}f(x)\in \text{AC}[a,b]\bigg\}.\](9)
\[ {^{C}}{D_{{a^{+}}}^{\alpha }}f(t)=\frac{1}{\varGamma (m-\alpha )}{\int _{a}^{t}}{(t-s)^{m-1-\alpha }}\frac{{\text{d}^{m}}}{\text{d}{s^{m}}}f(s)\text{d}s.\](10)
\[ {^{C}}{D_{{0^{+}}}^{\alpha }}{p_{k}}(t)=-\lambda {p_{k}}(t)+\lambda {p_{k-1}}(t),\hspace{1em}k\ge 0,\]2.2 Mixed Poisson processes up to a time transformation (mPp-utt)
We start describing mixed Poisson processes (mPp), first introduced by J. Dubourdieu in 1938 [15]. For full details the reader can refer to the monograph by J. Grandell [20].
Consider a unit-rate homogeneous Poisson process $N={({N_{t}})}_{t\ge 0}$. A point process ${({\widetilde{M}_{t}})}_{t\ge 0}$ is an mPp if and only if ${\widetilde{M}_{t}}={N_{Wt}}$ in distribution, where W is an almost surely non-negative random variable independent of N. Common choices for the mixing random variable W are the Gamma distribution, leading to the Pólya process, or the uniform distribution on $[0,c)$, $c\in {\mathbb{R}_{+}}$.
Clearly, if W is degenerate on w, then an mPp coincides with a homogeneous Poisson process of rate w.
An mPp is characterized by the so-called property P which means that, conditional on ${\widetilde{M}_{t}}-{\widetilde{M}_{0}}=k$, the random jump times $\{{t_{1}},{t_{2}}\dots ,{t_{k}}\}$ are distributed as the order statistics of k iid uniform random variables on $[0,t]$ [17].
This result, first appeared in [32], can be further extended considering a deterministic time-change, leading to the class of point processes with the OS property (order statistics property) [10, 17]: a point process ${({K_{t}})}_{t\ge 0}$ with unit steps is said to have the OS property if and only if, conditional on ${K_{t}}-{K_{0}}=k$, the random jump times $\{{t_{1}},{t_{2}}\dots ,{t_{k}}\}$ are distributed as the order statistics of k iid random variables supported on $[0,t]$ with distribution function ${F_{t}}(x)=q(x)/q(t)$, where $q(t)=\mathbb{E}({K_{t}})-\mathbb{E}({K_{0}})$ is continuous and non-decreasing. In this respect, property P is also called uniform OS property.
Notably, K.S. Crump proved that point processes with the OS property are Markovian (see [10], Theorem 2).
Taking into account the results presented in [10, 17], and [39], we recall the following theorem due to P.D. Feigin [17]:
Theorem 2.1.
Let M be a point process with the $OS$ property relative to the distribution function ${F_{t}}(x)=q(x)/q(t)$, where $q(t)=\mathbb{E}({M_{t}})-\mathbb{E}({M_{0}})$ is a continuous and non-decreasing function. Then there exists a unit-rate homogeneous Poisson process N and an independent non-negative random variable W defined on the same probability space, such that ${M_{t}}={N_{Wq(t)}}$ almost surely.
Notice also that Theorem 2.1 implies $\mathbb{E}(W)=1$. Furthermore, the above theorem does not exclude the case of bounded q, that is when ${\lim \nolimits_{t\to \infty }}q(t)=\gamma <\infty $. Processes in that subclass are usually called mixed sample processes. To gain more insight on them, the reader can consult [39], Section 2, in which an interesting example is described (see also [11]).
It seems clear that the class of point processes with the OS property contains that of Mixed Poisson processes up to the time transformation $q(\cdot )$ (mPp-utt). We will consider in the following only the subclass of mPp-utt.
Let us first present a didactic example of a member of the class of mPp-utt, the Yule process. Being an mPp-utt, the Yule process exhibits the OS property. The reader may refer to [10] for more details.
Example 2.1.
Let M be a Yule process starting with a single individual, shifted downwards by one and with individual splitting rate λ. Let N be a unit-rate homogeneous Poisson process and let W be independent of N and exponentially distributed with mean one. Set $q(t)={e^{\lambda t}}-1$. Then, we construct the mPp-utt representation of M, i.e. ${M_{t}}={N_{Wq(t)}}$, $t\ge 0$. Note that the state space of ${N_{Wq(t)}}$ is $\{0,1,\dots \hspace{0.1667em}\}$. The distribution of M can be derived easily by conditioning:
The Yule process has the OS property with ${F_{t}}=\frac{{e^{\lambda x}}-1}{{e^{\lambda t}}-1}$, $x\in [0,t]$. The OS property of the Yule process has been implicitly used in many papers, starting from the seminal paper by Yule [45] (see also [9, 27, 28])
(11)
\[\begin{aligned}{}\mathbb{P}({M_{t}}=k)& ={\int _{0}^{\infty }}\frac{{[w({e^{\lambda t}}-1)]^{k}}}{k!}{e^{-w({e^{\lambda t}}-1)}}{e^{-w}}\text{d}w\\ {} & =\frac{{({e^{\lambda t}}-1)^{k}}}{k!}{\int _{0}^{\infty }}{w^{k}}{e^{-w{e^{\lambda t}}}}\text{d}w\\ {} & ={e^{-\lambda t}}{\big(1-{e^{-\lambda t}}\big)^{k}},\hspace{1em}k\ge 0.\end{aligned}\]In the following example we define a specific mPp-utt which is connected with the tfPp by means of its marginal distribution.
Example 2.2.
Let $M={({M_{t}})}_{t\ge 0}$ be such that ${M_{t}}={N_{Wq(t)}}$ where $q(t)=\lambda {t^{\nu }}/\varGamma (1+\nu )$, $\lambda >0$, $\nu \in (0,1)$, and W is a unit-mean non-negative random variable with probability density function
The above density is written in terms of the Wright function (see [22] for details)
(12)
\[ {f_{W}}(w)=\phi \bigg(-\nu ,1-\nu ;\frac{-w}{\varGamma (1+\nu )}\bigg)\bigg/\varGamma (1+\nu ),\hspace{1em}w\in {\mathbb{R}_{+}}.\]Plainly, $\mathbb{E}({M_{t}})=\lambda {t^{\nu }}/\varGamma (1+\nu )$. Now, let us derive the marginal distribution of the mPp-utt M:
By letting $\xi =w{t^{\nu }}/\varGamma (1+\nu )$ we have
This last step is justified by the time-change construction of the tfPp and by the fact that the marginal density function of the inverse stable subordinator ${({E_{t}})}_{t\ge 0}$ is exactly ${f_{{E_{t}}}}(\xi )={t^{-\nu }}\phi (-\nu ,1-\nu ;-\xi {t^{-\nu }})$, $\xi \in {\mathbb{R}_{+}}$ (see e.g. [14], Section 2).
(14)
\[\begin{aligned}{}\mathbb{P}({M_{t}}=k)& ={\int _{0}^{\infty }}\mathbb{P}({N_{wq(t)}}=k){f_{W}}(w)\text{d}w\\ {} & ={\int _{0}^{\infty }}\frac{{(w\frac{\lambda {t^{\nu }}}{\varGamma (1+\nu )})^{k}}}{k!}{e^{-w\frac{\lambda {t^{\nu }}}{\varGamma (1+\nu )}}}\phi \bigg(-\nu ,1-\nu ;\frac{-w}{\varGamma (1+\nu )}\bigg)\frac{\text{d}w}{\varGamma (1+\nu )}.\end{aligned}\](15)
\[\begin{aligned}{}\mathbb{P}({M_{t}}=k)& ={\int _{0}^{\infty }}\frac{{(\lambda \xi )^{k}}}{k!}{e^{-\lambda \xi }}{t^{-\nu }}\phi \big(-\nu ,1-\nu ;-\xi {t^{-\nu }}\big)\text{d}\xi \\ {} & =\mathbb{P}({\widetilde{N}_{{E_{t}}}}=k)=\mathbb{P}({\mathcal{N}_{t}}=k).\end{aligned}\]It is worthy of note that the tfPp $\mathcal{N}$ and the mPp-utt M share the same marginal distribution. This entails that the probabilities $\mathbb{P}({M_{t}}=k)$, $k\ge 0$, solve the equations (10).
3 Generalized Yule model
We proceed now to the analysis of a generalization of the Yule model in the sense we have anticipated in the introductory section. The focus here is to construct a model in which the arrival in time of genera is driven by an mPp-utt and the process describing the evolution of species for each different genus is a tfPp or a nonlinear time-fractional pure birth process (see [28, 33]). Hence what we drop here is the deterministic constant intensity assumption for genera evolution. The genera arrival is instead described by random intensities. For the sake of clarity, before describing the generalization of the Yule model, let us recall the definition of the nonlinear time-fractional pure birth process. Analogously as for the tfPp, the construction of the nonlinear time-fractional pure birth process is by time-change with the inverse stable subordinator ${({E_{t}})}_{t\ge 0}$. Thus, consider the nonlinear pure birth process ${({Y_{t}})}_{t\ge 0}$, starting with a single progenitor, with nonlinear rates ${\lambda _{k}}>0$, $k\ge 1$, and being independent of ${({E_{t}})}_{t\ge 0}$. The time-changed process $\mathcal{Y}={({\mathcal{Y}_{t}})}_{t\ge 0}={({Y_{{E_{t}}}})}_{t\ge 0}$ is called a nonlinear time-fractional pure birth process. Interestingly enough, when ${\lambda _{k}}=\lambda $ for each $k\ge 1$, the nonlinear time-fractional pure birth process reduces to the tfPp of parameters λ and ν, shifted upwards by one.
Let us now consider the following model.
Definition 3.1 (Generalized Yule model).
The generalized Yule model represents the growth of a population which evolves according to:
Then, for each time $t\in {\mathbb{R}_{+}}$ we define the random variable ${_{t}}\mathfrak{N}$ measuring the number of species belonging to a genus chosen uniformly at random. With respect to the classical Yule model this random variable is linked to the degree distribution of a vertex chosen uniformly at random in the Barabási–Albert model [36]. Our aim is to investigate the distribution of ${_{t}}\mathfrak{N}$ for the generalized Yule model. To do so, it is enough to condition on the random creation time T of the selected genus, obtaining
Notice that, due to the $OS$ property satisfied by the considered mPp-utt, the distribution function of T is ${F_{t}}(\cdot )$ (see Section 2.2).
(16)
\[ \mathbb{P}({_{t}}\mathfrak{N}=k)={\mathbb{E}_{T}}\mathbb{P}({\mathcal{Y}_{t}}=k|{\mathcal{Y}_{T}}=1),\hspace{1em}k\ge 1.\]We specialize now the model by choosing the process of Example 2.2 for the random arrival of genera.
In this case the distribution function of T reads
with density
Figure 1 shows the shapes of the distribution function (17) and the density function (18) for different values of the characterizing parameter ν. Notice the rather different behaviour for values of ν strictly less than 1.
Fig. 1.
Distribution function (17) (top) and density function (18) in linear plot (middle) and loglog plot (bottom). The parameter ν is set to $\nu =(1/4,1/2,3/4,1)=(\text{blue},\text{orange},\text{green},\text{red})$ and $t=1$. Note how, for $\nu \in (0,1)$ (in contrast with the classical case $\nu =1$), the concentration of probability mass near zero makes the appearance of genera more likely to occur in the very early evolution of the process
Regarding the evolution of the number of species for each genus, the fractional exponent of the process $\mathcal{Y}$ will be denoted by $\beta \in (0,1)$. We suppose now that the nonlinear rates of $\mathcal{Y}$ are all different and recall that in this case
with the convention that empty products equal unity. We obtain
Now we make use of Corollary 2.3 of [21] and arrive at
(19)
\[ \mathbb{P}({\mathcal{Y}_{t}}=k)={\prod \limits_{j=1}^{k-1}}{\lambda _{j}}{\sum \limits_{m=1}^{k}}\frac{{E_{\beta }}(-{\lambda _{m}}{t^{\beta }})}{{\textstyle\prod _{l=1,l\ne m}^{k}}({\lambda _{l}}-{\lambda _{m}})},\hspace{1em}k\ge 1,\](20)
\[\begin{aligned}{}\mathbb{P}({_{t}}\mathfrak{N}=k)& ={\int _{0}^{t}}\mathbb{P}({\mathcal{Y}_{t-x}}=k)\frac{\nu {x^{\nu -1}}}{{t^{\nu }}}\text{d}x\\ {} & =\frac{\nu }{{t^{\nu }}}{\prod \limits_{j=1}^{k-1}}{\lambda _{j}}{\sum \limits_{m=1}^{k}}\frac{1}{{\textstyle\prod _{l=1,l\ne m}^{k}}({\lambda _{l}}-{\lambda _{m}})}{\int _{0}^{t}}{E_{\beta }}\big[-{\lambda _{m}}(t-x)\big]{x^{\nu -1}}\text{d}x.\end{aligned}\]A special case of interest is when the rates are linear, ${\lambda _{k}}=\lambda k$, $k\ge 1$. In this case, from (21) we obtain easily the following probabilities:
Figure 2 shows how the above probability mass function changes with respect to parameter ν, taking a constant $\beta =1$, that is, considering a classical behaviour for species.
(22)
\[ \mathbb{P}({_{t}}\mathfrak{N}=k)=\varGamma (\nu +1){\sum \limits_{j=1}^{k}}\left(\genfrac{}{}{0.0pt}{}{k-1}{j-1}\right){(-1)^{j-1}}{E_{\beta ,\nu +1}}\big(-\lambda j{t^{\beta }}\big),\hspace{1em}k\ge 1.\]Recalling that in the linear rates case $\mathbb{E}{\mathcal{Y}_{t}}={E_{\beta }}(\lambda {t^{\beta }})$ and $\mathbb{E}{\mathcal{Y}_{t}^{2}}=2{E_{\beta }}(2\lambda {t^{\beta }})-{E_{\beta }}(\lambda {t^{\beta }})$, we derive the first two moments for the random variable ${_{t}}\mathfrak{N}$ and its variance:
(23)
\[\begin{aligned}{}\mathbb{E}\hspace{0.2778em}{_{t}}\mathfrak{N}& =\frac{\nu }{{t^{\nu }}}{\int _{0}^{t}}{E_{\beta }}\big[\lambda {(t-x)^{\beta }}\big]{x^{\nu -1}}\mathrm{d}x=\varGamma (\nu +1){E_{\beta ,\nu +1}}\big(\lambda {t^{\beta }}\big),\end{aligned}\](24)
\[\begin{aligned}{}\mathbb{E}\hspace{0.2778em}{_{t}}{\mathfrak{N}^{2}}& =\frac{\nu }{{t^{\nu }}}{\int _{0}^{t}}\mathbb{E}{\mathcal{Y}_{t-x}^{2}}{x^{\nu -1}}\mathrm{d}x\\ {} & =2\varGamma (\nu +1){E_{\beta ,\nu +1}}\big(2\lambda {t^{\beta }}\big)-\varGamma (\nu +1){E_{\beta ,\nu +1}}\big(\lambda {t^{\beta }}\big),\end{aligned}\](25)
\[\begin{aligned}{}\mathbb{V}\text{ar}\hspace{0.2778em}{_{t}}\mathfrak{N}& =2\varGamma (\nu +1){E_{\beta ,\nu +1}}\big(2\lambda {t^{\beta }}\big)\\ {} & \hspace{1em}-\varGamma (\nu +1){E_{\beta ,\nu +1}}\big(\lambda {t^{\beta }}\big)\big(1+\varGamma (\nu +1){E_{\beta ,\nu +1}}\big(\lambda {t^{\beta }}\big)\big).\end{aligned}\]Fig. 2.
The probability mass function of ${_{t}}\mathfrak{N}$ in the linear rate case (formula (22)). Top: linear plot, bottom: loglog plot, with $t=10$, $\beta =1$ (which corresponds to classical linear birth processes for species evolution), $\lambda =1$, and $\nu =(0.4,0.6,0.8,1)$ from bottom to top
When the nonlinear rates are actually constant and all equal (i.e. ${\lambda _{k}}=\lambda $, $\forall k\ge 1$) we cannot make use of a specialized form of formula (19). In this case however, the nonlinear time-fractional pure birth process reduces to the tfPp suitably shifted upwards by one. Hence, recalling formula (2), the distribution of ${_{t}}\mathfrak{N}$ can be written by conditioning as
Then we have
The above integral is known and can be calculated by using Corollary 2.3 of [21]. We finally obtain
(26)
\[\begin{aligned}{}\mathbb{P}({_{t}}\mathfrak{N}=k)& ={\mathbb{E}_{T}}\mathbb{P}({\mathcal{N}_{t}}+1=k|{\mathcal{N}_{T}}+1=1)\\ {} & ={\int _{0}^{t}}\mathbb{P}({\mathcal{N}_{t-x}}=k-1)\hspace{0.2778em}\frac{\nu {x^{\nu -1}}}{{t^{\nu }}}\text{d}x,\hspace{1em}k\ge 1.\end{aligned}\](27)
\[\begin{aligned}{}\mathbb{P}({_{t}}\mathfrak{N}=k)& ={\int _{0}^{t}}{\big[\lambda {(t-x)^{\beta }}\big]^{k-1}}{E_{\beta ,\beta (k-1)+1}^{k}}\big[-\lambda {(t-x)^{\beta }}\big]\hspace{0.1667em}\frac{\nu {x^{\nu -1}}}{{t^{\nu }}}\text{d}x\\ {} & =\frac{\nu {\lambda ^{k-1}}}{{t^{\nu }}}{\int _{0}^{t}}{(t-x)^{\beta k-\beta }}{E_{\beta ,\beta k-\beta +1}^{k}}\big[-\lambda {(t-x)^{\beta }}\big]{x^{\nu -1}}\text{d}x.\end{aligned}\]In the classical Yule model, ${\lim \nolimits_{t\to \infty }}{_{t}}\mathfrak{N}=\mathfrak{N}$, where $\mathfrak{N}$ is a non-degenerate limiting random variable. The distribution of $\mathfrak{N}$ is known and is called Yule–Simon distribution. Its main feature is the characteristic right tail which slowly decays as a power-law. In our cases, however, the random variable ${_{t}}\mathfrak{N}$ has a different behaviour at ∞. This can be observed by considering the asymptotic expansion of the Prabhakar function [19]. We have from formula (28), for $t\to \infty $,
for each finite value of $k\ge 1$.
(29)
\[ \mathbb{P}({_{t}}\mathfrak{N}=k)\sim \frac{\varGamma (\nu +1)}{\lambda }\frac{{t^{-\beta }}}{\varGamma (-\beta k)}\longrightarrow 0,\]3.1 A critical macroevolutionary model with species deletion
We introduce here a model of macroevolution in which the possibility of extinction of genera is taken into consideration. To achieve this, the species dynamics is described by independent critical birth-death processes (of parameter $\lambda >0$), each starting with a single species, while the genera appearance follows the mPp-utt of Example 2.2. Figure 3 shows a possible realization of the superimposed processes counting the number of species belonging to each existent genus. Extinction of genera is represented by squares while their births by circles.
Fig. 3.
A possible realization of the superimposed processes counting the number of species belonging to each existent genus (shown in different colours). Birth of genera, governed by the mPp-utt of Example 2.2, is represented in the figure by circles while their extinction time is marked by squares
Recurring to the integral representation of the Gauss hypergeometric function,
for $c>b>0$, we derive the exact form of the transient distribution of ${_{t}}\mathfrak{N}$:
and
(30)
\[ {_{2}}{F_{1}}(a,b;c;z)=\frac{\varGamma (c)}{\varGamma (b)\varGamma (c-b)}{\int _{0}^{1}}{y^{b-1}}{(1-y)^{c-b-1}}{(1-yz)^{-a}}\mathrm{d}y,\](31)
\[\begin{aligned}{}\mathbb{P}({_{t}}\mathfrak{N}=0)& =\frac{\nu }{{t^{\nu }}}{\int _{0}^{t}}{x^{\nu -1}}\frac{\lambda (t-x)}{1+\lambda (t-x)}\mathrm{d}x\\ {} & =\frac{\lambda t}{\nu +1}{_{2}}{F_{1}}(1,2;\nu +2;-\lambda t),\end{aligned}\](32)
\[\begin{aligned}{}\mathbb{P}({_{t}}\mathfrak{N}=k)& =\frac{\nu }{{t^{\nu }}}{\int _{0}^{t}}{x^{\nu -1}}\frac{{\lambda ^{k-1}}{(t-x)^{k-1}}}{{(1+\lambda (t-x))^{k+1}}}\mathrm{d}x\\ {} & ={(\lambda t)^{k-1}}\frac{\varGamma (k)\varGamma (\nu +1)}{\varGamma (\nu +k)}{_{2}}{F_{1}}(k+1,k;\nu +k;-\lambda t),\hspace{1em}k\ge 1.\end{aligned}\]Fig. 5.
The probability mass function (32) depicted for $\nu =(0.1,0.5,1)=(\text{blue},\text{orange},\text{green})$ with $\lambda =1$
In Figure 4 the above probabilities (for $k=20$) are pictured with respect to time. The probability mass function concentrates on zero for $t\to \infty $ as it should be. Notably, it exhibits an exponential tail (see in Figure 5), differently from the case without deletion (for example compare it with (22), see Figure 2). The derivation of the moments of the random variable ${_{t}}\mathfrak{N}$ is simpler in this model. Recalling that each species process has mean 1 and $\mathbb{E}{\mathcal{Y}_{t}^{2}}=2\lambda t+1$ we obtain that the expectation of ${_{t}}\mathfrak{N}$ is also 1 and that
