A branching process model driven by geometric reproduction of particles was introduced in [12]. 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 (1) 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 [12] as the unique solution to the Kolmogorov equations [2, 7, 10] after the Lagrange inversion method following the classical theory [5, 8, 9].

The factorial moments in the subcritical case are computed in [12]. 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 [12].

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 [12]. 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.

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.

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, [2, 7, 10]. 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 (2) F ( t , s ) = ∑ k = 0 ∞ s k P ( X ( t ) = k | X ( 0 ) = 1 ) satisfies the backward Kolmogorov equation (3) ∂ ∂ 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 (4) ∫ 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 (5) 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 composite inverse function of V ( G ( s ) ) is obtained in [12] 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 (6) 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$]]>, [10]. 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}}).\]]]>

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, (7) ∂ 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.

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 (8) 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 (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 . The coefficients a n , k are obtained recurrently as a direct result from (8), (10) 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 , … .

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 (11) 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 (12) 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: (13) 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 [3] 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 [14, 4]. The complete implementation is summarised in the following lemma. Lemma 2.

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: (14) 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 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).

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 [12]. 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 ) [14, 4] 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 .