On Order Statistics from the Gompertz–Makeham Distribution and the Lambert W Function

Abstract The aim of this paper is twofold. First, we show that the expected value of the minimum order statistic from the Gompertz-Makeham distribution can be expressed in closed form in terms of the incomplete gamma function. We also give a general formula for the moments of the minimum order statistic in terms of the generalized integro-exponential function. As a consequence, the moments of all order statistics from this probability distribution can be more easily evaluated from the moments of the minimum order statistic. Second, we show that the maximum and minimum order statistics from the Gompertz-Makeham distribution are in the domains of attraction of the Gumbel and Weibull distributions, respectively. Lambert W function plays an important role in solving these problems.


Introduction
The Gompertz-Makeham distribution was introduced by the British actuary William M. Makeham [16] in the second half of the 19th century as an extension of the Gompertz distribution. Since then, this probability model has been successfully used in actuarial science, biology and demography to describe mortality patterns in numerous species, including humans. Marshall and Olkin [17,Chap. 10] provides a comprehensive review of the history and theory of this probability distribution and its practical importance is strongly highlighted in Golubev [9]. Some more recent mathematical contributions can also be found in Feng et al. [7], Jodrá [12], Lagerås [14] and Teimouri and Gupta [20].
To be more precise, let X be a continuous random variable having a Gompertz-Makeham distribution with positive real parameters α, β and λ, that is, the cumulative distribution function F (x; α, β, λ) := P (X ≤ x) is given by F (x; α, β, λ) = 1 − exp −λx − (α/β) e βx − 1 , x > 0, α, β, λ > 0. (1.1) The parameters α and β are usually interpreted as the initial mortality and the mortality increase when advancing age, respectively, whereas the parameter λ represents the risk of death due to causes which do not depend on age such as accidents and acute infections, among others.
Jodrá [12,Theorem 2] expresses the inverse of the cumulative distribution function of the Gompertz-Makeham distribution in closed form, more specifically, where log denotes the natural logarithm and W 0 denotes the principal branch of the Lambert W function, which is briefly described in the following paragraph. Throughout this paper, for notational simplicity we use F −1 (u) instead of the more cumbersome notation F −1 (u; α, β, λ).
For the sake of completeness, we recall that the Lambert W function is a multivalued complex function whose defining equation is W(z)e W(z) = z, where z is a complex number. This function has attracted a great deal of attention after the seminal paper by Corless et al. [6], which summarizes a wide variety of applications in mathematics and physics together with the main properties of W. In this regard, if z is a real number such that z ≥ −1/e then the Lambert W function has only two real branches. The real branch taking on values in [−1, ∞) (resp. (−∞, −1]) is called the principal (resp. negative) branch and denoted in the literature by W 0 (resp. W −1 ). Both real branches are depicted in Figure 1. In this paper we use only the principal branch W 0 and it is known that W 0 (z) is increasing as z increases, W 0 (−1/e) = −1, W 0 (0) = 0, W 0 (ze z ) = z and, in addition, the first derivative of W 0 is given by This paper is mainly motivated by the observation in Marshall and Olkin [17, p. 380] that the moments of the Gompertz-Makeham distribution, that is, E[X k ] := ∞ 0 x k dF (x; α, β, λ), k = 1, 2, . . . , cannot be given in closed form. Here, we show that the expected value of the minimum order statistic from the Gompertz-Makeham distribution can be expressed in closed form in terms of the incomplete gamma function. This result can be derived using Eq. (1.2) in order to highlight the importance of the Lambert W function. We also see that the kth moment of the minimum order statistic can be expressed explicitly in terms of the generalized integro-exponential function for any k = 1, 2 . . . . Then, as a particular case, we have a closed-form expression for the moments of the Gompertz-Makeham distribution, E[X k ], for any k = 1, 2, . . . . Additionally, using Eq. (1.2), we also determine the limit distributions of the maximum and minimum order statistics from this probability distribution. The remainder of this paper is organized as follows. In Section 2, we provide a closed-form expression for an integral involving the principal branch of the Lambert W function, specifically, 1 0 x p W 0 (s/x q ) dx, where p, q and s are positive real numbers. In Section 3, we see that the above integral arises when we compute the expected value of the minimum order statistic from the Gompertz-Makeham distribution and, as a consequence, the expected value can be expressed in closed form in terms of the incomplete gamma function. In this section, we also see that the moments of the minimum order statistic can be expressed in terms of the generalized integro-exponential function. These explicit expressions corresponding to the minimum are useful to evaluate the moments of all order statistics from the Gompertz-Makeham distribution. A numerical example illustrates the results. Finally, in Section 4, we identify the asymptotic behavior of the maximum and minimum order statistics from the Gompertz-Makeham distribution. More precisely, we establish that the domains of attraction corresponding to the maximum and minimum order statistics are the Gumbel and Weibull distributions, respectively.

An Integral Involving the Principal Branch of the Lambert W Function
In the following result, we provide a closed-form expression for an integral involving the principal branch of the Lambert W function. Denote by Γ (a, z) the (upper) incomplete gamma function (cf. Abramowitz and Stegun [1, p. 260]), that is, Lemma 1. For any positive real number q, we have Proof. First, we make the change of variable u = 1/x q in the integral in Eq. (2.2) and we obtain By setting w = W 0 (u) on the right-hand side in the above equation, which implies u = w e w and also du = (1 + w)e w dw by virtue of Eq. (1.3), we get Now, on the right-hand side in Eq. (2.3) we take into account Eq. (2.1) together with the fact that Γ (a, ∞) = 0, and thus we obtain It is well-known that the incomplete gamma function satisfies the recurrence formula Γ (a+1, z) = aΓ (a, z)+z a e −z (see Abramowitz where in the last equality we have used the fact that W 0 (z) e W0(z) = z. This completes the proof.
Essentially the same reasoning as in Lemma 1 leads to the result stated in Lemma 2 below, so we omit the proof here.
Lemma 2. For any positive real numbers p, q and s, we have In the next section, we apply Lemma 2 to obtain a closed-form expression for the expected value of the minimum order statistic from the Gompertz-Makeham distribution.

Expected Value of Order Statistics from the Gompertz -Makeham Distribution
We first introduce some notation and terminology. Denote by N the set of non-negative integers and by N * := N {0}. For any n ∈ N * , let X 1 , . . . , X n be a random sample of size n from the Gompertz-Makeham distribution with parameters α, β and λ, and let X 1:n , X 2:n , . . . , X n:n be the corresponding order statistics obtained by arranging X i , i = 1, . . . , n, in non-decreasing order of magnitude. The rth element of this sequence, X r:n , is called the rth order statistic and, in particular, the minimum and maximum order statistics, that is, X 1:n = min{X 1 , . . . , X n } and X n:n = max{X 1 , . . . , X n }, the so-called extremes, are very important in practical applications. For any n ∈ N * and k ∈ N * , it is known that the kth moment of X r:n , r = 1, . . . , n, can be computed as follows (cf. Balakrishnan and Rao [3, p. 7]) For the Gompertz-Makeham distribution, even the problem of computing the expected value E[X r:n ] by direct numerical integration of Eq. (3.1) is not efficient from a computational viewpoint, and, in particular, the problem becomes intractable when the sample size n is relatively large (see the example at the end of this section).
Our first aim in this section is to compute the expected value E[X r:n ] in a more efficient manner, that is, avoiding the numerical integration of Eq. (3.1), and, to this end, the minimum order statistic X 1:n plays a crucial role. More precisely, the expected value of X 1:n can also be obtained by means of the following formula involving F −1 (cf. Arnold et al. [2,Chapter 5] for further details) where in our case F −1 is given by Eq. (1.2). By using Eq. (3.2) together with Lemma 2, in the following result we give an analytical expression for the expected value of the minimum order statistic X 1:n , for n = 1, 2, . . . . The special case n = 1 corresponds to the expected value of X.
. . , X n be a random sample of size n from a Gompertz-Makeham distribution with parameters α, β and λ. Then, the expected value of X 1:n is given by Proof. For any n ∈ N * , from Eq. (3.2) together with Eq. (1.2) we have Now, by taking into account Lemma 2 in the above equation, we get where in the last equality we have used the fact that W 0 (ze z ) = z. The proof is completed.
Here, we point out the importance of the closed-form expression given in Corollary 1. On the one hand, from a theoretical point of view, the expected value of the minimum order statistic can be used to determine if two distributions with finite expected values are identical (cf. Chan [5] and also Huang [11]). On the other hand, from a numerical point of view, in order to calculate E[X 1:n ] by virtue of Corollary 1, computer algebra systems such as Maple or Mathematica, among others, can be used to directly evaluate the incomplete gamma function. Moreover, the analytical expression of E[X 1:n ] in Corollary 1 can also be used to compute E[X r:n ], for r = 2, . . . , n, avoiding the numerical integration of Eq. (3.1). To this end, for any n ∈ N * and k ∈ N * we can use the following relation (cf. Balakrishnan and Rao [3, p. 156 By setting k = 1 in Eq. (3.3) and in view of Corollary 1, we get the following result.
Corollary 2. Let X 1 , . . . , X n be a random sample of size n from a Gompertz-Makeham distribution with parameters α, β and λ. Then, for r = 2, . . . , n, we have Unfortunately, we have not found simple analytical expressions for the moments E[X k 1:n ] for integers k ≥ 2. However, we provide a general formula for the moments E[X k 1:n ], which can be used for computing these moments for integers k ≥ 2 in a more efficient manner than using Eq. (3.1). We first give the following preliminary result.
Lemma 3. Let X 1 , . . . , X n be a random sample of size n from a Gompertz-Makeham distribution with parameters α, β and λ. Then, the minimum order statistic X 1:n has a Gompertz-Makeham distribution with parameters nα, β and nλ.

P. Jodrá
Proof. For any n ∈ N * , recall that the cumulative distribution function of X 1:n , F 1:n (x; α, β, λ) := P (X 1:n ≤ x), can be computed as follows Then, by using Eq. (1.1) in the above formula we obtain that is, F 1:n (x; α, β, λ) = F (x; nα, β, nλ) which implies the result. Now, we are in a position to express the moments E[X k 1:n ] in terms of the generalized integro-exponential function. We recall that the generalized integro-exponential function can be defined by the following integral representation (cf. Milgram [18] for further details) where m = 0, 1, . . . and s is a real number.
Proposition 1. Let X 1 , . . . , X n be a random sample of size n from a Gompertz -Makeham distribution with parameters α, β and λ. Then, for the minimum order statistic X 1:n we have Proof. By virtue of Lemma 3, together with Eq. (1.1), we get Now, the change of variable u = e βx leads to the following From the above equation and the definition in Eq. (3.4), for any k ∈ N * we have Finally, by taking into account the following recursion formula (cf. Milgram  It is interesting to note that if we take k = 1 in Proposition 1, and using the fact that E 0 s (z) = z s−1 Γ (1 − s, z) (cf. Milgram [18, Eq. (2. 2)]), we have an alternative derivation of Corollary 1. Moreover, it is clear that if we take n = 1 in Proposition 1, then we have an explicit expression for the moments of the Gompertz-Makeham distribution, E[X k ], for any k = 1, 2, . . . . We highlight that Lenart [15] has recently given explicit expressions for the moments of the Gompertz distribution in terms of the generalized integro-exponential function. In fact, if we take n = 1 and λ = 0 in Eq. (3.5), then we obtain the same expression provided by Lenart for the moments of the Gompertz distribution.
As a consequence of Proposition 1 together with Eq. (3.3), we give a general formula for the moments of order statistics from the Gompertz-Makeham distribution, as stated below.
Corollary 3. Let X 1 , . . . , X n be a random sample of size n from a Gompertz-Makeham distribution with parameters α, β and λ. Then, for any k ∈ N * and for r = 2, . . . , n, we have To complete this section, we provide an example in which we compute the two first moments corresponding to order statistics from the Gompertz-Makeham distribution. Unless differently specified, all of the computations in the example below were performed with the software Mathematica 8.0, on an Intel Core2 Quad Q8200 at 2.33GHz with 4GB RAM.
Numerical Example. Promislow and Haselkorn [19] have studied the aging process for different species of flies. In particular, it is considered a sample of size n = 990 from a Gompertz-Makeham distribution X with parameters α = 0.00049, β = 0.071 and λ = 0.00092. In this case, the random variable X represents the lifetime (in days) of females of Drosophila melanogaster (fruit fly).
With the help of the computer algebra system Mathematica and applying Corollary 1, in a preprocessing step we have first calculated and stored the expected values E[X

Domains of Attraction of the Extreme Order Statistics
Let X 1 , . . . , X n be a random sample of size n from a Gompertz-Makeham distribution with parameters α, β and λ. Let us consider the maximum and minimum order statistics, X n:n and X 1:n , respectively, and denote by F n:n and F 1:n the corresponding cumulative distribution functions. It is known that F n:n (x) = F n (x) and F 1:n (x) = 1 − (1 − F (x)) n , x > 0, where F is given by Eq. (1.1). For notational simplicity, in this section we use F (x) instead of F (x; α, β, λ) (similarly for F 1:n and F n:n ). In order to identify the asymptotic distributions of X n:n and X 1:n when the sample size n increases to ∞, it is well-known that the limits of F n:n and F 1:n take only values 0 and 1, which means that the limit distributions are degenerate. To avoid degeneracy, it is necessary to look for linear transformations to find the asymptotic non-degenerate distributions. If there exist non-degenerate cumulative distribution functions H and L such that lim n→∞ F n (a n + b n x) = H(x), lim where a n , b n > 0, c n and d n > 0 are normalization constants depending on the sample size n, then it is said that F belongs to the maximum -resp. minimumdomain of attraction of the limit distribution H -resp. L-. Moreover, it is wellestablished in the statistical literature that the limit distributions H and L can only be either a Fréchet, a Gumbel or a Weibull distribution. For further details on the asymptotic theory of extremes see, for instance, Castillo at al. [4,Chap. 9], Kotz and Nadarajah [13, Chap. 1] and also Galambos [8].
In this section, we determine the limit distribution of the extremes of the Gompertz-Makeham distribution, that is, we identify H and L in Eqs. (4.1) together with the normalizing constants. With the previous notations, below we state the result and the proof is given in the remainder of this section.
The proof of Theorem 1 is based on applying Theorems 9.5 and 9.6 given in Castillo et al. [4, pp. 203-205]. Both theorems require that the inverse of the cumulative distribution function corresponding to the continuous random variable involved can be expressed analytically and, by virtue of Eq. (1.2), this is the case for the Gompertz-Makeham distribution.
Before proceeding with the proof of Theorem 1, we need some auxiliary results concerning asymptotic properties of the principal branch of the Lambert W function. Denote by where log denotes the natural logarithm. Hoorfar and Hassani [10, Theorem 2.7] have shown that L W0 (x) and U W0 (x) are bounds for W 0 (x), more specifically, where the equality is attained only in the case x = e. Now, we establish the following results. Proof. For any positive real numbers a and b, from Eq. (4.4) we can bound the difference W 0 (ax) − W 0 (bx) as follows

P. Jodrá
As a consequence, the statement of this lemma is obtained by applying the squeeze theorem to Eq. (4.5).
Lemma 5. For any a > 0, we have Proof. For any positive real number a, we can bound (1+W 0 (ax))/(1+W 0 (x)) by using Eq. (4.4) as follows It can be seen that lim x→∞ U W0 (ax) = ∞ and also that lim x→∞ L W0 (ax) = ∞. Then, we can compute lim x→∞ (1 + U W0 (ax))/(1 + L W0 (x)) via l'Hôpital's rule. Denote by U W0 (resp. L W0 ) the first derivative of U W0 (resp. L W0 ). In addition, denote by After a bit of algebra, U W0 (ax)/L W0 (x) can be written as below Thereby, we have On the other hand, by similar arguments, it can also be verified that Finally, from Eqs. (4.8) and (4.9) the statement of this lemma is obtained by applying the squeeze theorem to Eq. (4.6).
With the previous results, we are now in a position to prove Theorem 1.
We can now apply l'Hôpital's rule to obtain lim ε→0 H 1 (ε)/H 2 (ε). Denote by H i the first derivative of H i , i = 1, 2. Then, we have where W 0 denotes the first derivative of W 0 with respect to ε and the second equality in Eq. (4.13) is obtained by taking into account Eq. (1.3). Finally, we apply Lemmas 4 and 5 to calculate the limit in Eq. (4.13) and then we get lim ε→0 H 1 (ε)/H 2 (ε) = 1, which implies that Eq. (4.10) holds. The normalizing constants a n and b n given in Eq. where F −1 is given by Eq. (1.2). We claim that c = 1 and, if this is the case, then the Gompertz-Makeham distribution belongs to the minimum domain of attraction of the Weibull distribution. To see this, we denote by L k (ε) := P. Jodrá F −1 (kε) − F −1 (2kε), k = 1, 2. We need to compute lim ε→0 L 1 (ε)/L 2 (ε). From Eq. (1.2), we have L k (ε) = 1 λ log 1 − 2kε 1 − kε where α, β and λ are positive real parameters. Now, by taking into account that W 0 (xe x ) = x it is easy to check that lim ε→0 L k (ε) = 0 for k = 1, 2. Then, we can apply l'Hôpital's rule to compute lim ε→0 L 1 (ε)/L 2 (ε), so that, denoting by L i the first derivative of L i , i = 1, 2, we have to compute where (F −1 ) denotes the first derivative of F −1 with respect to ε. For k = 1, 2, 4, it can be checked the following . which implies that c = 1 in Eq. (4.14), as it was claimed. Finally, the normalizing constants c n and d n in Eq. (4.3) are chosen according to Castillo et al. [4,Theorem 9.6]. This concludes the proof of part (ii). The proof of Theorem 1 is completed.