Geometric branching reproduction Markov processes

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.


Introduction
Any time-homogeneous Markov branching process (MBP) X(t), t > 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, 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, is expressed by the values of the Lambert-W function at the point x = e Kt+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) = xe x , 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, 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. The reproduction is defined by the integer-valued random variable η representing the offspring numbers as follows, The parameter m > 0 defines the mean of the offspring numbers. In this notation the p.g.f. of the reproduction law is, The infinitesimal generating function f (s) present in the Kolmogorov equations is defined by the constant K and p.g.f. h(s) as and has the following derivatives, for k = 1, 2, 3, . . ., The p.g.f. of the number of particles alive at the time t > 0 yields the backward Kolmogorov equation The equation h(s) = s has two solutions, s 1 = 1/m and s 2 = 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 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, In the subcritical case, we study in parallel the p.g.f. F (t, s) of the branching process X(t), t > 0, defined by (2) and equation (3), and the function satisfying the following non-linear equation: In order to solve these equations (3) and (6), when f (1) < 0, we define respectively the function and the function see [14,1]. We remark that in the subcritical case The solution of the indefinite integral is Knowing the primitive (8) of the equation (3) we formulate the following lemma without proof.

the composite inverse functions, and
The mathematical expectation of the number of particles alive at the time t > 0 in the subcritical case has the following exponential decreasing behaviour, The extinction probability to the positive time P (X(t) = 0) = F (t, 0), t > 0, is expressed by the inversion function Respectively, the survival probability is given as follows, We remark that at the point s = 1/2 the functions A 0 (1/2) = A(1/2) and the respective vertical asymptotes are left and right symmetric to s = 1/2. The function A 0 (s) is decreasing, concave and all its derivatives A In details, denoting the falling and rising factorials as follows: and we write the derivatives in the following form, In particular, for k = 1, 2, . . ., at the points s = 0 and s = 1 we have and Obviously, the following relation follows from (11) and (12) 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/m m+1 , 0 < m < 1.

Lagrange inversion and special functions
Now, we look for the properties of the composite inverse of the functions A 0 (s) and A(s). The function A(s) and its inverse A −1 (s) are some particular cases of the Wright function [11,16] and the Gauss hypergeometric function [8,3] Obviously, knowing the notation θ = m 1−m > 0, 0 < m < 1, we see that the Using the previous definitions, it will be proved in the following theorem that for |x| < 1/m m+1 and 0 < m < 1 the inverse function Theorem 1. Let us consider the function A 0 (s) given by

Then for its composite inverse function is valid
Proof. First of all, we confirm the relation The function A(s) is represented as an exponential generating function in the neighbourhood of zero |s| To apply the Lagrange inversion method to it, we use the following representation, and obviously g(1) = 1, see [3,8]. Denote the inverse function in its series expansion as Then the coefficients b k are given by the derivatives of the function (g(s)) k at the point s = 0 as follows, The Taylor series expansion of the function is given by the binomial coefficients as follows, The derivative of the order j to this function at the point s = 0 is It is enough to take j = k − 1 in order to obtain the coefficient b k in (15). Then Applying the definition of decreasing factorials by Gamma function (9, 10), The change of variable k − 1 = j leads to the representation (14). We remark that the coefficients b k are either positive, zero or negative, due to the decreasing factorials. It is confirmed, when the reflection formulas for Gamma function are applied on the representation Obviously, for k(1 − m) = j, j = 2, 3, . . ., the value b k = 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 is calculated with a root test based on and Stiling's formula for the Gamma function. When k → ∞ we have We extract the multiple Then, obviously, And finally, 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, 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 , Respectively, their radiuses of convergence are as follows In the aim to describe the behaviour of the process X(t), t > 0, we consider the factorial moments E[X(t)] n↓ = F (n) s (t, 1). The derivatives of the p.g.f. F (t, s) are expressed by the partial Bell polynomials defined as see [21]. The sum is over all partitions of n into k parts, that is over all nonnegative integer solutions (k 1 , k 2 , . . . , k n ) of the equations: k 1 + 2k 2 + · · · + nk n = n, k 1 + k 2 + · · · + k n = k.
The following theorem is based on the inequality (17) showing that the convergence interval of the series ∞ k=1 b k s k /k! contains the point s = 1.
Proof. As we have seen previously, the solution of the equation (6) is 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 Then, exchanging the order of summation we obtain from (18), 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 Finally, the important process behaviour like conditional limit probability is obtained in the following theorem.

It means that all factorial moments F
The factorial moments E[ξ ] k↓ of the limit random variable is Proof. The p.g.f. of the conditional probability is written as follows Knowing the differentiability of the function A −1 in the neighborhood of zero, we apply the theorem of l'Hospital to find the limit To obtain the p.m.f. we consider the Taylor series expansion of the p.g.f. F * (s) in the neighbourhood of zero, The change of variable j + 1 = k in the second sum leads to, We remark that The relation between derivatives of the functions A and A 0 at the points zero and s = 1 is obvious. In this way, we recognise the factorial moments f k as given previously in (11), and also their relation with the p.m.f. p k 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 p k 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, . . .., This ratio (19) shows that the Panjer recursion (even generalized) does not take place. The inequality ϒ k (m) < 1 is equivalent to the following where (m + 1 − 2m 2 ) > 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 The histogram of the p.m.f. can be generated by the consecutive multiplication with ϒ k (m) < 1, where the maximum is given by p 1 = 1 − m 2 . Some results for ϒ k (m) and computed from them probabilities p k 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, The central moments, respectively, are as follows The main parameter for applications is the variance-to-mean ratio (VMR) or index of dispersion equal respectively to 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 with a precision higher than 10 −5 . The conditional limit probability is over-dispersed when the branching process X(t), t > 0, is near critical, m * < m < 1, and underdispersed when 0 < m < m * < 1. The values of VMR are shown in Figure 2. In the class of Panjer probability distributions, the VMR = 1 for the Poisson distribution, respectively, VMR > 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 "underdispersed" 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.

Remark 1.
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 (1).

The infinitesimal generating function is
We denote the corresponding functions expressing the solutions of the Kolmogorov equations by the letters D 0 and D, 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 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,

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. Hence, the p.g.f.
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, Let us introduce the composite function C(s) = V (G(s)), s = 1, such that, 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)) = e Kt C(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 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, 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 lim x→∞ 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. The linear-fractional function G(s) has the following derivatives for s = 1,

Lemma 2. 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,
We have, left of the point s = 1 − , s < 1, In particular, at the point s = 0 the corresponding derivatives are Namely, The Taylor series expansion confirms the representation (23) with c n = C (n) (0).
A direct outcome of from the lemma is the explanation of how the inverse function C −1 (x), x > 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, The purpose of the lemma is that if the coefficients c k 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 where d 1 c 1 = 1 and the coefficients d n , for n = 2, 3, . . ., are given by 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 The inversion of the composite function V (G(s)) is defined as follows, We consider only interval x > 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, 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, based on the logarithmic rate of increasing of the Lambert-W function [4], formula (28). In particular, The family of functions {F (t, s), 0 < s < 1, t > 0, } converges to the constant function y(s) = 1 uniformly by |s| < 1 and increasing by t > 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 where A n (x) is a polynomial of degree (n − 1) specified by the recursion, A 1 (x) = 1, and for n ≥ 1, For x = 0 it is convenient to write the derivatives as follows, The relation between inversion of the composite function V (G(s)) and Lambert-W function is defined and proved in the following theorem.

Theorem 4. Let the composite function C(s) = V (G(s)) satisfies
where the sequence (w • ) = (w n = W (n) (e), n = 1, 2, . . .) is defined by the derivatives of the Lambert-W function at the point x = e.
Proof. In this case, knowing the representation (26), the coefficients d n can be calculated directly by derivatives of the representation at the point x > 0. The function G −1 (x) has the following derivatives, 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 (e Kt+1 ). In our case, we used the function lambertW from package V GAM (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 d n as follows, S(j, n)e n d n .
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, is Thus, after applying the Theorem (4) we obtain e n d n (e Kt − 1) n n! .
The properties of exponential generating function and definition of the partial exponential Bell polynomials, as it is shown in [21], allow writing where B j,n (1 • ) = S(j, n), 1 • = (1, 1 2 , 1 3 , . . . .) are the Stirling numbers of the second kind. In this way the extinction probability is obtained as expansion by double sum 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 r j with any computational tool. Then, knowing these coefficients r n the approximation of the extinction probability is expanded as (Kt) 4 4! + 193 2 9 .
(Kt) 5 On the other side, the series expansion of the Lambert-W function [4], see formula (28) with the radius of convergence Applied on (27), it gives the following representation for the extinction probability The division of polynomials in the increasing order confirms the previous approximation obtained by the coefficients r n from (31) and agreement of two series expansions. The series expansion obtained by the Lagrange inversion method converges very slowly with growing t > 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 (e Kt+1 ), calculated with the corresponding software packages.
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 (e Kt+1 ) in the following theorem.
Theorem 5. Let X(t), t > 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 where the sequence (W • ) is given by the derivatives of the Lambert-W function at the point e Kt+1 .
Proof. We start with the representation F (t, s) = C −1 (e Kt C(s)).
Then the n th derivative by s at the point s = 0 is given by 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 Knowing the derivatives of the linear-fractional function at the point x = e Kt+1 , we obtain for the derivatives of the composite function the following representation   , F s (t, 0) = (W − 1)(3W + 1)

Conclusions
The infinitesimal geometric branching reproduction describes the models wherein at any time t > 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.