VMSTA Modern Stochastics: Theory and Applications 2351-6054 2351-6046 2351-6046 VTeXMokslininkų g. 2A, 08412 Vilnius, Lithuania VMSTA163 10.15559/20-VMSTA163 Research Article Geometric branching reproduction Markov processes TchorbadjieffAssenatchorbadjieff@math.bas.bg MaysterPenkapenka.mayster@math.bas.bg∗ Institute of Mathematics and Informatics, Bulgarian Academy of Sciences, Acad. G. Bonchev street, Bloc 8, 1113 Sofia, Bulgaria Corresponding author. 2020 309202074357378 762020 2092020 2092020 © 2020 The Author(s). Published by VTeX2020 Open access article under the CC BY license.

We present a model of a continuous-time Markov branching process with the infinitesimal generating function defined by the geometric probability distribution. It is proved that the solution of the backward Kolmogorov equation is expressed by the composition of special functions – Wright function in the subcritical case and Lambert-W function in the critical case. We found the explicit form of conditional limit distribution in the subcritical branching reproduction. In the critical case, the extinction probability and probability mass function are expressed as a series containing Bell polynomial, Stirling numbers, and Lah numbers.

 Branching process Lagrange inversion Gauss hypergeometric Wright Lambert-W functions extinction probability 11B73 33C05 60J80 Bulgarian National Science Fund KP-06-H22/3 This research has been partially supported by the Bulgarian National Science Fund, grant No KP-06-H22/3. Introduction

Any time-homogeneous Markov branching process (MBP) X(t),t>0 0$]]>, is defined and studied by its probability generating function (p.g.f.) F(t,s),|s|<1 , as an unique solution of the Kolmogorov equations [1, 7, 19]. The backward Kolmogorov equation is nonlinear, of the type separate differentials, due to the time-homogeneous property. First of all, we find the indefinite integral primitive. Then, in order to respect the initial conditions, it is essential to define the inverse function of the primitive. We apply the Lagrange inversion method as it is shown in many books and articles [3, 8, 12].

The correspondence between MBP and integer-valued positive Lévy processes is noted in the books [9], page 236, and [1], page 126. The geometric distribution is present in many application models, mainly in biology and epidemiology. However, the most important implementation is in statistical physics, where it is known as Bose-Einstein distribution [20].

In the present text, we are focused on the subcritical and critical infinitesimal geometric branching reproductions. We prove that the p.g.f. F(t,s) in the subcritical case is expressed as a composition of Gauss hypergeometric and Wright functions. Many special cases of the Wright function 1Ψ1(α,a;β,b;z) are considered in [11, 16]. The numerical evaluation of the Wright function is studied in [10]. The Laplace transform pairs related to the Mittag-Leffler functions and the reduced Wright functions is surveyed in [11]. We remark that in our model, for the subcritical MBP with mean of reproduction 0<m<1 , the parameters of the Wright function are α=a=m , and −1<β=m−1<0 , and b=m+1 . The principal parameter, defining the domain of convergence is equal to 1+β−α=1+(m−1)−m=0 [16]. The factorial moments of X(t),t>0 0$]]>, are presented by series containing the partial Bell polynomials [21]. The conditional limit distribution is defined in its explicit form. It is a new unimodal integer-valued distribution supported by (1,2,…) . Its index of dispersion depends on the solution of a transcendental equation.

In the critical case, the p.g.f. F(t,s) is defined by the composition of the Lambert-W function and linear-fractional one. The probability mass function (p.m.f.) of X(t) , t>0 0$]]>, is expressed by the values of the Lambert-W function at the point x=eKt+1 and can be calculated by using any corresponding software packages. The agreement between Lagrange inversion method and series expansion of the Lambert-W function is confirmed at the approximation of the extinction probability.

The Lambert-W function is considered as the real-valued composite inverse function of the function V(x)=xex,x≥−1 . The properties and many applications of the Lambert-W function are provided in [5]. The sequence of series and derivatives of W(x) are studied in [4]. The new probability property such that W(x) is a Bernstein function is developed in [15]. The computation of Lambert W-functions without transcendental function evaluations is considered in [6]. A historical review is given in [2].

The main convenience of the obtained results based on special functions is computation simplicity. It is not a novelty and there are many applications. The most important ones are related to the multi-body computational problems, firstly noticed in the numeric calculation of the eigenstates of hydrogen molecular ion H+ . The effective solution was obtained by using the Lambert-W function, which becomes a common method for solving similar problems in Physics [13]. The best advantage of using Lambert-W is the effective approach to the solution of slowly convergent series expansions [17]. It is possible due to very well studied properties and well-developed software packages, noted in [5, 6, 22].

We began our work on the infinitesimal geometric branching mechanism following the methods of branching processes theory and Lagrange inversion. Using them, we extended our results to obtain a computational effective generalisation based on Lambert-W and Wright functions. Some of the improvements are demonstrated by easily computed numerical applications.

 Subcritical geometric branching

The studied geometric branching reproduction process X(t),t>0 0$]]>, is a time-homogeneous MBP starting with one particle as initial condition. Its Markov property is guaranteed by the assumption that the lifetime of particles is an exponentially distributed random variable with a constant parameter K>0 0$]]>. The reproduction is defined by the integer-valued random variable η representing the offspring numbers as follows, P(η=k)=mk(1+m)k+1,k=0,1,…,m>0. 0.\]]]> The parameter m>0 0$]]> defines the mean of the offspring numbers. In this notation the p.g.f. of the reproduction law is, h(s)=11+m−ms,h′(s)=m(1+m−ms)2,h″(s)=2m2(1+m−ms)3. The infinitesimal generating function f(s) present in the Kolmogorov equations is defined by the constant K and p.g.f. h(s) as f(s)=K(h(s)−s)=K11+m−ms−s=K(1−s)(1−ms)1+m(1−s), and has the following derivatives, for k=1,2,3,… , f′(s)=f′(1)−Km(1−s)(2+m−ms)(1+m−ms)2,f(k)(s)=Kmkk!(1+m−ms)k+1. The p.g.f. of the number of particles alive at the time t>0 0$]]> F(t,s)= ∑k=0∞skP(X(t)=k|X(0)=1) yields the backward Kolmogorov equation ∂∂tF(t,s)=fF(t,s),F(0,s)=s. The equation h(s)=s has two solutions, s1=1/m and s2=1 , where m=h′(1) . The value s=1/m is the fixed point for the p.g.f. h(s) and for its first derivative h′(s) . The branching reproduction is classified as subcritical if 0<m<1 and critical if m=1 . This classification is in accordance with the notion of ultimate extinction probability q=1 in both cases in consideration, subcritical and critical one. In particular, for the subcritical reproduction when 0<m<1 , the first derivative f′(s) is negative in the interval 0≤s≤1 , and f′(0)=−K(1+m+m2)(1+m)2<0,f′(1)=K(m−1)<0. All derivatives for k=2,3,… , at the points s=0 and s=1 are related by the constant (1+m)k+1 as follows, f(k)(0)=Kk!mk(1+m)k+1,f(k)(1)=Kk!mk,f(k)(0)=f(k)(1)(1+m)k+1. In the subcritical case, we study in parallel the p.g.f. F(t,s) of the branching process X(t),t>0 0$]]>, defined by (2) and equation (3), and the function R(t,s)=1−F(t,s), satisfying the following non-linear equation: ∂R(t,s)∂t=−f(1−R(t,s)),R(0,s)=1−s. In order to solve these equations (3) and (6), when f′(1)<0 , we define respectively the function A0(s)=expf′(1)∫0sdxf(x),s∈[0,1),A0(0)=1,A0(1)=0. and the function A(s)=expf′(1)∫0s−dxf(1−x),s∈[0,1),A(0)=0,A(1)=1, where A0′(s)A0(s)=f′(1)f(s),A′(s)A(s)=−f′(1)f(1−s), see [14, 1]. We remark that in the subcritical case 1f(s)=−1f′(1)11−s−m21−ms,f′(1)=K(m−1)<0. The solution of the indefinite integral is ∫dxf(x)=∫−1f′(1)11−x−m21−mxdx=1f′(1)log|1−x||1−mx|m. Knowing the primitive (8) of the equation (3) we formulate the following lemma without proof.

Let X(t),t>0 0$]]>, be a subcritical time-homogeneous MBP with branching mechanism given by (1) where 0<m<1 . The solution of the backward Kolmogorov equation (3) satisfying F(t,1)=1 , is F(t,s)=A0−1(ef′(1)tA0(s)),A0(s)=1−s(1−ms)m,s<1/m, respectively, the solution of the equation (6) with R(t,1)=0 , is R(t,s)=A−1(ef′(1)tA(1−s)),A(s)=s(1−m+ms)m,1−1/m<s. where A0−1(x) and A−1(x) are the composite inverse functions, and A0(s)=A(1−s),A0−1(x)=1−A−1(x),−∞<x<∞.

The mathematical expectation of the number of particles alive at the time t>0 0$]]> in the subcritical case has the following exponential decreasing behaviour, E[X(t)]=exp{f′(1)t}=eMt<1,M=f′(1)=K(m−1)<0,m=h′(1). The extinction probability to the positive time P(X(t)=0)=F(t,0),t>0 0$]]>, is expressed by the inversion function F(t,0)=A0−1(eMtA0(0))=A0−1(eMt),eMt<1,A0(0)=1. Respectively, the survival probability is given as follows, R(t,0)=A−1(eMtA(1))=A−1(eMt)=1−A0−1(eMt). We remark that at the point s=1/2 the functions A0(1/2)=A(1/2) and the respective vertical asymptotes are left and right symmetric to s=1/2 . The function A0(s) is decreasing, concave and all its derivatives A0(k)(s) are negative in the interval, s<1/m . The function A(s) is increasing, concave, the consecutive derivatives are alternating positive and negative in the interval, s>1−1/m 1-1/m$]]>. Namely, A0(2k)(s)=A(2k)(1−s)<0,A0(2k−1)(s)=−A(2k−1)(1−s)<0. In details, denoting the falling and rising factorials as follows: [x]n↓=x(x−1)...(x−n+1)=Γ(x+1)Γ(x+1−n), and [x]n↑=x(x+1)...(x+n−1)=Γ(x+n)Γ(x), we write the derivatives in the following form, A(k)(s)=(−1)k−1[m](k−1)↑θk−1(k+ms)(1−m)m(1+θs)m+k,θ=m1−m>0. 0.\]]]> In particular, for k=1,2,… , at the points s=0 and s=1 we have ak=A(k)(0)=(−1)k−1k[m](k−1)↑θk−1(1−m)m=(−1)kA0(k)(1), and A(k)(1)=(−1)k−1[m](k−1)↑θk−1(k+m)(1−m)m(1+θ)m+k=(−1)kA0(k)(0). Obviously, the following relation follows from (11) and (12) A(k)(1)=A(k)(0)k+mk(1+θ)m+k. The range of the function A(s) is the whole real line (−∞,∞) . It is considered as the domain of definition for the A−1(x) . Knowing the derivatives (11), we find the function A−1(x) in its representation as an exponential generating function in the interval |x|<1/mm+1,0<m<1 .

Lagrange inversion and special functions

Now, we look for the properties of the composite inverse of the functions A0(s) and A(s) . The function A(s) and its inverse A−1(s) are some particular cases of the Wright function [11, 16] 1Ψ1(α,a;β,b;z):= ∑n=0∞Γ(αn+a)Γ(βn+b)znn!, and the Gauss hypergeometric function [8, 3] 2F1(c,d;g;z):= ∑k=0∞[c]k↑[d]k↑[g]k↑zkk!,2F1c,d;d;z=11−zc,|z|<1. Obviously, knowing the notation θ=m1−m>0,0<m<1 0,0, we see that the A(s)=s(1−m)m(1+θs)m=2F1(m,1;1;−θs)s(1−m)m,|θs|<1. Using the previous definitions, it will be proved in the following theorem that for |x|<1/mm+1 and 0<m<1 the inverse function A−1(x)=1Ψ1m,m;m−1,m+1;mx(1−m)1−mm(1−m)x.

Let us consider the function A0(s) given by A0(s)=A(1−s)=1−s(1−ms)m,0<m<1,A0(0)=1,A0(1)=0. Then for its composite inverse function is valid A0−1(x)=1−A−1(x)=1− ∑k=1∞bkxkk!,|x|<1/mm+1, where bk=(1−m)mkm1−mk−1[mk](k−1)↓.

First of all, we confirm the relation A0−1(A0(s))=1−A−1(A0(s))=1−A−1(A(1−s))=1−(1−s)=s. The function A(s) is represented as an exponential generating function in the neighbourhood of zero |s|<1/θ=(1−m)/m,0<m<1 (13), with coefficients ak=A(k)(0) (11) as follows A(s)= ∑k=1∞k(−1)k−1[m](k−1)↑(1−m)mm1−mk−1skk!,a1=1(1−m)m>1. 1.\]]]> To apply the Lagrange inversion method to it, we use the following representation, A(s)=sg(s),g(s)=(1−m)m1+ms1−mm,g(0)=(1−m)m>0, 0,\]]]> and obviously g(1)=1 , see [3, 8]. Denote the inverse function in its series expansion as A−1(x)= ∑k=1∞bkxkk!. Then the coefficients bk are given by the derivatives of the function (g(s))k at the point s=0 as follows, bk=dk−1dsk−1(g(s))ks=0. The Taylor series expansion of the function (g(s))k=(1−m)mk1+ms1−mmk,mk>0, 0,\]]]> is given by the binomial coefficients as follows, (g(s))k=(1−m)mk ∑j=0∞ms1−mj[mk]j↓j!. The derivative of the order j to this function at the point s=0 is djdsj[(g(s))k](s=0)=(1−m)mkm1−mj[mk]j↓. It is enough to take j=k−1 in order to obtain the coefficient bk in (15). Then A−1(x)= ∑k=1∞(1−m)mkm1−mk−1[mk](k−1)↓xkk!,b1=(1−m)m<1. Applying the definition of decreasing factorials by Gamma function (9, 10), Γ(mk+1)Γ(mk+1−k+1)=mkΓ(mk)Γ((m−1)k+2)=mkΓ(m(k−1)+m)Γ((m−1)(k−1)+m+1) we find A−1(x)=m(1−m)x ∑k=1∞m(1−m)1−mk−1Γ(m(k−1)+m)Γ((m−1)(k−1)+m+1)xk−1(k−1)!. The change of variable k−1=j leads to the representation (14).

We remark that the coefficients bk are either positive, zero or negative, due to the decreasing factorials. It is confirmed, when the reflection formulas for Gamma function Γ(z)Γ(−z)=−πzsin(πz),Γ(z)Γ(1−z)=πsin(πz),z≠0,±1,±2,…, are applied on the representation [mk](k−1)↓=Γ(mk+1)Γ(1+(1−(1−m)k))=Γ(mk+1)Γ((1−m)k)sin(π(1−m)k)π(1−(1−m)k). Obviously, for k(1−m)=j,j=2,3,… , the value bk=0 . By this reason, in the considered case, the function A−1(x) is not a p.g.f..

The radius of convergence of the series expansion A−1(x)= ∑k=1∞bkxkk!,bk=m(1−m)1−mk(1−m)[mk](k−1)↓m, is calculated with a root test based on lim supk→∞|bk|k!k=m(1−m)1−mlim supk→∞(1−m)|[mk](k−1)↓|mk!k and Stiling’s formula for the Gamma function. When k→∞ we have Γ(mk+1)∼2πmkmkemk, and k!∼2πkkek,Γ((1−m)k)∼2π(1−m)k(1−m)k(1−m)ke(1−m)k. We extract the multiple mkem(1−m)ke(1−m)ekk={mm(1−m)1−m}k. Then, obviously, lim supk→∞2π(1−m)mk2k=1,lim supk→∞|sin(π(1−m)k)|k=1. And finally, lim supk→∞|bk|k!k=m(1−m)1−m{mm(1−m)1−m}=m1+m. We obtain the radius of convergence of the series expansion to the inverse function A−1(x) in dependence of the parameter m as follows, R(A−1)(m)=1m1+m>1,0<m<1. 1,\hspace{1em}0  □

In particular, when the parameter 0<m<1 is a rational number then the inverse function can be calculated as a solution of corresponding algebraic equation. For example, when m=1/2 then A(s)=2s1+s,A−1(x)=x(x+x2+8)4, and if m=1/3 then A(s)=33s2+s3,A−1(x)=x331+1−x3813+1−1−x3813, and if m=2/3 then A(s)=93s(1+2s)23, and A−1(x)=x183237x6+16x327+1+235x3+13+x183237x6+16x327+1−235x3+13+4x327. Respectively, their radiuses of convergence are as follows R(A−1)(1/2)=22,R(A−1)(1/3)=333,R(A−1)(2/3)=32943.

In the aim to describe the behaviour of the process X(t),t>0 0$]]>, we consider the factorial moments E[X(t)]n↓=Fs(n)(t,1) . The derivatives of the p.g.f. F(t,s) are expressed by the partial Bell polynomials defined as Bn,k(a∙)=∑(k1,k2,…,kn)n!a1k1...anknk1!(1!)k1...kn!(n!)kn,a∙=(a1,a2,…), see [21]. The sum is over all partitions of n into k parts, that is over all nonnegative integer solutions (k1,k2,…,kn) of the equations: k1+2k2+⋯+nkn=n,k1+k2+⋯+kn=k. The following theorem is based on the inequality (17) showing that the convergence interval of the series ∑k=1∞bksk/k! contains the point s=1 .

Let X(t),t>0 0$]]>, be a subcritical time-homogeneous MBP with branching mechanism given by (1) where 0<m<1 . Then the factorial moments E[X(t)]n↓=Fs(n)(t,1) are expressed by the mathematical expectation eMt=E[X(t)]=Fs′(t,1) and partial Bell polynomials Bn,k(a∙) as follows, Fs(n)(t,1)=(−1)n+1 ∑k=1nBn,k(a∙)bk(eMt)k,M<0, where the sequences a∙=(a1,a2,…) and b∙=(b1,b2,…) are defined as (11) and (15).

As we have seen previously, the solution of the equation (6) is 1−F(t,s)=A−1(eMtA(1−s))= ∑k=1∞bk(eMt)k(A(1−s))kk!. The change of variable z=1−s,s→1,z→0 , gives the opportunity to work with the representation (16) of A as the exponential generating function in the neighbourhood of zero. The powers of the exponential generating functions, see [21], are given by (A(z))kk!=1k! ∑j=1∞ajzjj!k= ∑n=k∞Bn,k(a∙)znn!. Then, exchanging the order of summation we obtain from (18), 1−F(t,s)= ∑k=1∞bk(eMt)k ∑n=k∞Bn,k(a∙)(1−s)nn!= ∑n=1∞ ∑k=1nBn,k(a∙)bk(eMt)k(1−s)nn!.  □

Note, that applying the Faa Di Bruno formula for derivatives of the composite functions A(A−1(x))=x and A−1(A(s))=s gives B1,1(a∙)b1=1, ∑k=1nBn,k(a∙)bk=0,n=2,3,… It means that all factorial moments Fs(k)(t,1) of order k=2,3,… contain the multiple eMt(1−eMt) .

Finally, the important process behaviour like conditional limit probability is obtained in the following theorem.

Let X(t),t>0 0$]]>, be a subcritical time-homogeneous MBP with branching mechanism given by (1) where 0<m<1 . Then the conditional limit probability exists and is given by limt→∞P(X(t)=n|X(t>0))=pn=P(ξ=n),n=1,2,…, 0))={p_{n}}=P(\xi =n),\hspace{1em}n=1,2,\dots ,\]]]> with p.g.f. F∗(s)= ∑n=1∞pnsn=1−1(1−ms)m+s(1−ms)m. The p.m.f. is defined by the increasing factorials (10), [m](k−1)↑ as follows pk=mk−1(1−m)(m+k)[m](k−1)↑k!>0,0<m<1,k=1,2,…. 0,\hspace{1em}0 The factorial moments E[ξ]k↓ of the limit random variable is fk=dkF∗dsk(1)=k[m](k−1)↑(1−m)mm1−mk−1>0,0<m<1,k=1,2,…. 0,\hspace{1em}0

The p.g.f. of the conditional probability is written as follows ∑n=1∞snP(X(t)=n|X(t)>0)=F(t,s)−F(t,0)1−F(t,0)=1−R(t,s)R(t,0). 0)=\frac{F(t,s)-F(t,0)}{1-F(t,0)}=1-\frac{R(t,s)}{R(t,0)}.\]]]> Knowing the differentiability of the function A−1 in the neighborhood of zero, we apply the theorem of l’Hospital to find the limit limt→∞R(t,s)R(t,0)=limt→∞A−1(eMtA(1−s))A−1(eMt)=A(1−s)=A0(s)=1−s(1−ms)m. To obtain the p.m.f. we consider the Taylor series expansion of the p.g.f. F∗(s) in the neighbourhood of zero, F∗(s)=1−A0(s)=1− ∑k=0∞[m]k↑(ms)kk!+s ∑j=0∞[m]j↑(ms)jj!. The change of variable j+1=k in the second sum leads to, F∗(s)= ∑k=1∞−[m]k↑(m)kk!+[m](k−1)↑(m)k−1(k−1)!sk= ∑k=1∞[m](k−1)↑(m)k−1(k−1)!1−m(m+k−1)ksk=(1−m) ∑k=1∞mk−1sk(m+k)[m](k−1)↑k!. We remark that 1−m(m+k−1)k=k−m2−mk+mk=m(1−m)+k(1−m)k. The relation between derivatives of the functions A and A0 at the points zero and s=1 is obvious. In this way, we recognise the factorial moments fk as given previously in (11), and also their relation with the p.m.f. pk in (12).  □

In a more general view, the conditional limit distribution for the studied process can be also considered as a mixture of shifted Negative-Binomial distributions, closely connected to the family of delta Lagrangian probability distribution [3, 8]. The probabilities pk of this process can be computed numerically after direct implementation of the Theorem 3. The computation process is recurrent and relies on the following ratio, where for k=1,2,… ., Υk(m)=pk+1pk=m1+mk+11−1k+m<1,0<m<1. This ratio (19) shows that the Panjer recursion (even generalized) does not take place. The inequality Υk(m)<1 is equivalent to the following k2(1−m)+k(m+1−2m2)+m(2−m2)>0,k=1,2,…, 0,\hspace{1em}k=1,2,\dots ,\]]]> where (m+1−2m2)>0 0$]]> when −1/2<m<1 . It is important to note that the ratio (19) is increasing as a function of the parameter m, namely Υk(m)<Υk(m+ε),0<m<m+ε<1. The histogram of the p.m.f. can be generated by the consecutive multiplication with Υk(m)<1 , where the maximum is given by p1=1−m2 . Some results for Υk(m) and computed from them probabilities pk are shown in Figure 1.

Mean-Variance-Skewness-Kurtosis can be easily expressed by the factorial moments. In particular, for several factorial moments we have, E[ξ]=f1=1(1−m)m,f2=2m2(1−m)m+1,f3=3m3(m+1)(1−m)m+2. The central moments, respectively, are as follows Var[ξ]=f2+f1−(f1)2=1(1−m)m2m21−m+1−1(1−m)m, and E[ξ−f1]3=f3+3f2+f1+2(f1)3−3f1(f1+f2). The main parameter for applications is the variance-to-mean ratio (VMR) or index of dispersion equal respectively to Var[ξ]E[ξ]=2m21−m+1−1(1−m)m,0<m<1. The measure VMR of dispersion obtains the value equal to one for the threshold parameter m∗ , which is estimated as 0.58905<m∗<0.589058 . This result is obtained after the numeric solution of the equation 2m2=(1−m)1−m, with a precision higher than 10−5 . The conditional limit probability is over-dispersed when the branching process X(t),t>0 0$]]>, is near critical, m∗<m<1 , and under-dispersed when 0<m<m∗<1 . The values of VMR are shown in Figure 2.

The graphics shows Υk(m) (left) and related values for pk (right) computed for m=1/3 (red), m=1/2 (green) and m=2/3 (blue)

Index of Dispersion (VMR)

In the class of Panjer probability distributions, the VMR =1 for the Poisson distribution, respectively, VMR >1 1$]]> for the geometric and Negative-Binomial distribution, and VMR <1 for the binomial distribution. This enables the index of dispersion to assess whether observed data can be modeled using a Poisson process. The “under-dispersed” distribution corresponds to the relatively regular randomness. If the index of dispersion is larger than 1, a data-set can correspond to the existence of clusters of occurrences.

In its general setting, the problem of conditional limit behaviour of MBP is studied in the book [19]. It is interesting to compare the infinitesimal geometric branching reproduction with a reproduction defined by the quadratic p.g.f. having the same offspring’s mean and variance. Let the p.g.f. of the offspring’s numbers for the process Y(t) be given by u(s)=(m2−m+1)+s(m−2m2)+m2s2,u′(1)=m,u″(1)=2m2. The infinitesimal generating function is v(s)=K(m−1)(s−1)+Km2(s−1)2,v′(1)=f′(1)=M,v″(1)=f″(1). Then the decomposition v′(1)v(s)=1s−1−ϱϱs−1,ϱ=m2m2+(1−m),1−ϱ=(1−m)m2+(1−m). We denote the corresponding functions expressing the solutions of the Kolmogorov equations by the letters D0 and D , D0(s)=∫0sv′(1)dxv(x)=1−s1−ϱs=(1−s)1−ϱ+ϱ(1−s),D(s)=s1−ϱ+ϱs. Then D0−1(s)=D0(s),D−1(s)=1−D0(s)=(1−ϱ)s1−ϱs=D∗(s). We see that the conditional limit distribution is a shifted geometric distribution with parameter ϱ<1 . The solution of the equation corresponding to (5, 6) is given by r(t,s)=D−1(eMtD(1−s))=eMt(1−s)1+(1−eMt)(1−s)ϱ/(1−ϱ), where ϱ1−ϱ=m21−m=f″(1)−2f′(1). The transient phenomena when the process X(t) is near critical is introduced by Sevastyanov, B. A. in 1959 [19, 18]. It take place when (t→∞,m→1) as follows, R(t,s)=r(t,s)(1+ε(t,s)),ε(t,s)→0,(t→∞,m→1).

Critical geometric branching

Let the reproduction of particles be of mean m=1 , therefore q=1 and mathematical expectation E[X(t)]=1 for any t>0 0$]]>. Hence, the p.g.f. h(s)=12−s,f(s)=K(1−s)22−s,Kf(s)=1(1−s)2+11−s. The backward Kolmogorov equation is the same as (3), but with a new parameter (4), namely f′(1)=0 . The analytical deviation (difference) of critical and subcritical MBP starts with decomposition (7) and hence (8). Now, in the critical case, we have only similarly (analogously), the following indefinite integral primitive, for x≠1 , ∫Kdxf(x)=∫1(1−x)2+11−xdx=log1|1−x|exp11−x. Let us introduce the composite function C(s)=V(G(s)),s≠1 , such that, C(s)=11−sexp11−s,C′(s)C(s)=Kf(s),V(x)=xex,G(s)=11−s. Then obviously, knowing (21), the solution F(t,s) of the backward Kolmogorov equation in the critical case is expressed by the following relation, C(F(t,s))=eKtC(s),|s|<1. The function C(s) (22) has a vertical asymptote at the point s=1 . First of all, we study its behaviour left and right of the vertical s=1 , especially the domain of increasing, and then we look for its inverse. The first derivative C′(s)=2−s(1−s)3exp11−s,s≠1,C(0)=e,C′(0)=2e, shows that C(s) is increasing in the following two intervals −∞<s<1 and 2<s<∞ . Note that the V(G(s)) is negative right of the point s=1 , decreasing in the interval 1<s<2 and has minimum at the point s=2 . Only the interval −∞<s<1 , where C(s) is positive and increasing, is convenient for definition of inverse function, C−1(x)=s,x>0 0$]]>, if and only if x=C(s),s<1 . Obviously, the inverse function C−1(x) has a vertical asymptote at the point x=0 and limx→∞C−1(x)=1 .

An important advantage of this composite function C(s)=V(G(s)) is that all higher order derivatives at zero can be expressed by the Lah numbers [21], as it is shown in the following lemma.

Let the composite function C(s)=V(G(s)),s≠1 satisfies (22). Then the Taylor series expansion of the function C(s) in the neighbourhood of zero is expressed by the Lah numbers as follows, C(s)=e+C0(s),C0(s)= ∑k=1∞ckskk!,cn=e ∑k=1n(k+1)L(n,k), where the Lah numbers L(n,k) are defined by the partial Bell polynomials over the sequence (∙!)=(n!,n=1,2,…),L(n,k)=Bn,k(∙!)= ∑k=1nn−1k−1n!k!.

At first, directly from the definition we have C(0)=V(1)=e . Now, let us consider the derivatives of the function V(x)=xex,V(n)(x)=V(x)+nex. Then, at the point x=1 all derivatives are the multiples of the number e, as follows, V(1)=e,V′(1)=2e,…,V(n)(1)=(n+1)e,n=2,3,…, and at the point x=−1 as follows V(−1)=−e−1,V′(−1)=0,…,V(n)(−1)=(n−1)e−1>0,n=2,3,…. 0,\hspace{1em}n=2,3,\dots .\]]]> We remark that, the first derivative V′(x)<0 in the interval (−∞,−1) and the function V(x) is decreasing in this interval. But V(n)(−n)=0 and the function V(x) changes the convexity at these points. Moreover, limx→∞V(n)(x)=+∞,limx→−∞V(n)(x)=0. The linear-fractional function G(s) has the following derivatives for s≠1 , G(s)=11−s,G(n)(s)=n!(1−s)n+1,G(n)(0)=n!. We have, left of the point s=1−,s<1 , lims→1−G(n)(s)=+∞,n=0,1,2,…, and right of the point s=1+,s>1 1$]]>, lims→1+G(2n)(s)=−∞,lims→1+G(2n+1)(s)=+∞,n=0,1,2,…. The n−th derivative of the composite function C(s)=V(G(s)) at any point s≠1 is given by the Faa Di Bruno formula, C(n)(s)= ∑k=1nV(k)(G(s))Bn,k(G∙),(G∙)=(G(n)(s),n=1,2,…). In particular, at the point s=0 the corresponding derivatives are G(0)=1,V(k)(1)=e(k+1),G(j)(0)=j!,Bn,k(∙!)=L(n,k). Namely, C(n)(0)= ∑k=1nV(k)(1)Bn,k(g∙),(g∙)=(gn=G(n)(0),n=1,2,…). The Taylor series expansion confirms the representation (23) with cn=C(n)(0) .  □

A direct outcome of from the lemma is the explanation of how the inverse function C−1(x),x>0 0$]]>, approaches its horizontal asymptote. The conclusion follows from (24), implying the asymptotic behaviour of derivatives in the neighbourhood of the vertical asymptote, left and right, lims→1−C(n)(s)=+∞,lims→1+C(n)(s)=0,C(n)(1+1/k)=0,k=2,3,…. The purpose of the lemma is that if the coefficients ck in the Taylor series expansion (23) of the function C(s) are known, we can apply the Lagrange inversion method in its general setting to obtain the representation C0(s)= ∑k=1∞ckskk!,C0−1(x)= ∑k=1∞dkxkk!,C−1(x)=C0−1(x−e),x>0, 0,\]]]> where d1c1=1 and the coefficients dn , for n=2,3,… , are given by (2e)ndn= ∑k=1n−1(−1)k[n]k↑Bn−1,k(q∙),(q∙)=(qn=cn+1/(n+1)c1). However, this approach expects applications to perform large iterative computations in aim to obtain highly precise results. For this reason, we looked for an inversion based on the composite function V(G(s)) by introducing the Lambert-W function (considered as inverse of V(x)), to obtain another more effective solution.

Lambert-W function and Lagrange inversion

In general, the Lambert-W function is a complex-valued function with branching point z=−e−1 . We consider only its principal branch as the real-valued inversion function of V(s)=ses,s≥−1,W(x)=V−1(x),x≥−e−1. The inversion of the composite function V(G(s)) is defined as follows, C(s)=V(G(s)),s<1,C−1(x)=G−1(V−1(x))=G−1(W(x)),x>0. 0.\]]]> We consider only interval x>0 0$]]> because the function G−1 has a vertical asymptote at zero, W(0)=0 . Then, the solution F(t,s) of the backward Kolmogorov equation is given by, 11−F(t,s)=WeKt(1−s)exp11−s,|s|<1. This representation of the p.g.f. F(t,s) is very convenient to realize the polynomial rate of increasing with time parameter t>0 0$]]>, based on the logarithmic rate of increasing of the Lambert-W function [4], formula (28) . In particular, F(t,0)=1−1W(eKt+1)<F(t,s)<1,0<s<1,t>0. 0.\]]]> The family of functions {F(t,s),0<s<1,t>0,} 0,\}$]]> converges to the constant function y(s)=1 uniformly by |s|<1 and increasing by t>0 0$]]> [19].

The expression for the derivatives of Lambert-W function on the domain of differentiability, −e−1<x<∞ is given in [4, 15] as dndxnW(x)=(−1)n−1An(W(x))exp(−nW(x))(1+W(x))2n−1,n=1,2,…, where An(x) is a polynomial of degree (n−1) specified by the recursion, A1(x)=1 , and for n≥1 , An+1(x)=(3n−1+nx)An(x)−(1+x)An′(x). For x≠0 it is convenient to write the derivatives as follows, W(n)(x)=(−1)n−1An(W(x))(1+W(x))2n−1W(x)xn,n=1,2,…. Then the consecutive derivatives of the Lambert-W function at the point x=e are easily calculated as, W(n)(e)=(−1)n−1An(1)22n−11en=(−1)n−1An(1)2n−1(2e)n,n=1,2,…. The relation between inversion of the composite function V(G(s)) and Lambert-W function is defined and proved in the following theorem.

Let the composite function C(s)=V(G(s)) satisfies C(s)=11−sexp11−s=e+C0(s),C0(s)= ∑k=1∞ckskk!,s<1, then C0−1(x)= ∑k=1∞dkxkk!,C−1(x)=C0−1(x−e),x>0, 0,\]]]> and dn= ∑k=1n(−1)k−1k!Bn,k(w∙), where the sequence (w∙)=(wn=W(n)(e),n=1,2,…) is defined by the derivatives of the Lambert-W function at the point x=e .

In this case, knowing the representation (26), the coefficients dn can be calculated directly by derivatives of the representation C0−1(x)=C−1(x+e)=1−1W(x+e)=G−1(W(x+e)), at the point x>0 0$]]>. The function G−1(x) has the following derivatives, G−1(x)=1−1x,1−1x(k)=(−1)k−1k!xk+1,(G−1)(k)(1)=(−1)k−1k!. At the point x=0 the function W(x+e) take values W(e)=1 . Applying the Faa Di Bruno formula on the composite function (30) we obtain, dn= ∑k=1n(−1)k−1k!Bn,k(w∙),(w∙)=(wn=W(n)(e),n=1,2,…), We remark the equivalence of (25) obtained by the Lagrange inversion method and (29) expressing dn by the values of the Lambert-W function.  □

The assumption for computational simplification is justified during the calculation of the very important for every branching process extinction probability. It can be computed from (27) applying any software packages giving the values of the function W(eKt+1) . In our case, we used the function lambertW from package VGAM (see [22]) in R environment for statistical computing, where the Lambert-W functionality properties are implemented according [5]. The computed results are shown in Figure 3.

For comparison, the series expansion of F(t,0) is expressed by the Stirling numbers of second kind S(j,n) and coefficients dn as follows, F(t,0)= ∑j=1∞rj(Kt)jj!,rj= ∑n=1jS(j,n)endn. It is easy to notice that the value F(t,0)=P(X(t)=0) , representing the extinction probability to the positive time t>0 0$]]>, is F(t,0)=C−1(eKt+1)=C0−1(eKt+1−e),C(0)=e. Thus, after applying the Theorem (4) we obtain C0−1(eKt+1−e)=C0−1(e(eKt−1))= ∑n=1∞endn(eKt−1)nn!. The properties of exponential generating function and definition of the partial exponential Bell polynomials, as it is shown in [21], allow writing (eKt−1)nn!=1n! ∑j=1∞(Kt)jj!n= ∑j=n∞Bj,n(1∙)(Kt)jj!, where Bj,n(1∙)=S(j,n),1∙=(1,12,13,….) are the Stirling numbers of the second kind. In this way the extinction probability is obtained as expansion by double sum F(t,0)= ∑n=1∞endn ∑j=n∞S(j,n)(Kt)jj!. The change of summation order n≥1,j≥n into j≥1,1≤n≤j leads to (31).

For the practical implementation of this result, it is easier to calculate several terms of the coefficients rj with any computational tool. Then, knowing these coefficients rn the approximation of the extinction probability is expanded as F(t,0)=12.Kt−323.(Kt)22!+1125.(Kt)33!−4527.(Kt)44!+19329.(Kt)55!+⋯. On the other side, the series expansion of the Lambert-W function [4], see formula (28) with the radius of convergence 4+π2 , is W(ez)=1+(z−1)2+(z−1)216−(z−1)3192−(z−1)43072+13(z−1)561440.... Applied on (27), it gives the following representation for the extinction probability 11−F(t,0)=1+(Kt)2+(Kt)216−(Kt)3192−(Kt)43072+13(Kt)561440........ The division of polynomials in the increasing order confirms the previous approximation obtained by the coefficients rn from (31) and agreement of two series expansions.

The series expansion obtained by the Lagrange inversion method converges very slowly with growing t>0 0$]]> requiring expansion to larger order number n. This issue implies the already known computational problem in similar applications, in particular, for calculation of eigenstates of hydrogen molecular ion H+ [13, 17]. To solve this computational issue, the representation of F(t,s) and its derivatives at the point s=0 is redefined in terms of the values W(eKt+1) , calculated with the corresponding software packages.

The probability mass function of X(t),t>0 0$]]>, is defined by the derivatives of the p.g.f. F(t,s) at zero as follows P(X(t)=n)=1n!Fs(n)(t,0),n=0,1,2,…. Obviously, for n=0 we have the extinction probability. Applying (28) on the derivatives of the Lambert-W function, we give the representation of p.m.f. in terms of the values W(eKt+1) in the following theorem.

Let X(t),t>0 0$]]>, be a critical time-homogeneous MBP with branching mechanism given by (20). The consecutive derivatives of the p.g.f. F(t,s) at zero are Fs(n)(t,0)= ∑k=1n(eKt)kBn,k(c∙) ∑j=1k(−1)j−1j!Bk,j(W∙)(W(eKt+1))j+1. where the sequence (W∙) is given by the derivatives of the Lambert-W function at the point eKt+1 .

We start with the representation F(t,s)=C−1(eKtC(s)). Then the n′th derivative by s at the point s=0 is given by Fs(n)(t,0)= ∑k=1n(C−1)(k)(eKtC(0))Bn,k(eKtc∙). To obtain the derivatives of the inverse function C−1 at the point C(0)=e we apply the Faa Di Bruno formula on the representation C−1(x)=G−1(W(x)). Knowing the derivatives of the linear-fractional function at the point x=eKt+1 , G−1(x)=1−1x,(G−1)(k)(x)=(−1)k−1k!xk+1, we obtain for the derivatives of the composite function the following representation (C−1)(k)(eKtC(0))=(C−1)(k)(eKt+1)= ∑j=1k(−1)j−1j!Bk,j(W∙)(W(eKt+1))j+1, where the sequence (W∙) is defined by derivatives of the Lambert-W function at the point eKt+1 as follows (W∙)=(W(k)(eKt+1),k=1,2,…).  □

As a demonstrative example, several derivatives presenting probabilities P(X(t)=n) are computed and their values are shown in Figure 3. Fs′(t,0)=2W(1+W),Fs″(t,0)=(W−1)(3W+1)W(1+W)3,F(3)(t,0)=2(W−1){4W3+2W2+5W+1}W(1+W)5,F(4)(t,0)=(W−1)(30W5+28W4+91W3+3W2+35W+5)W(1+W)7.

Probabilities for P(X(t)=n) . When process is in extinction, i.e. n=0 , the graphics is shown in black colour. Otherwise, when n=1 is shown in brown, n=2 in red, n=3 in green and n=4 in blue

 Conclusions

The infinitesimal geometric branching reproduction describes the models wherein at any time t>0 0$]]> there is a family of particles from all generations (with positive probability). Similar physical problems are very difficult and time consuming to solve analytically, leading to the preference of numerical methods. However, with probabilistic methods, we found explicit solutions for subcritical and critical cases.

 Acknowledgments

The authors are very thankful to the anonymous referees for the valuable comments on the paper.

 References Athreya, K.B., Ney, P.E.: Branching Processes. Springer, New York (1972). MR0373040 Brito, P., Fabiao, F., Staubyn, A.: Euler, Lambert, and the Lambert W-function today. Math. Sci. 33, 127–133 (2008). MR2474739 Consul, P.C., Famoye, F.: Lagrangian Probability Distributions. Birkhäuser, Boston (2006). MR2209108 Corless, R.M., Jeffrey, D.J., Knuth, D.E.: A sequence of series for the Lambert-W function. In: Proceedings of the 1997 International Symposium on Symbolic and Algebraic Computation. ISSAC ’97, pp. 197–204. ACM, New York (1997). MR1809988. https://doi.org/10.1145/258726.258783 Corless, R.M., Gonnet, G.H., Hare, D., Jeffrey, D.J., Knuth, D.E.: On the Lambert W function. Adv. Comput. Math. 5, 329–359 (1996). MR1414285. https://doi.org/10.1007/BF02124750 Fukushima, T.: Precise and fast computation of Lambert W-functions without transcendental function evaluations. J. Comput. Appl. Math. 244, 77–89 (2013). MR3005760. https://doi.org/10.1016/j.cam.2012.11.021 Harris, T.E.: The Theory of Branching Processes. Springer, Berlin (1963). MR0163361 Johnson, N.L., Kemp, A.W., Kotz, S.: Univariate Discrete Distributions. John Wiley & Sons (2005). MR2163227. https://doi.org/10.1002/0471715816 Kovalenko, I.N., Kuznetchov, N.J., Shurenkov, V.M.: Random Processes. Naukova dumka, Kiev (1983). MR0751235 Luchko, Y., Trujillo, J., Velasco, M.P.: The Wright function and its numerical evaluation. Int. J. Pure Appl. Math. 64 (4), 567–575 (2010). MR2796971 Mainardi, F., Mura, A., Pagnini, G.: The M-Wright Function in Time-Fractional Diffusion Processes: A tutorial survey. Int. J. Differ. Equ. 2010, 104505 (2010). MR2592742. https://doi.org/10.1155/2010/104505 Merlini, D., Sprugnoli, R., Verri, M.C.: Lagrange inversion: When and How. Acta Appl. Math. 94, 233–249 (2006). MR2290868. https://doi.org/10.1007/s10440-006-9077-7 Mező, I., Keady, G.: Some physical applications of generalized Lambert functions. Eur. J. Phys. 37 (6), 065802 (2016). https://doi.org/10.1088/0143-0807/37/6/065802 Pakes, A.G.: Characterization of discrete laws via mixed sums and Markov branching processes. Stoch. Process. Appl. 55 (2), 285–300 (1995). MR1313024. https://doi.org/10.1016/0304-4149(94)00049-Y Pakes, A.G.: Lambert’s W, infinite divisibility and Poisson mixtures. J. Math. Anal. Appl. 378, 480–492 (2011). MR2773259. https://doi.org/10.1016/j.jmaa.2011.01.066 Paris, R.B., Vinogradov, V.: Asymptotic and structural properties of special cases of the Wright function arising in probability theory. Lith. Math. J. 56, 377–409 (2016). MR3530225. https://doi.org/10.1007/s10986-016-9324-1 Scott, T.C., Fee, G., Grotendorst, J.: Asymptotic series of generalized Lambert-W function. ACM Commun. Comput. Algebra 47 (3-4), 75–83 (2014). MR3167418. https://doi.org/10.1145/2576802.2576804 Sevastyanov, B.A.: Transient phenomena in branching stochastic processes. Theory Probab. Appl. 4 (2), 113–128 (1959). MR0121854. https://doi.org/10.1137/1104011 Sevastyanov, B.A.: Branching Processes. Nauka, Moscow (1971). MR0345229 Snyder, D.L., Miller, M.I.: Random Point Processes in Time and Space, Second Edition. Springer, New York (1991). https://doi.org/10.1007/978-1-4612-3166-0 Wang, W., Wang, T.: General identities on Bell polynomials. Comput. Math. Appl. 58, 104–118 (2009). MR2535972. https://doi.org/10.1016/j.camwa.2009.03.093 Yee, T.W.: The vgam package for categorical data analysis. J. Stat. Softw. 32 (10), 1–34 (2010). https://doi.org/10.18637/jss.v032.i10