VMSTA Modern Stochastics: Theory and Applications 2351-60542351-60462351-6046 VTeXMokslininkų g. 2A, 08412 Vilnius, Lithuania VMSTA201 10.15559/22-VMSTA201 Research Article Factorial moments of the critical Markov branching process with geometric reproduction of particles https://orcid.org/0000-0001-9322-262X 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. 2022 722022922292443062021261120212412022 © 2022 The Author(s). Published by VTeX2022 Open access article under the CC BY license.

The factorial moments of any Markov branching process describe the behaviour of its probability generating function F ( t , s ) in the neighbourhood of the point s = 1. They are applied to solve the forward Kolmogorov equation for the critical Markov branching process with geometric reproduction of particles. The solution includes quickly convergent recurrent iterations of polynomials. The obtained results on factorial moments enable computation of statistical measures as shape and skewness. They are also applicable to the comparison between critical geometric branching and linear birth-death processes.

Factorial moments geometric reproduction branching process harmonic numbers Lambert-W function Stirling numbers 60J80 60K05 Bulgarian National Science FundKP-06-H22/3This research has been partially supported by the Bulgarian National Science Fund, grant No. KP-06-H22/3.
Introduction

A branching process model driven by geometric reproduction of particles was introduced in . By design, it is a time-homogeneous Markov branching process X ( t ) , t > 00$]]>, with probability generating function F ( t , s ) , | s | < 1. The geometric branching mechanism is defined by the integer-valued random variable η representing the particles offspring number with probability mass function (p.m.f.) P ( η = n ) = ( 1 ϱ ) ϱ n , n = 0 , 1 , 2 , , 0 < ϱ < 1 . The lifetime of particles is an exponentially distributed random variable with a constant parameter K > 00$]]>. The initial condition is X ( 0 ) = 1. The time interval to the first splitting is exponentially distributed with density K e K t . We remark that at time t > 00$]]>, when the number of particles in the system is X ( t ) = n, the time interval to the next splitting is exponentially distributed with density n K e n K t as the minimum of n independent exponentially distributed random variables. The probability generating function (p.g.f.) h ( s ) : = E [ s η ], | s | 1, and the mathematical expectation m = E [ η ] are h ( s ) = E [ s η ] = 1 ϱ 1 ϱ s , | s | 1 , m = h ( 1 ) = ϱ 1 ϱ > 0 . 0.\]]]> The reproduction is classified as critical when the mean of progeny m = 1, i.e. the parameter ϱ = 1 / 2. Then the p.g.f. h ( s ) is reduced to h ( s ) = 1 2 s , h ( s ) = 1 ( 2 s ) 2 , h ( s ) = 2 ( 2 s ) 3 . The p.g.f. F ( t , s ) , | s | < 1, is determined in  as the unique solution to the Kolmogorov equations  after the Lagrange inversion method following the classical theory . The factorial moments in the subcritical case are computed in . They are used to obtain the probability mass function, central moments and variance-to-mean ratio (VMR) for the limiting random variable. In critical case, the p.g.f. F ( t , s ) is defined by the composition of the Lambert-W and linear-fractional functions. The probability mass function of X ( t ) , t > 00$]]>, is expressed by the values of the Lambert-W function at the point x = e K t + 1 . The agreement between Lagrange inversion method and series expansion of the Lambert-W function is confirmed by the approximation of extinction probability. The detailed explanation can be found in .

The graph of the Lambert-W function W ( x ) for 1 / e x 5. The real principal branch of solution W 0 is shown in red. The second branch of real solution W 1 in the interval 1 / e x < 0 is drawn in blue. The used software implementation in R is available in 

However, the use of the Lambert-W function on its principal branch W 0 ( x ) of x 1 / e is limited to obtaining values of F ( t , s ) only for s < 1 . Thus, major statistical moments as variance and skewness remain undefined. The graph of the Lambert-W function is shown in Figure 1. Very precise approximations for F ( t , s ) , s 1, can be computed by any dedicated software packages. But, the real positive values of the p.g.f. F ( t , s ) for s > 11$]]> can be derived only from the W 1 ( x ) branch of 1 / e x < 0, as it is shown in Figure 2. To find analytical solution, the series expansion of the p.g.f. F ( t , s ) in the neighbourhood of s = 1 is obtained. The results follow from the solution to the forward Kolmogorov equation and rely on the factorial moments of the critical infinitesimal geometric branching reproduction. The graphs of F ( t , s ) computed with the Lambert-W function at the point x = e K t C ( s ) (see (5), (6)) for different values of K t (computation details can be found in ). The values of F ( t , s ) for s > 11$]]> are computed with the W 1 branch of the Lambert-W function

In the next Section 2, the description of the proposed model is given through the solution to the backward Kolmogorov equation. We outline there the problem of continuity of the p.g.f. F ( t , s ) at the point s = 1 for s 1. Section 3 is devoted to the factorial moments of the number of particles alive to time t > 00$]]>. The recurrent relation is established. The main result, concerning the representation of the p.g.f. F ( t , s ) as a convergent series in the domain | s 1 | < 1 / K t, is proved in Section 4. Finally, the statistical inference is introduced in Section 5. The series expansion of many special functions is based on the combinatorics relations, especially for the Lambert-W function. In Appendix we state some useful formulas for quick orientation on the outcome of the summation by diagonals introduced in Section 4. Backward Kolmogorov equation for critical geometric branching mechanism Suppose that the lifetime of particles is exponentially distributed with constant rate K > 00$]]>. Then with positive probability, there is a family of particles from all generations at time t > 00$]]>. The number of particles alive at time t > 00$]]> is described by a Markov branching process (MBP) X ( t ) , t 0, X ( 0 ) = 1, . Only the critical case, defined by (1) with parameter ϱ = 1 / 2, is considered. Then the reproduction of particles is of mean m = 1, therefore the ultimate extinction probability q = lim t F ( t , 0 ) = 1 and the mathematical expectation E [ X ( t ) ] = 1 for t > 00$]]>, see [10, Theorem 4, page 53]. For the reproduction p.g.f. introduced in (1) the infinitesimal generating function is defined as f ( s ) : = K ( h ( s ) s ) = K ( 1 s ) 2 2 s , K f ( s ) = 1 ( 1 s ) 2 + 1 1 s . The probability generating function of the process X ( t ) , t > 00$]]>, defined by F ( t , s ) = k = 0 s k P ( X ( t ) = k | X ( 0 ) = 1 ) satisfies the backward Kolmogorov equation t ( F ( t , s ) ) = f ( F ( t , s ) ) , F ( 0 , s ) = s . Equation (3) is nonlinear with separable differentials due to the property of time-homogeneity. It can be written as d F f ( F ) = d t , f ( s ) = K ( 1 s ) 2 2 s , f ( 0 ) = K 2 . It means that the indefinite integral obtains the form K d x f ( x ) = 1 ( 1 x ) 2 + 1 1 x d x = log 1 | 1 x | exp 1 1 x . Now, for the implicit solution to (3) in the critical case the following equation is valid: 1 1 F ( t , s ) log ( 1 F ( t , s ) ) = 1 1 s log ( 1 s ) + K t , | s | < 1 . Once we have (4), we can find the explicit solution to (3). It is defined by the composite function C ( s ) = V ( G ( s ) ), for s 1, introduced in the form C ( s ) = 1 1 s exp 1 1 s , C ( s ) C ( s ) = K f ( s ) , V ( x ) = x e x , G ( s ) = 1 1 s . Then, from Equation (4), the solution F ( t , s ) (see (2)) to the backward Kolmogorov equation in the critical case is C ( F ( t , s ) ) = e K t C ( s ), | s | < 1. The function C ( s ) from Equation (5) has a vertical asymptote at the point s = 1 (see Figure 3). Its first derivative C ( s ) = 2 s ( 1 s ) 3 exp 1 1 s , s 1 , C ( 0 ) = e , C ( 0 ) = 2 e , shows that C ( s ) is increasing in the intervals < s < 1 and 2 < s < . It is negative right to the point s = 1, decreases in the interval 1 < s < 2, and has a minimum at the point s = 2. Only the interval < s < 1, where C ( s ) is positive and increasing, is convenient for definition of the inverse function, C 1 ( x ) = s, if and only if, x = C ( s ).

The graphs of e K t C ( s ) = e K t V ( G ( s ) ), s 1, as a function of s for different K t. The horizontal red dashed line marks the min C ( s ) = C ( 2 ) = 1 / e. The result lines crossing dashed line are inflated by e K t , but they preserve the abscissa of minimum min ( e K t C ( s ) ) = e K t C ( 2 ) = e K t / e < 1 / e

The composite inverse function of V ( G ( s ) ) is obtained in  by introduction of the Lambert-W function. As the Lambert-W function is multivalued with a branching point at z = e 1 , the solution is limited only to the principal branch W 0 ( x ), shown in Figure 1, where V ( s ) = s e s , s 1 , W ( x ) = V 1 ( x ) , x e 1 . The inversion of the composite function V ( G ( s ) ) for s < 1 is defined as C ( s ) = V ( G ( s ) ) , C 1 ( x ) = G 1 ( V 1 ( x ) ) = G 1 ( W ( x ) ) , x > 0 . 0.\]]]> Then, the solution F ( t , s ) to the backward Kolmogorov equation is given by 1 1 F ( t , s ) = W e K t ( 1 s ) exp 1 1 s , | s | < 1 . This form of the p.g.f. F ( t , s ) is very convenient to realise the polynomial rate of increase with time parameter t > 00$]]>. It is a result derived from the rising logarithmic rate of the Lambert-W function [6, formula (28)]. In particular, the extinction probability P ( X ( t ) = 0 ) by time t > 00$]]> is directly computed from (6) as F ( t , 0 ) = 1 1 W ( e K t + 1 ) < F ( t , s ) < 1 , 0 < s < 1 , t > 0 . 0.\]]]> The convergence of the family of functions F ( t , s ), t > 00$]]>, | s | < 1, to the constant function equal to one on the interval 0 s 1 is monotone by t > 00$]]> and uniform by | s | 1 ε < 1, ε > 00$]]>, . Moreover, the following limit exists: lim s 1 F ( t , s ) = 1. The conditional mathematical expectation, knowing the nonextinction, is given by E [ X ( t ) | X ( t ) > 0 ] = 1 1 F ( t , 0 ) = W ( e K t + 1 ) . 0]=\frac{1}{1-F(t,0)}=W({e^{Kt+1}}).\]]]> Forward Kolmogorov equation and factorial moments The inconsistent part of the solution (5), (6) for the p.g.f. F ( t , s ) is when s > 11$]]>. The asymptotic behaviour of derivatives of the function C ( s ) in the neighbourhood of the vertical asymptote, left and right, are lim s 1 C ( n ) ( s ) = + , lim s 1 + C ( n ) ( s ) = 0 , C ( n ) ( 1 + 1 / k ) = 0 , k = 2 , 3 , (see again Figure 3). This is a serious obstacle for studying the continuity of F ( t , s ) at the point s = 1 (left and right).

A possible resolution to this problem is the method of factorial moments. It begins with the definition of the forward Kolmogorov equation, F ( t , s ) t = f ( s ) F ( t , s ) s , F ( 0 , s ) = s , f ( s ) = K ( s 1 ) 2 2 s , with partial derivatives of the p.g.f. F ( t , s ) denoted as follows: n F ( t , s ) s n = F s ( n ) ( t , s ) , | s | < 1 , F s ( n ) ( t , 1 ) = lim s 1 F s ( n ) ( t , s ) , s < 1 . Next, the required n-th factorial moments can be expressed by E [ X ( t ) ] n : = F s ( n ) ( t , 1 ). For them, it is valid that if f ( n ) ( 1 ) < then F s ( n ) ( t , 1 ) < , directly following well-known results from [10, Theorem 2, page 33]. Now, the factorial moments of a critical geometric branching process are well defined and can be obtained by recurrence, as it is proved in the following lemma.

Let X ( t ) be a critical MBP with geometric branching mechanism. Then all factorial moments exist, E [ X ( t ) ] n : = E [ X ( t ) ( X ( t ) 1 ) ( X ( t ) n + 1 ) ] = F s ( n ) ( t , 1 ) and they satisfy the recurrence relation F s ( n + 1 ) ( t , 1 ) = ( n + 1 ) F s ( n ) ( t , 1 ) + K n ( n + 1 ) x = 0 t F s ( n ) ( x , 1 ) d x , n = 1 , 2 , .

We consider the consecutive derivatives of the forward Kolmogorov equation (7) for any critical MBP, knowing that F ( t , 1 ) = F s ( t , 1 ) = 1 and F s ( t , 1 ) = f ( 1 ) t, see : F s ( n ) ( t , s ) t = j = 0 n n ! j ! ( n j ) ! f ( j ) ( s ) F s ( n + 1 j ) ( t , s ) , | s | < 1 , n = 1 , 2 , . Let s = 1. Then for the Markov process with geometric branching mechanism it is valid that f ( 1 ) = 0, f ( 1 ) = 0, f ( j ) ( 1 ) = j ! K for j = 2 , 3 , 4 , . Thus, F s ( n ) ( t , 1 ) t = j = 2 n n ! K j ! j ! ( n j ) ! F s ( n + 1 j ) ( t , 1 ) , n = 2 , 3 , . It is equivalent to F s ( n ) ( t , 1 ) t = K n ! F s 0 ! + F s 1 ! + + F s ( n 2 ) ( n 3 ) ! + F s ( n 1 ) ( n 2 ) ! , n = 2 , 3 , . The next derivative for n = 1 , 2 , becomes F s ( n + 1 ) ( t , 1 ) t = K ( n + 1 ) n ! F s 0 ! + F s 1 ! + + F s ( n 1 ) ( n 2 ) ! + K ( n + 1 ) n ! F s ( n ) ( n 1 ) ! . It leads us to F s ( n + 1 ) ( t , 1 ) t = ( n + 1 ) F s ( n ) ( t , 1 ) t + K ( n + 1 ) n F s ( n ) ( t , 1 ) , n = 1 , 2 , . Now, the following recurrence relation is yielded after integration: F s ( n + 1 ) ( t , 1 ) = ( n + 1 ) F s ( n ) ( t , 1 ) + K n ( n + 1 ) x = 0 t F s ( n ) ( x , 1 ) d x , n = 1 , 2 , . In particular, F ( t , 1 ) = 1 , F s ( t , 1 ) = 1 , F s ( t , 1 ) = 2 K t , F s ( 3 ) ( t , 1 ) = 3 ! ( K t + ( K t ) 2 ) .  □

Obviously, the obtained recurrent solution for factorial moments enables a form more convenient for computation polynomial, F s ( n ) ( t , 1 ) = n ! Q n 1 ( x ) , x = K t , Q 0 ( x ) = 1 , with the following recurrence relation Q n ( x ) = Q n 1 ( x ) + n Q n 1 ( x ) , Q 1 ( x ) = x , n = 1 , 2 , . Consequently, we obtain the remarkable property Q n ( x ) = 1 + 2 Q 1 ( x ) + + ( n 1 ) Q n 2 ( x ) + n Q n 1 ( x ) , n = 1 , 2 , . The polynomials Q n with their coefficients are denoted by Q n ( x ) = x + a n , 2 x 2 + a n , 3 x 3 + a n , 4 x 4 + + a n , n 1 x n 1 + x n . The coefficients a n , k are obtained recurrently as a direct result from (8), a n , k = a n 1 , k + n k a n 1 , k 1 , n = 1 , 2 , , k = 1 , 2 , , with the following initial and boundary values, a 0 , 0 = 1 , a n , 0 = a n , n + 1 = 0 , a n , 1 = a n , n = 1 , n = 1 , 2 , .

Summation by diagonals

The application of the Taylor theorem for the p.g.f. F ( t , s ) in the neighbourhood of s = 1, where F ( t , 1 ) = 1, yields the series expansion F ( t , s ) = 1 + n = 1 F s ( n ) ( t , 1 ) ( s 1 ) n n ! = 1 + n = 1 Q n 1 ( K t ) ( s 1 ) n . This series expansion (11) in the neighbourhood of the point s = 1 has a very little radius of convergence. This prevents the calculation of values of p.g.f. F ( t , s ) in the neighbourhood of zero by (11). To overcome this disadvantage we proceed with solution using summation of diagonal terms. The following convenient diagonal form for F ( t , s ) is obtained after placing (9) in (11) F ( t , s ) = 1 + ( s 1 ) + n = 2 ( s 1 ) n j = 1 n 1 a n 1 , j ( K t ) j . When the order of summation is exchanged, the p.g.f. takes the form F ( t , s ) = 1 + ( s 1 ) j = 0 ( s 1 ) j n = j a n , n j ( K t ( s 1 ) ) n j . The main diagonal is obtained for j = 0, where a 0 , 0 = 1, a n , n = 1, and ( s 1 ) n = 0 ( K t ( s 1 ) ) n = ( s 1 ) 1 K t ( s 1 ) . To avoid a confusion with the boundary values, the formula (12) is rewritten as follows: F ( t , s ) = 1 + ( s 1 ) 1 K t ( s 1 ) + ( s 1 ) j = 1 ( s 1 ) j n = j + 1 a n , n j ( K t ( s 1 ) ) n j . The summation by the second diagonal (the first after the main one) is computed for j = 1 and it is based on the recurrence relation (10), taking following form, a n , n 1 = a n 1 , n 1 + n n 1 a n 1 , n 2 , n = 2 , 3 , . It can be easily rewritten by using harmonic numbers H n  as follows, a n , n 1 = n ( H n 1 ) = n 1 2 + 1 3 + + 1 n . After some technical work, the summation on the second diagonal for | s 1 | < 1 / K t is generalised as ( s 1 ) 2 n = 2 a n , n 1 ( K t ( s 1 ) ) n 1 = s 1 1 K t ( s 1 ) 2 ( log ( 1 K t ( s 1 ) ) ) . Proceeding with the summation by the following diagonals, j 2, as it is shown in Appendix, leads us to a generalised solution for coefficients a n , n j and series expansion (13). Not surprisingly, working with diagonal forms implies the inclusion of combinatorial recurrent relations. In this case, the generalised solution and its triangular form involve the usage of the unsigned Stirling numbers of the first kind | s ( n , k ) | and equivalently of harmonic numbers . The complete implementation is summarised in the following lemma.

The coefficients of polynomials (9), Q n ( x ) = x + a n , 2 x 2 + a n , 3 x 3 + a n , 4 x 4 + + a n , n 1 x n 1 + x n , are expressed by the unsigned Stirling numbers of the first kind as follows: a n , n j = 1 ( n j ) ! k = 1 j + 1 ( 1 ) j + 1 k | s ( n + 1 , k ) | , n = j + 1 , j + 2 , .

The representation (14) is valid for the first and second diagonals ( j = 0 , j = 1 ), and for the first vertical column ( j = n 1 ), a n , n = 1 n ! | s ( n + 1 , n ) | = 1 , a n , 1 = 1 1 ! k = 1 n ( 1 ) n k | s ( n + 1 , k ) | = 1 , a n , n 1 = 1 ( n 1 ) ! { | s ( n + 1 , 1 ) | + | s ( n + 1 , 2 ) | } = 1 ( n 1 ) ! { n ! + n ! H n } . Moreover, for j = n 2 we have a n , 2 = 1 2 ! k = 1 n 1 ( 1 ) n 1 k | s ( n + 1 , k ) | = 1 2 ! { 1 + | s ( n + 1 , n ) | } = ( n + 2 ) ( n 1 ) 4 , where | s ( n + 1 , n ) | = n + 1 2 , see .

Now we verify whether the representation (14) yields the main recurrence relation (10) on the coefficients of the polynomials Q n ( x ), rewritten as follows, a n , n j = a n 1 , n j + n n j a n 1 , n 1 j , j = 1 , 2 , , n j > 0 . 0.\]]]> The proof is based on the main relation of the unsigned Stirling numbers of the first kind represented as | s ( n + 1 , k ) | = n | s ( n , k ) | + | s ( n , k 1 ) |. Now, it is computed from (14) and (15) that a n 1 , n j = 1 ( n j ) ! k = 1 j ( 1 ) j k | s ( n , k ) | = 1 ( n j ) ! k = 1 j + 1 ( 1 ) j + 1 k | s ( n , k 1 ) | , where the number of summands j = ( n 1 ) ( n j ) + 1. We remark that for k = 1, and n > 11$]]> the | s ( n , k 1 ) | = | s ( n , 0 ) | = 0. In the same way, n n j a n 1 , n 1 j = n n j 1 ( n 1 j ) ! k = 1 j + 1 ( 1 ) j + 1 k | s ( n , k ) | , where the number of summands j + 1 = ( n 1 ) ( n 1 j ) + 1. Then a n 1 , n j + n n j a n 1 , n 1 j = 1 ( n j ) ! k = 1 j + 1 ( 1 ) j + 1 k { | s ( n , k 1 ) | + n | s ( n , k ) | } . □ Solution of the Kolmogorov equation by the diagonal summation The computed coefficients a n , n j (14) allow to obtain the solution to the forward Kolmogorov equation of F ( t , s ). This can be done with a composite recurrent procedure of polynomials Q n . The sequence of iterations yields the following functions as a part of this solution: L ( t , s ) : = log ( 1 K t ( s 1 ) ) , A ( t , s ) : = ( s 1 ) 1 K t ( s 1 ) . Their consecutive derivatives by s (respectively by t) will be denoted as n F ( t , s ) s n = F s ( n ) ( t , s ) , n L ( t , s ) s n = L s ( n ) ( t , s ) , n A ( t , s ) s n = A s ( n ) ( t , s ) . Obviously, the following generalised forms are valid: L s ( t , s ) = K t 1 K t ( s 1 ) , L t ( t , s ) = K ( s 1 ) 1 K t ( s 1 ) = K A ( t , s ) , A s ( t , s ) = 1 ( 1 K t ( s 1 ) ) 2 , A t ( t , s ) = K ( s 1 ) 2 ( 1 K t ( s 1 ) ) 2 = K A 2 ( t , s ) . We remark also the following identities: A t ( t , s ) = K ( s 1 ) 2 A s ( t , s ) , L t ( t , s ) = ( s 1 ) L s ( t , s ) t . Now, the solution to the forward Kolmogorov equation (7) can be defined and proved by the following theorem. The p.g.f. of a Markov process with geometric branching F ( t , s ) = 1 + n = 0 ( s 1 ) 1 K t ( s 1 ) n + 1 Q n ( log ( 1 K t ( s 1 ) ) ) yields the forward Kolmogorov equation for | s 1 | < 1 / K t, F ( t , s ) t = K ( s 1 ) 2 2 s F ( t , s ) s . After denoting L : = L ( t , s ), A : = A ( t , s ), we rewrite (19) as F ( t , s ) = 1 + n = 1 { A ( t , s ) } n Q n 1 ( L ( t , s ) ) , Q 0 ( L ) = 1 . The first derivatives by t and s, respectively, in (20) are F t ( t , s ) = n = 1 { n A n 1 ( t , s ) Q n 1 ( L ) A t ( t , s ) + A n ( t , s ) Q n 1 ( L ) L t ( t , s ) } and F s ( t , s ) = n = 1 { n A n 1 ( t , s ) Q n 1 ( L ) A s ( t , s ) + A n ( t , s ) Q n 1 ( L ) L s ( t , s ) } . Then, taking into account the relations (16)–(18) and (8), the first derivatives by t and s are F t ( t , s ) = K A 2 ( t , s ) n = 1 A n 1 ( t , s ) Q n ( L ) = A t n = 1 A n 1 ( t , s ) Q n ( L ) and F s ( t , s ) = A s ( t , s ) n = 1 A n 1 ( t , s ) { n Q n 1 ( L ) + K t ( s 1 ) Q n 1 ( L ) } , or equivalently, F s ( t , s ) = A s ( t , s ) n = 1 A n 1 ( t , s ) { Q n ( L ) Q n 1 ( L ) + K t ( s 1 ) Q n 1 ( L ) } . The forward Kolmogorov equation can be written as ( 1 ( s 1 ) ) F t ( t , s ) = K ( s 1 ) 2 F s ( t , s ) . The left-hand side of (21) is K A 2 ( t , s ) n = 1 A n 1 ( t , s ) Q n ( L ) ( s 1 ) n = 1 A n 1 ( t , s ) Q n ( L ) . The right-hand side of (21) is K A 2 n = 1 A n 1 ( t , s ) Q n ( L ) ( 1 K t ( s 1 ) ) A ( t , s ) n = 1 A n 2 ( t , s ) Q n 1 ( L ) = K A 2 n = 1 A n 1 ( t , s ) Q n ( L ) ( s 1 ) n = 2 A n 2 ( t , s ) Q n 1 ( L ) , because, Q 0 = 1 and Q 0 = 0. Then, the change of the summation index gives n = 1 A n 1 ( t , s ) Q n ( L ) = n = 2 A n 2 ( t , s ) Q n 1 ( L ) . □ Statistical inferences The polynomial recurrent relation of the form (8), (9) is an important computation enhancement enabling a tabular representation, as is shown in Table 1. This recurrent relation converges quickly due to minimising of the fraction Q n / Q n 1 (see Figure 4). For practical implementations, only about n = 15 recurrent computations are required to obtain values for the p.g.f. F ( t , s ) as close as about 10 4 to the solution obtained through the Lambert-W function (see Figure 2). Q n ( x ) = Q n 1 ( x ) + n Q n 1 ( x )  Polynomials Q n Polynomials Q n ′ Q 0 ( x ) = 1 Q 0 ′ ( x ) = 0 Q 1 ( x ) = x Q 1 ′ ( x ) = 1 Q 2 ( x ) = x + x 2 Q 2 ′ ( x ) = 1 + 2 x Q 3 ( x ) = x + 5 x 2 2 + x 3 Q 3 ′ ( x ) = 1 + 5 x + 3 x 2 Q 4 ( x ) = x + 9 x 2 2 + 13 x 3 3 + x 4 Q 4 ′ ( x ) = 1 + 9 x + 13 x 2 + 4 x 3 Q 5 ( x ) = x + 14 x 2 2 + 71 x 3 3 ! + 77 x 4 12 + x 5 Q 5 ′ ( x ) = 1 + 14 x + 71 x 2 2 + 77 x 3 3 + 5 x 4 The graphs of the recurrent proportion Q n 1 / Q n for different K t The probabilities P ( X ( t ) = n ) of the critical geometric branching process can be computed directly by the usage of any available software tool for the Lambert-W function, as it is shown in . However, the statistical measures as variance, skewness, and kurtosis remain undefined relying on the solution based on the Lambert-W function. Their expressions follow directly from the representation of the n-th moments. It is given by the Stirling numbers of the second kind S ( n , k )  and factorial moments as E [ X ( t ) ] n = k = 1 n S ( n , k ) E [ X ( t ) ] k , n = 1 , 2 , . Thus, the variance of particles X ( t ) when E [ X ( t ) ] = 1 is V a r [ X ( t ) ] = E [ X ( t ) 1 ] 2 = E X ( t ) ( X ( t ) 1 ) E [ X ( t ) 1 ] = F s ( 2 ) ( t , 1 ) , t > 0 . 0.\]]]> The shape and skewness are also defined by the factorial moments: E [ ( X 1 ) 3 ] = E [ X ( X 1 ) ( X 2 ) ] + E [ X 1 ] = F s ( 3 ) ( t , 1 ) , and S k e w [ X ] = E [ ( X 1 ) 3 ] ( E [ ( X 1 ) 2 ] ) 3 / 2 = F s ( 3 ) ( t , 1 ) ( F s ( 2 ) ( t , 1 ) ) 3 / 2 . Finally, the function A ( t , s ) in the series expansion (19), (20) for geometric branching reproduction is of exactly twice greater intensity β = 2 K than the p.g.f. of the critical birth-death process, see . In this case the lifetime parameter K > 00$]]> (called killing rate) controls the intensity of reproduction and the size of population X ( t ) , t > 00\$]]>.

This connection conveniently enables comparison between different models for statistical inferences. For instance, it is not difficult to observe in parallel the models of critical geometric branching (GB), critical linear birth-death processes (BD) and the models with the number of progeny particles concentrated in ( 0 , 2 , 3 , 4 ). The reproduction p.g.f. for the new proposed models are h B D ( s ) = E [ s η ] = 1 + s 2 2 , h 3 ( s ) = 7 12 + s 2 4 + s 3 6 , h 4 ( s ) = 2 3 + s 2 8 + s 3 12 + s 4 8 . The intensities of reproduction are respectively K, β, α, γ. The corresponding infinitesimal generating functions and derivatives at s = 1 are, respectively, f G B ( s ) = K ( s 1 ) 2 2 s , f G B ( s ) = 2 ! K , , f G B ( n ) ( s ) = n ! K , f B D ( s ) = β ( s 1 ) 2 2 , f B D ( s ) = β , f B D ( 3 ) ( s ) = 0 , f 3 ( s ) = α 12 ( 1 s ) 2 ( 2 s + 7 ) , f 3 ( 1 ) = 3 α 2 , f 3 ( 3 ) ( 1 ) = α , f 3 ( 4 ) ( 1 ) = 0 . f 4 ( s ) = γ 24 ( 1 s ) 2 ( 3 s 2 + 8 s + 16 ) , f 4 ( 1 ) = 9 γ 4 , f 4 ( 3 ) ( 1 ) = 7 γ 2 , f 4 ( 4 ) ( 1 ) = 3 γ . In order to keep equal variances for all models, we take β = 2 K , α = 4 K / 3 , γ = 8 K / 9 . Under this condition, the measures of asymmetry are ordered as follows: S k e w ( B D ) < S k e w ( h 3 ) < S k e w ( h 4 ) < S k e w ( G B ) . The series expansion for reproduction defined by f 3 is approximated as F ( t , s ) = 1 + ( s 1 ) 1 3 α t 4 ( s 1 ) + 2 9 ( s 1 ) 1 3 α t 4 ( s 1 ) 2 log 1 3 α t 4 ( s 1 ) + . The explicit solution to the Kolmogorov equation for f 3 is F ( t , s ) = C 1 ( e 3 α t / 4 C ( s ) ) , C ( x ) = 2 x + 7 x 1 2 / 9 e 1 / ( 1 x ) , | 1 x | < 4 3 α t .

Furthermore, the intensities for similar branching processes could be computed using a higher order approximation of the p.g.f. F ( t , s ).

Conclusion

The continuity of the p.g.f. F ( t , s ) in the neighbourhood of the point s = 1 for critical branching process with geometric reproduction is proved analytically by the representations (19), (20) based on the factorial moments. This is an important extension to the already available solution based on the Lambert-W function and provides a computational tool for MBP with geometric reproduction.

Appendix

In this Appendix, we formally state and prove some formulas mentioned in the main text. This additional part is useful to give technical details for the computations of summation by diagonals introduced in Section 4.

The general combinatorics formulas are derived from the generating and exponential generating functions of the corresponding sequences . For example, log ( 1 z ) = n = 1 z n n = n = 1 | s ( n , 1 ) | z n n ! , log ( 1 z ) 1 z = n = 1 H n z n . The formula (28) in  for the series expansion of the Lambert-W function is based on the iterated logarithmic function L 2 : = log ( log ( x ) ) and the unsigned Stirling numbers of the first kind | s ( n , k ) |. It is natural that the same notions take place in the representations (19), (20).

The representation of the composite logarithmic and linear-fractional functions is given by the exponential partial Bell polynomials B n , k ( x ) over the sequence x = ( x 1 , x 2 , x n , ) as ( log ( 1 z ) ) r ( 1 z ) m + 1 = n = m j = 1 r B r , j ( x ) n m z n m . The sequence x is expressed by the powers of the harmonic numbers H n m and the generalised harmonic numbers H n ( m ) , H 0 ( m ) = 0, m = 1 , 2 , , , defined as follows: H n m = 1 + 1 2 + + 1 n m , H n ( m ) = 1 + 1 2 m + + 1 n m . The elements of the sequence x are x = ( { H n H m } , ( 1 ) 2 1 { H n ( 2 ) H m ( 2 ) } , , ( 1 ) r 1 ( r 1 ) ! { H n ( r ) H m ( r ) } ) . The relation between the Stirling and harmonic numbers reads : | s ( n + 1 , r + 1 ) | = n ! r ! k = 1 r B r , k ( H n , 1 ! H n ( 2 ) , , ( 1 ) r 1 ( r 1 ) ! H n ( r ) ) . The summation by the several first diagonals is obvious. To prove that the fourth diagonal is given by A 4 ( t , s ) Q 3 ( L ( t , s ) ), we calculate, applying (22), the series expansion of ( log ( 1 z ) ) 3 + 5 2 ( log ( 1 z ) ) 2 + ( log ( 1 z ) ) ( 1 z ) 4 . The coefficient a n , n 3 calculated from (23) and (14) is a n , n 3 = 1 ( n 3 ) ! { | s ( n + 1 , 4 ) | | s ( n + 1 , 3 ) | + | s ( n + 1 , 2 ) | | s ( n + 1 , 1 ) | } = n ( n 1 ) ( n 2 ) { H n 1 1 2 ( H n 2 H n ( 2 ) ) + 1 3 ! ( ( H n ) 3 3 H n H n 1 ( 2 ) + 2 H n ( 3 ) ) } . The comparison of the terms at z n 3 leads to the equality ( H n H 3 ) 3 3 ( H n H 3 ) ( H n ( 2 ) H 3 ( 2 ) ) + 2 ( H n ( 3 ) H 3 ( 3 ) ) + 5 2 { ( H n H 3 ) 2 ( H n ( 2 ) H 3 ( 2 ) ) } + ( H n H 3 ) = 3 ! ( H n 1 ) 3 ! 2 ( H n 2 H n ( 2 ) ) + { ( H n ) 3 3 H n H n ( 2 ) + 2 H n ( 3 ) } . The coefficients at H n 3 and H n ( 3 ) are evident. The coefficient at H n 2 is calculated as 3 H 3 + 5 / 2 = 3. The coefficient at H n ( 2 ) is calculated as 3 ( H n H 3 ) 5 / 2 = 3 3 H n .

At the end of combinatoric remarks, we cite the formula ( 28 ) of . When we denote W = W ( x ) , L 1 = log ( x ) , L 2 = log ( log ( x ) ) , the following expansion takes place: W = L 1 L 2 + k = 0 m = 1 ( 1 ) k | s ( k + m , k + 1 ) | m ! L 2 L 1 m 1 L 1 k . We remark that the software programme for numerical evaluating of the Lambert-W function is created basically on this formula.

References Adler, A.: lamW: Lambert-W Function (2015). https://CRAN.R-project.org/package=lamW Athreya, K.B., Ney, P.E.: Branching Processes. Springer, New York (1972) MR0373040 Choi, J., Srivastava, H.M.: Some summation formulas involving harmonic numbers and generalized harmonic numbers. Math. Comput. Model. 54, 22202234 (2011) MR2834625. https://doi.org/10.1016/j.mcm.2011.05.032 Connon, D.: Various applications of the (exponential) complete Bell polynomials. arXiv 2010, arXiv:1001.2835. 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. 197204. ACM, New York (1997) MR1809988. https://doi.org/10.1145/258726.258783 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, Hoboken (2005) MR2163227. https://doi.org/10.1002/0471715816 Merlini, D., Sprugnoli, R., Verri, M.C.: Lagrange inversion: when and how. Acta Appl. Math. 94, 233249 (2006) MR2290868. https://doi.org/10.1007/s10440-006-9077-7 Sevastyanov, B.A.: Branching Processes. Nauka, Moscow (1971) MR0345229 Spies, J.: Some identities involving harmonic numbers. Math. Comput. 55, 839863 (1990) MR1023769. https://doi.org/10.2307/2008451 Tchorbadjieff, A., Mayster, P.: Geometric branching reproduction Markov processes. Mod. Stoch. Theory Appl. 7(4), 357378 (2020) MR4195641. https://doi.org/10.15559/20-vmsta163 Tchorbadjieff, A., Mayster, P.: Models induced from critical birth death process with random initial conditions. J. Appl. Stat. 47(13-15), 28622878 (2020) MR4149585. https://doi.org/10.1080/02664763.2020.1732309 Wang, W., Wang, T.: General identities on Bell polynomials. Comput. Math. Appl. 58, 104118 (2009) MR2535972. https://doi.org/10.1016/j.camwa.2009.03.093