A new lifetime distribution - Semantic Scholar
Computational Statistics & Data Analysis 51 (2007) 4497 – 4509 www.elsevier.com/locate/csda
A new lifetime distribution Co¸skun Ku¸s∗ Department of Statistics, Faculty of Art and Science, Selcuk University, 42031, Konya, Turkey Received 10 November 2005; received in revised form 26 June 2006; accepted 12 July 2006 Available online 4 August 2006
Abstract A new two-parameter distribution with decreasing failure rate is introduced. Various properties of the introduced distribution are discussed. The EM algorithm is used to determine the maximum likelihood estimates and the asymptotic variances and covariance of these estimates are obtained. Simulation studies are performed in order to assess the accuracy of the approximation of the variances and covariance of the maximum likelihood estimates and investigate the convergence of the proposed EM scheme. Illustrative examples based on real data are also given. © 2006 Elsevier B.V. All rights reserved. Keywords: Compounding; Decreasing failure rate; EM algorithm; Exponential distribution; Lifetime distributions; Maximum likelihood estimation; Zero truncated Poisson distribution
1. Introduction In general, a population is expected to exhibit decreasing failure rate (DFR) when its behavior over time is characterized by ‘work-hardening’ (in engineering terms) or ‘immunity’ (in biological terms); sometimes the broader term ‘infant mortality’ is used to denote the DFR phenomenon (Adamidis and Loukas, 1998). The distributions with DFR are discussed in the works of Lomax (1954), Proschan (1963), Barlow et al. (1963), Barlow and Marshall (1964, 1965), Marshall and Proschan (1965), Cozzolino (1968), Dahiya and Gurland (1972), McNolty et al. (1980), Saunders and Myhre (1983), Nassar (1988), Gleser (1989), Gurland and Sethuraman (1994) and Adamidis and Loukas (1998). This paper is organized as follows: In Section 2, a new DFR distribution is obtained by mixing exponential and zero truncated Poisson distribution where mixing procedure was previously carried out by Adamidis and Loukas (1998). In Section 3, various properties of the introduced distribution are discussed. In Section 4, the estimation of parameters is studied by the method of maximum likelihood. Finally, in Section 5, illustrative examples based on real data are provided to close the paper. 2. The distribution Let W1 , W2 , . . . , WZ be a random sample from distribution with density f (w; ) = e−w , w, ∈ R+ , Z is a −1 , z ∈ N, ∈ R+ zero truncated Poisson variable with probability function P (z; ) = e− z −1 (z + 1) 1 − e− ∗ Tel.: +90 3322231354; fax: +90 3322410106.
E-mail address: [email protected]. 0167-9473/$ - see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.csda.2006.07.017
4498
C. Ku¸s / Computational Statistics & Data Analysis 51 (2007) 4497 – 4509
.
.
.
.
. .
.
.
.
Fig. 1. Probability density functions of the EP distribution for = 1 and = 1, 2, 3, identified, respectively, by the y-axis intercepts appearing in increasing order of magnitude.
.
.
.
.
.
.
.
.
Fig. 2. Probability density functions of the EP distribution for = 1 and = 1, 2, 3, identified, respectively, by the y-axis intercepts appearing in increasing order of magnitude.
where (·) is the Gamma function and Z and W s are independent. Let us define X = min (W1 , W2 , . . . , WZ ). Then, f (x|z; ) = ze−zx and marginal probability density function of X is e−−x+ exp(−x) , 1 − e−
f (x; ) =
x, , ∈ R+ ,
(1)
where = (, ). In the sequel, distribution of X will be referred to as the exponential–Poisson distribution (EP) which is customary for such names to be given to the distributions arising via the operation of compounding (mixing) in −1 the literature. It can be seen that the EP density function is monotone decreasing with modal value 1 − e− at x = 0. EP probability density functions are displayed in Figs. 1 and 2 for selected parameter values. For all values of parameters, the density is strictly decreasing in x and tending to zero as x → ∞. Its graph resembles those of the EG, Pareto II and exponential distributions. As approaches zero, the EP leads to exponential distribution with parameter . The model is obtained under the concept of population heterogeneity (through the process of compounding). An interpretation of the proposed model is as follows: a situation where failure (of a device for example) occurs due to the presence of an unknown number, Z, of initial defects of same kind (a number of semiconductors from a defective lot, for example). The Ws represent their lifetimes and each defect can be detected only after causing failure, in which case
C. Ku¸s / Computational Statistics & Data Analysis 51 (2007) 4497 – 4509
4499
it is repaired perfectly (Adamidis and Loukas, 1998). Then the distributional assumptions given earlier lead to the EP distribution for modelling the time to the first failure X. 3. Properties of the distribution 3.1. Distribution function and moments The distribution function is given by −1 F (x; ) = e exp(−x) − e 1 − e
(2)
−1 −1 . and hence, the median is obtained by log log 2−1 e + 1 Since the EP is not part of the exponential family, the evaluations of moments become cumbersome. The raw moments of X are determined from (1) by direct integration. For r ∈ N, raw moments are given by (r + 1) Fr+1,r+1 ([1, 1, . . . , 1], [2, 2, . . . , 2], ), E Xr ; = e − 1 r
r ∈ N,
where Fp,q (n, d, ) is the generalized hypergeometric function. This function is also known as Barnes’s extended hypergeometric function. The definition of Fp,q (n, d, ) is Fp,q (n, d, ) =
∞
k=0
(k
p
−1 (ni ) i=1 (ni + k) , q + 1) i=1 (di + k) −1 (di )
k
where n = [n1 , n2 , . . . , np ], p is the number of operands of n, d = [d1 , d2 , . . . , dq ] and q is the number of operands of d. Generalized hypergeometric function is quickly evaluated and readily available in standard software such as Maple. Hence the mean and variance of the EP distribution are given, respectively, by F2,2 ([1, 1], [2, 2], ), E(X; ) = e −1
2 F ([1, 1], [2, 2], ) . 2F3,3 ([1, 1, 1], [2, 2, 2], ) − Var(X; ) = e − 1 2,2 e − 1 2
(3)
3.2. The survivor, hazard and mean residual life functions Using (1) and (2), survival function (also known reliability function) and hazard function (also known failure rate function) of the EP distribution are given, respectively, by −1 s(x; ) = 1 − F (x; ) = 1 − e exp(−x) 1 − e , −−x+ exp(−x) 1 − e f (x; ) e . h(x; ) = = s(x; ) 1 − e− 1 − e exp(−x) Both functions have simple forms in contrast to those of some DFR distributions such as Gamma with DFR. The hazard function is decreasing because the DFR property follows from the results of Barlow et al. (1963) on mixtures. The EP hazard functions for selected parameter values are displayed in Fig. 3. −1 and h(∞; ) = , respectively. They are both finite The initial and long-term hazards are h(0; ) = e e − 1 in contrast to those of Weibull distribution with h(0; ) = ∞ and h(∞; ) = 0.
4500
C. Ku¸s / Computational Statistics & Data Analysis 51 (2007) 4497 – 4509
.
.
.
.
.
.
.
Fig. 3. Hazard functions of the EP distribution for = 1 and = 0.5, 1.0, 1.5, identified, respectively, by the y-axis intercepts appearing in increasing order of magnitude.
The mean residual lifetime of the EP distribution is given by m (x0 ; ) = E (X − x0 |X x0 ; ) =
e−x0 F2,2 [1, 1], [2, 2], e−x0 , e exp(−x0 ) − 1
where Fp,q (n, d, ) is the generalized hypergeometric function. Because the EP is a DFR distribution m (x0 ; ) increases monotonically up to its limit at infinity. Clearly, m (x0 ; ) → E(X; ) as x0 → 0 and m (x0 ; ) → −1 as x0 → ∞. 4. Estimation of the parameters 4.1. Estimation by maximum likelihood The log-likelihood function based on the observed sample size of n, yobs = (xi ; i = 1, 2, . . . , n), from the EP distribution is given by (; yobs ) = n log() −
n
xi − n +
i=1
n
e−xi − n log 1 − e−
i=1
and subsequently the associated gradients are found to be n −1 j (; yobs ) −xi , = e + n −1 − 1 − e− j
(4)
i=1
j (; yobs ) = n−1 − xi − xi e−xi . j n
n
i=1
i=1
The equation j (; yobs ) /j = 0 could be solved exactly for , namely
n −1 1
−1 Xi exp −Xi = −X n
(5)
(6)
i=1
conditional upon the value of , where and are the maximum likelihood estimates (MLEs) for the parameters and , respectively. Eq. (6) might be used for building a profile likelihood. In the following, Theorems 1 and 2 express when the other parameter is given (or known), what the conditions are there to obtain the existence and uniqueness of the MLE.
C. Ku¸s / Computational Statistics & Data Analysis 51 (2007) 4497 – 4509
4501
in (4), where is the true value of the Theorem 1. Let h (; , yobs ) denote the function on the RHS of the expression , is unique for n/2 < ni=1 e−xi . parameter. Then, the root of h (; , yobs ) = 0, Theorem 2. Let g (; , yobs ) denote the function on the RHS of the expression in (5) and x = n−1 ni=1 xi , where 2 , lies in the interval is the true value of the parameter. Then, for a given ∈ 0, e the root of g (; , yobs ) = 0, −1 −1 x( + 1) , x and is unique. 4.2. An EM algorithm The MLEs of and must be derived numerically. Newton–Raphson algorithm is one of the standard methods to determine the MLEs of the parameters. To employ the algorithm, second derivatives of the log-likelihood are required for all iteration. EM algorithm is a very powerful tool in handling the incomplete data problem (Dempster et al., 1977; McLachlan and Krishnan, 1997). It is an iterative method by repeatedly replacing the missing data with estimated values and updating the parameter estimates. It is especially useful if the complete data set is easy to analyze. As pointed out by Little and Rubin (1983), the EM algorithm will converge reliably but rather slowly (as compared to the Newton–Raphson method) when the amount of information in the missing data is relatively large. Recently, EM algorithm has been used by several authors such as Adamidis and Loukas (1998), Adamidis (1999), Ng et al. (2002), Karlis (2003) and Adamidis et al. (2005). To start the algorithm, hypothetical complete-data distribution is defined with density function f (x, z; ) = −1 , x > 0, z = 1, 2, . . . , , > 0. Thus, it is straightforward to verify that the Eze−zx e− z −1 (z + 1) 1 − e− step of an EM cycle requires the computation of the conditional expectation of Z|X; (h) , where (h) = (h) , (h) is the current estimate of . Using that P (z|x; ) = e−xz− exp(−x)+x z−1 −1 (z), z ∈ N, this is found to be E(Z|X; ) = 1 + e−x . The EM cycle is completed with M-step, which is complete data maximum likelihood over , with the missing Z’s replaced by their conditional expectations E(Z|X; ) (Adamidis and Loukas, 1998). Thus, an EM iteration, taking (h) into (h+1) , is given by
(h+1)
=n
n
−1
xi 1 +
(h) −(h) xi
e
,
i=1
(h+1)
=n
−1
1−e
−(h+1)
n
1 + (h) e−
(h)
xi
.
(7)
i=1
It can be seen that only a one-dimensional search such as Newton–Raphson is required for M-step of an EM cycle. Ng et al. (2002) used the same method for estimation of the parameters of the Weibull distribution based on progressively type-II right censored sample. 4.3. Asymptotic variances and covariance of the MLEs Applying the usual large sample approximation, the MLE of can be treated as being approximately bivariate normal with mean and variance–covariance matrix, which is the inverse of the expected information matrix J () = E(I ; ), where I = I (; yobs ) is the observed information matrix with elements Iij = −j2 /ji jj with i, j = 1, 2 and the expectation is to be taken with respect to the distribution of X. Differentiating (4) and (5), the elements of the symmetric, second-order observed information matrix are found to be −2 , I11 = n 1 + e2 − 2 e − 2e 1 − e− I12 =
n
i=1
xi e
−xi
,
I22 = n
−2
−
n
i=1
xi2 e−xi .
4502
C. Ku¸s / Computational Statistics & Data Analysis 51 (2007) 4497 – 4509
The elements of the expected information matrix, J (), are calculated by taking the expectations of Iij , i, j = 1, 2 with respect to the distribution of X, i.e. following expectations are required: E Xe−X =
e− F2,2 ([2, 2], [3, 3], ), 4 1 − e−
E X2 e−X =
e− F3,3 ([2, 2, 2], [3, 3, 3], ). 4 1 − e− 2
Thus J () are derived in the form −2 , J11 = n 1 + e2 − 2 e − 2e 1 − e− J12 =
ne− F2,2 ([2, 2], [3, 3], ), 4 1 − e−
J22 = n−2 −
n2 e− F3,3 ([2, 2, 2], [3, 3, 3], ). 42 1 − e−
The inverse of J (), evaluated at provides the asymptotic variance–covariance matrix of the MLEs. Alternative estimates can be obtained from the inverse of the observed information matrix since it is a consistent estimator of J −1 ().
4.4. Simulation study Simulation results are obtained based on the assumption of Theorem 1, to guarantee the existence and uniqueness of the solution in each iteration of (h+1) in EM algorithm. No restriction has been imposed on the maximum number of iterations and convergence is assumed when the absolute differences between successive estimates are less than 10−5 . For each value of , the parameter estimates have been found by employing (7) with different initial values. In order to assess the accuracy of the approximation of the variances and covariance of the MLEs determined from the information a simulation matrix, study (based on 10 000 simulations) has been carried out. The simulated values of Var , Var and Cov , as well as the approximate values determined by averaging the corresponding values obtained from the expected and observed information matrix are presented in Table 1. From Table 1, it is observed that the approximate values determined from expected and observed information matrix are quite close to the simulated values for large values of n. Furthermore, it is noted that the approximation becomes quite accurate as n increases. In addition, variances and covariances of the MLEs obtained from the observed information matrix is quite close to the that of expected information matrix for large values of n. Also simulations have been performed to investigate the convergence of the proposed EM scheme. Ten thousand samples of size 100 and 500 each of which are randomly sampled from the EP for each of five values of are generated. The results from simulated data sets are reported in Table 2, which gives the averages of the 10 000 MLEs av and averages numbers of iterations to convergence av(t) together with standard errors se and se(t). Table 2 indicates the following results: convergence has been achieved in all cases, even when the starting values are poor and this emphasizes the numerical stability of the EM algorithm. The values of av and se suggest the EM estimates performed consistently. Standard errors of the MLEs decrease when sample size increases. MLEs of the parameters and have bias. However, bias of the MLEs are reduced for large sample size in five cases of parameters in Table 2.
C. Ku¸s / Computational Statistics & Data Analysis 51 (2007) 4497 – 4509
4503
Table 1 Variances and covariance of the MLEs
n
Simulated Var
Var
, Cov
From expected information , Var Var Cov
From observed information , Var Var Cov
50
(0.5, 1) (1, 0.5) (0.5, 2) (2, 0.5) (2, 2)
1.5599 1.8175 1.6120 2.2499 2.2971
0.0662 0.0242 0.2688 0.0502 0.7984
−0.2564 −0.1715 −0.5273 −0.2817 −1.1291
3.9604 5.2131 4.0185 8.5655 8.4945
0.1507 0.0545 0.5990 0.1178 1.8565
−0.6135 −0.4319 −1.2233 −0.8435 −3.3073
3.1078 3.5839 3.1327 4.5887 4.4832
0.1291 0.0422 0.5208 0.0744 1.1730
−0.4909 −0.3052 −0.9892 −0.4571 −1.7868
100
(0.5, 1) (1, 0.5) (0.5, 2) (2, 0.5) (2, 2)
1.0886 1.4367 1.1664 2.0511 2.1500
0.0429 0.0173 0.1747 0.0385 0.6070
−0.1781 −0.1352 −0.3744 −0.2477 −1.0097
1.3737 2.0329 1.4191 4.2078 4.2756
0.0661 0.0245 0.2655 0.0566 0.8995
−0.2389 −0.1811 −0.4825 −0.4132 −1.6525
1.3121 1.8245 1.2934 2.7439 2.8582
0.0651 0.0224 0.2579 0.0404 0.6671
−0.2344 −0.1605 −0.4625 −0.2685 −1.1210
500
(0.5, 1) (1, 0.5) (0.5, 2) (2, 0.5) (2, 2)
0.1126 0.2821 0.1119 0.9741 0.9941
0.0098 0.0049 0.0391 0.0139 0.2241
−0.0289 −0.0330 −0.0574 −0.1080 −0.4383
0.1333 0.2063 0.1328 0.7133 0.7222
0.0118 0.0045 0.0471 0.0116 0.1851
−0.0355 −0.0273 −0.0708 −0.0811 −0.3263
0.1325 0.2145 0.1323 0.7444 0.7630
0.0118 0.0045 0.0469 0.0114 0.1825
−0.0354 −0.0280 −0.0706 −0.0082 −0.3292
1000
(0.5, 1) (1, 0.5) (0.5, 2) (2, 0.5) (2, 2)
0.0570 0.0949 0.0575 0.5030 0.4862
0.0054 0.0022 0.0215 0.0074 0.1176
−0.0155 −0.0135 −0.0311 −0.0573 −0.2239
0.0631 0.0919 0.0631 0.3081 0.3055
0.0058 0.0022 0.0231 0.0059 0.0941
−0.0171 −0.0131 −0.0342 −0.0393 −0.1570
0.0633 0.0918 0.0636 0.3437 0.3528
0.0058 0.0022 0.0232 0.0061 0.0977
−0.0171 −0.0131 −0.0344 −0.0422 −0.1700
10 000
(0.5, 1) (1, 0.5) (0.5, 2) (2, 0.5) (2, 2)
0.0061 0.0086 0.0062 0.0236 0.0233
0.0006 0.0002 0.0023 0.0006 0.0092
−0.0017 −0.0012 −0.0034 −0.0035 −0.0142
0.0062 0.0087 0.0062 0.0232 0.0237
0.0006 0.0002 0.0023 0.0006 0.0093
−0.0017 −0.0013 −0.0034 −0.0036 −0.0143
0.0062 0.0087 0.0062 0.0237 0.0237
0.0006 0.0002 0.0023 0.0006 0.0093
−0.0017 −0.0013 −0.0034 −0.0036 −0.0143
5. Illustrative examples The fit of the EP distribution of real data is examined by graphical methods using MLEs. It is also compared with the EG, Weibull and Gamma models with respective densities −2 f1 (x; 1 ) = 1 1 − 1 e−1 x 1 − 1 e−1 x , x > 0, 1 > 0, 1 ∈ (0, 1),
2
f2 (x; 2 ) = 2 2 2 x 2 −1 e−(2 x ) ,
x > 0, 2 , 2 > 0,
−1 f3 (x; 3 ) = 3 3 x 3 −1 e−(3 x ) 3 , x > 0, 3 , 3 > 0, where j = j , j , j = 1, 2, 3. The first set consists of the number of successive failures for the air conditioning system of each member in a fleet of 13 Boeing 720 jet airplanes. The pooled data, yielding a total of 213 observations, were first analyzed by Proschan (1963) and further discussed in Dahiya and Gurland (1972), Gleser (1989) and Adamidis and Loukas (1998). In the second set, the data are 109 observations on the period between successive coal-mining disasters and can be found in Cox and Lewis (1978, p. 4) and further discussed in Adamidis and Loukas (1998). Now, we focus on the analysis of earthquakes in the last century in North Anatolia fault zone between 39.00◦ –42.00◦ North latitude and 30.00◦ –40.00◦ East longitude. In Table 3, the dates of the successive earthquakes with magnitudes greater than or equal to 6 Mw (moment magnitude), which are recorded with their exact locations, magnitudes and depths between the years 1900 and 2000. Also, their exact locations are shown on a partial map in Fig. 4. The third
4504
C. Ku¸s / Computational Statistics & Data Analysis 51 (2007) 4497 – 4509
Table 2 The means and standard errors of the EM estimators for = (, ) and iterations to convergence with initial value (0) = (0) , (0) from 10 000 samples of size 100 and 500 generated from the EP se av(t) se(t) av n (0) 100
500
(0.5, 1) (1, 0.5) (0.5, 2) (2, 0.5) (2, 2)
(0.5, 1) (1, 0.5) (0.5, 2) (2, 0.5) (2, 2)
(0.906, 0.930) (1.271, 0.488) (0.930, 1.854) (2.056, 0.544) (2.083, 2.162)
(0.010, 0.002) (0.012, 0.001) (0.011, 0.004) (0.014, 0.002) (0.015, 0.008)
626.980 646.143 622.885 732.388 740.541
6.087 7.785 6.332 8.135 8.463
(0.5, 1) (1, 0.5) (0.5, 2) (2, 0.5) (2, 2)
(0.1, 0.1) (0.1, 0.1) (0.1, 0.1) (0.1, 0.1) (0.1, 0.1)
(0.906, 0.930) (1.282, 0.490) (0.916, 1.860) (2.057, 0.545) (2.032, 2.190)
(0.010, 0.002) (0.012, 0.001) (0.011, 0.004) (0.014, 0.002) (0.014, 0.008)
645.543 691.886 653.821 784.597 786.981
6.412 7.560 6.175 8.094 8.779
(0.5, 1) (1, 0.5) (0.5, 2) (2, 0.5) (2, 2)
(0.5, 1) (1, 0.5) (0.5, 2) (2, 0.5) (2, 2)
(0.546, 0.994) (1.018, 0.503) (0.540, 1.988) (2.124, 0.504) (2.137, 2.007)
(0.003, 0.001) (0.005, 0.001) (0.003, 0.002) (0.010, 0.001) (0.010, 0.005)
497.997 443.625 501.344 837.308 857.118
2.720 3.554 2.712 10.862 11.224
(0.5, 1) (1, 0.5) (0.5, 2) (2, 0.5) (2, 2)
(0.1, 0.1) (0.1, 0.1) (0.1, 0.1) (0.1, 0.1) (0.1, 0.1)
(0.539, 0.993) (1.016, 0.503) (0.539, 1.987) (2.137, 0.503) (2.146, 2.003)
(0.003, 0.001) (0.005, 0.001) (0.003, 0.002) (0.010, 0.001) (0.014, 0.007)
543.856 543.028 538.514 992.273 1009.841
1.977 3.176 2.201 11.617 17.176
Table 3 Earthquakes in North Anatolia fault zones Date
Longitude
Latitude
Magnitude (Mw)
Depth (km)
04.12.1905 09.02.1909 25.06.1910 24.01.1916 18.05.1929 19.04.1938 26.12.1939 30.07.1940 20.12.1940 08.11.1941 11.12.1942 20.12.1942 20.06.1943 26.11.1943 01.02.1944 26.10.1945 13.08.1951 07.09.1953 20.02.1956 26.05.1957 22.07.1967 03.09.1968 13.03.1992 08.03.1997 12.11.1999
39 38 34 36.83 37.9 33.79 39.51 35.25 39.2 39.5 34.83 36.8 30.51 33.72 32.69 33.29 32.87 33.01 30.49 31 30.69 32.39 39.63 35.44 31.21
39 40 41 40.27 40.2 39.44 39.8 39.64 39.11 39.74 40.76 40.7 40.85 41.05 41.41 41.54 40.88 41.09 39.89 40.67 40.67 41.81 39.72 40.78 40.74
6.8 6.3 6.2 7.1 6.1 6.6 7.9 6.2 6 6 6.1 7 6.5 7.2 7.2 6 6.9 6.4 6.4 7.1 7.2 6.5 6.1 6 7.2
0 60 0 10 10 10 20 50 0 0 40 16 10 10 10 50 10 40 40 10 33 5 23 5 25
C. Ku¸s / Computational Statistics & Data Analysis 51 (2007) 4497 – 4509
30° 42°
4505
32°
34°
36°
38°
40°
32°
34°
36°
38°
40°
756 896 461 3709
409 8592 1821 979
40°
30°
Locations of successive earthquakes Centrums Fig. 4. The successive earthquakes occurred in North Anatolia fault zone.
Table 4 Time intervals of the successive earthquakes (To be read down the columns.) 1163 501 2039 4863
3258 616 217 143
323 398 9 182
159 67 633 2117
data set given in Table 4 includes the time intervals (in days) of the successive earthquakes mentioned above. The data and Fig. 4 are taken from University of Bosphoros, Kandilli Observatory and Earthquake Research Institute-National Earthquake Monitoring Center (KOERI-NEMC, web address: http://www.koeri.boun.edu.tr). The PP-plots are given in Figs. 5–7. As thefigures show, in all cases, four models are very similar and they provide good fits to the data sets. The MLEs = , and j = j , j , j = 1, 2, 3, Kolmogorov–Smirnov (K–S) statistics with the associated p-values under the EP, EG, Weibull and Gamma models are presented in Table 5. K–S test statistics takes the smallest values and the largest p-values for all the data sets under EP model with regard to the other models. The proposed model offers an attractive alternative to these three well-established models such as Weibull, Gamma and EG as in the analyzed real data sets. From the first data, we have J11 = 16.2796,
J22 = 2853969.0620 and by inverting the J , we have the covariance matrix of , J
−1
=
J12 = 6417.8254,
0.5412 −1.22 × 10−3
−1.22 × 10−3 3.09 × 10−6
.
Hence, approximate 0.95 confidence interval for the parameters and are given by (1.2333, 1.4309) and (7.27 × 10−3 , 7.75 × 10−3 ), respectively. From the second data, we have J11 = 5.7949,
J12 = 12009.5231,
J22 = 253888796.13
4506
C. Ku¸s / Computational Statistics & Data Analysis 51 (2007) 4497 – 4509
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
. .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Fig. 5. PP-plots for the first data set.
and by inverting the J , we have the covariance matrix of , J
−1
=
8.7642 −4.15 × 10−3
−4.15 × 10−3 2 × 10−6
.
0.95 confidence interval for the parameters and are given by (2.6590, 3.7706) and Hence, approximate 1.4 × 10−3 , 1.94 × 10−3 , respectively. From the third data, we have J11 = 1.4583,
J12 = 13312.5888, J22 = 125643909.8 and by inverting the J , we have the covariance matrix of , J
−1
=
20.9354 −2.22 × 10−3
−2.22 × 10−3 2.43 × 10−7
.
0.95 confidence interval for the parameters and are given by (0.8071, 4.4683) and Hence, approximate 1.59 × 10−4 , 5.53 × 10−4 , respectively. Acknowledgements The author would like to thank the Editor-in-Chief, the associate editor and referees for carefully reading the paper and for their great help in improving the article. The author also wishes to express his gratitude to Ms. ˙Ilknur Ku¸s for her support. This research has been funded in part by both S.Ü. BAP Office and TUBITAK.
C. Ku¸s / Computational Statistics & Data Analysis 51 (2007) 4497 – 4509
.
4507
. .
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
. .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Fig. 6. PP-plots for the second data set.
Appendix Proof of Theorem 1. Since lim→0 h (; , yobs )= ni=1 e−xi −n/2 then h (; , yobs ) > 0 when ni=1 e−xi > n/2. On the other hand, since lim→∞ h (; , yobs ) = ni=1 e−xi − n then h (; , yobs ) < 0 and, therefore, there is at least one root of h (; , yobs ) = 0. For the proof of theorem, we need to show that the function h is strictly decreasing in . To prove that the h is strictly decreasing in , the first derivative of h is considered and it is given by h (; , yobs ) = −n
exp(−2) − 2 exp(−) − 2 exp(−) + 1 2 (1 − exp(−))2
.
If z() = exp(−2) − 2 exp(−) − 2 exp(−) + 1 > 0, then h (; , yobs ) < 0 and, hence, h is strictly decreasing in . Since z() = exp(−2) − 2 exp(−) − 2 exp(−) + 1 = 4 exp(−)[sinh(/2) + /2][sinh(/2) − /2] and sinh(/2) − /2 > 0, then z() > 0. This completes the proof.
Proof of Theorem 2. Let w (; , yobs ) = ni=1 xi e−xi , then it is clear that w is strictly decreasing in and lim→∞ w (; , yobs ) = 0. It follows that g (; , yobs ) < n−1 −
n
i=1
xi − lim w (; , yobs ) = n−1 − →∞
n
i=1
xi .
4508
C. Ku¸s / Computational Statistics & Data Analysis 51 (2007) 4497 – 4509
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Fig. 7. PP-plots for the third data set.
Table 5 Parameter estimates, Kolmogorov–Smirnov statistics and p-values obtained from the fit of each of the four distribution for the data set 1,2 and 3 Data set 1 (n = 213)
Distribution EP EG Weibull Gamma
2 (n = 109)
EP EG Weibull Gamma
3 (n = 24)
EP EG Weibull Gamma
Estimates = 1.3321, 7.51 × 10−3 1 = 8.01 × 10−3 , 0.4276 2 = 1.12 × 10−2 , 0.9240 3 = 9.9 × 10−3 , 0.9220
= 1 = 2.65 × 10−3 , 0.5602 2 = 4.5 × 10−3 , 0.8650 3 = 3.44 × 10−3 , 0.8280 3.2148, 1.67 × 10−3
= 2.6377, 3.56 × 10−4 1 = 6.995 × 10−4 , 1.154 × 10−5 2 = 8.12 × 10−4 , 0.7854 3 = 4.98 × 10−4 , 0.7117
K–S
p-value
0.0470
0.7351
0.0494
0.6759
0.0509
0.6393
0.0634
0.3586
0.0625
0.7876
0.0761
0.5524
0.0773
0.5325
0.0852
0.4076
0.0972
0.9772
0.1839
0.3914
0.1004
0.9690
0.1239
0.8551
C. Ku¸s / Computational Statistics & Data Analysis 51 (2007) 4497 – 4509
and, hence, g (; , yobs ) < 0 when > x −1 . On the other hand, lim→0 w (; , yobs ) = g (; , yobs ) > n−1 −
n
xi − lim w (; , yobs ) = n−1 − ( + 1) →0
i=1
n
4509
n
i=1 xi
so that
xi .
i=1
−1 and, therefore, there is at least one root of g (; , yobs ) = 0 in the < x( + 1) Hence, g(; , yobs ) > 0 when −1 −1 ; this proof is analogous to that of Theorem 4.1 in Adamidis and Loukas (1998). interval x( + 1) , x To prove that the root is unique, we consider the first derivative of g, g (; , yobs ), given by
g (; , yobs ) = −n
−2
+
n
xi2 e−xi .
i=1
The function
g
is not always monotonic for all > 0. However, if ∗ is a solution of g (; , yobs ) = 0, then
n∗−1 = ∗
n
xi2 e−xi
i=1
and n
∗ ∗ g ∗ ; , yobs = xi e− xi ∗ xi − e xi / − 1 . i=1
For simplicity, denote ∗ xi − e
∗
xi / − 1
by m(p) = p − ep / − 1 with p = ∗ xi . Then it is not difficult to see that
m(p) log() − 2,
> 0. Therefore, m(p) < 0 for ∈ 0, e2 . Hence, g ∗ ; , yobs < 0 for all ∗ , xi > 0, i = 1, 2, . . . , n and ∈ 0, e2 . 2 This nimplies that g is negative at its stationary points for ∈ 0, e . Taking into account the lim→∞ g (; , yobs ) = − i=1 xi , the uniqueness of the root in the specified interval is proven. References Adamidis, K., 1999. An EM algorithm for estimating negative binomial parameters. Austral. New Zealand J. Statist. 41 (2), 213–221. Adamidis, K., Loukas, S., 1998. A life time distribution with decreasing failure rate. Statist. Probab. Lett. 39, 35–42. Adamidis, K., Dimitrakopoulou, T., Loukas, S., 2005. On an extension of the exponential–geometric distribution. Statist. Probab. Lett. 73, 259–269. Barlow, R.E., Marshall, A.W., 1964. Bounds for distributions with monotone hazard rate I and II. Ann. Math. Statist. 35, 1234–1274. Barlow, R.E., Marshall, A.W., 1965. Tables of bounds for distributions with monotone hazard rate. J. Amer. Statist. Assoc. 60, 872–890. Barlow, R.E., Marshall, A.W., Proschan, F., 1963. Properties of probability distributions with monotone hazard rate. Ann. Math. Statist. 34, 375–389. Cox, D.R., Lewis, P.A.W., 1978. The Statistical Analysis of Series of Events. Chapman & Hall, London. Cozzolino, J.M., 1968. Probabilistic models of decreasing failure rate processes. Naval Res. Logist. Quart. 15, 361–374. Dahiya, R.C., Gurland, J., 1972. Goodness of fit tests for the gamma and exponential distributions. Technometrics 14, 791–801. Dempster, A.P., Laird, N.M., Rubin, D.B., 1977. Maximum likelihood from incomplete data via the EM algorithm (with discussion). J. Roy. Statist. Soc. Ser. B 39, 1–38. Gleser, L.J., 1989. The gamma distribution as a mixture of exponential distributions. Amer. Statist. 43, 115–117. Gurland, J., Sethuraman, J., 1994. Reversal of increasing failure rates when pooling failure data. Technometrics 36, 416–418. Karlis, D., 2003. An EM algorithm for multivariate Poisson distribution and related models. J. Appl. Statist. 30 1, 63–77. Little, R.J.A., Rubin, D.B., 1983. Incomplete data. In: Kotz, S., Johnson, N.L. (Eds.), Encyclopedia of Statistical Sciences, vol. 4. Wiley, New York. Lomax, K.S., 1954. Business failures: another example of the analysis of failure data. J. Amer. Statist. Assoc. 49, 847–852. Marshall, A.W., Proschan, F., 1965. Maximum likelihood estimates for distributions with monotone failure rate. Ann. Math. Statist. 36, 69–77. McLachlan, G.J., Krishnan, T., 1997. The EM Algorithm and Extension. Wiley, New York. McNolty, F., Doyle, J., Hansen, E., 1980. Properties of the mixed exponential failure process. Technometrics 22, 555–565. Nassar, M.M., 1988. Two properties of mixtures of exponential distributions. IEEE Trans. Raliab. 37 (4), 383–385. Ng, H.K.T., Chan, P.S., Balakrishnan, N., 2002. Estimation of parameters from progressively censored data using EM algorithm. Comput. Statist. Data Anal. 39, 371–386. Proschan, F., 1963. Theoretical explanation of observed decreasing failure rate. Technometrics 5, 375–383. Saunders, S.C., Myhre, J.M., 1983. Maximum likelihood estimation for two-parameter decreasing hazard rate distributions using censored data. J. Amer. Statist. Assoc. 78, 664–673.
A new lifetime distribution Co¸skun Ku¸s∗ Department of Statistics, Faculty of Art and Science, Selcuk University, 42031, Konya, Turkey Received 10 November 2005; received in revised form 26 June 2006; accepted 12 July 2006 Available online 4 August 2006
Abstract A new two-parameter distribution with decreasing failure rate is introduced. Various properties of the introduced distribution are discussed. The EM algorithm is used to determine the maximum likelihood estimates and the asymptotic variances and covariance of these estimates are obtained. Simulation studies are performed in order to assess the accuracy of the approximation of the variances and covariance of the maximum likelihood estimates and investigate the convergence of the proposed EM scheme. Illustrative examples based on real data are also given. © 2006 Elsevier B.V. All rights reserved. Keywords: Compounding; Decreasing failure rate; EM algorithm; Exponential distribution; Lifetime distributions; Maximum likelihood estimation; Zero truncated Poisson distribution
1. Introduction In general, a population is expected to exhibit decreasing failure rate (DFR) when its behavior over time is characterized by ‘work-hardening’ (in engineering terms) or ‘immunity’ (in biological terms); sometimes the broader term ‘infant mortality’ is used to denote the DFR phenomenon (Adamidis and Loukas, 1998). The distributions with DFR are discussed in the works of Lomax (1954), Proschan (1963), Barlow et al. (1963), Barlow and Marshall (1964, 1965), Marshall and Proschan (1965), Cozzolino (1968), Dahiya and Gurland (1972), McNolty et al. (1980), Saunders and Myhre (1983), Nassar (1988), Gleser (1989), Gurland and Sethuraman (1994) and Adamidis and Loukas (1998). This paper is organized as follows: In Section 2, a new DFR distribution is obtained by mixing exponential and zero truncated Poisson distribution where mixing procedure was previously carried out by Adamidis and Loukas (1998). In Section 3, various properties of the introduced distribution are discussed. In Section 4, the estimation of parameters is studied by the method of maximum likelihood. Finally, in Section 5, illustrative examples based on real data are provided to close the paper. 2. The distribution Let W1 , W2 , . . . , WZ be a random sample from distribution with density f (w; ) = e−w , w, ∈ R+ , Z is a −1 , z ∈ N, ∈ R+ zero truncated Poisson variable with probability function P (z; ) = e− z −1 (z + 1) 1 − e− ∗ Tel.: +90 3322231354; fax: +90 3322410106.
E-mail address: [email protected]. 0167-9473/$ - see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.csda.2006.07.017
4498
C. Ku¸s / Computational Statistics & Data Analysis 51 (2007) 4497 – 4509
.
.
.
.
. .
.
.
.
Fig. 1. Probability density functions of the EP distribution for = 1 and = 1, 2, 3, identified, respectively, by the y-axis intercepts appearing in increasing order of magnitude.
.
.
.
.
.
.
.
.
Fig. 2. Probability density functions of the EP distribution for = 1 and = 1, 2, 3, identified, respectively, by the y-axis intercepts appearing in increasing order of magnitude.
where (·) is the Gamma function and Z and W s are independent. Let us define X = min (W1 , W2 , . . . , WZ ). Then, f (x|z; ) = ze−zx and marginal probability density function of X is e−−x+ exp(−x) , 1 − e−
f (x; ) =
x, , ∈ R+ ,
(1)
where = (, ). In the sequel, distribution of X will be referred to as the exponential–Poisson distribution (EP) which is customary for such names to be given to the distributions arising via the operation of compounding (mixing) in −1 the literature. It can be seen that the EP density function is monotone decreasing with modal value 1 − e− at x = 0. EP probability density functions are displayed in Figs. 1 and 2 for selected parameter values. For all values of parameters, the density is strictly decreasing in x and tending to zero as x → ∞. Its graph resembles those of the EG, Pareto II and exponential distributions. As approaches zero, the EP leads to exponential distribution with parameter . The model is obtained under the concept of population heterogeneity (through the process of compounding). An interpretation of the proposed model is as follows: a situation where failure (of a device for example) occurs due to the presence of an unknown number, Z, of initial defects of same kind (a number of semiconductors from a defective lot, for example). The Ws represent their lifetimes and each defect can be detected only after causing failure, in which case
C. Ku¸s / Computational Statistics & Data Analysis 51 (2007) 4497 – 4509
4499
it is repaired perfectly (Adamidis and Loukas, 1998). Then the distributional assumptions given earlier lead to the EP distribution for modelling the time to the first failure X. 3. Properties of the distribution 3.1. Distribution function and moments The distribution function is given by −1 F (x; ) = e exp(−x) − e 1 − e
(2)
−1 −1 . and hence, the median is obtained by log log 2−1 e + 1 Since the EP is not part of the exponential family, the evaluations of moments become cumbersome. The raw moments of X are determined from (1) by direct integration. For r ∈ N, raw moments are given by (r + 1) Fr+1,r+1 ([1, 1, . . . , 1], [2, 2, . . . , 2], ), E Xr ; = e − 1 r
r ∈ N,
where Fp,q (n, d, ) is the generalized hypergeometric function. This function is also known as Barnes’s extended hypergeometric function. The definition of Fp,q (n, d, ) is Fp,q (n, d, ) =
∞
k=0
(k
p
−1 (ni ) i=1 (ni + k) , q + 1) i=1 (di + k) −1 (di )
k
where n = [n1 , n2 , . . . , np ], p is the number of operands of n, d = [d1 , d2 , . . . , dq ] and q is the number of operands of d. Generalized hypergeometric function is quickly evaluated and readily available in standard software such as Maple. Hence the mean and variance of the EP distribution are given, respectively, by F2,2 ([1, 1], [2, 2], ), E(X; ) = e −1
2 F ([1, 1], [2, 2], ) . 2F3,3 ([1, 1, 1], [2, 2, 2], ) − Var(X; ) = e − 1 2,2 e − 1 2
(3)
3.2. The survivor, hazard and mean residual life functions Using (1) and (2), survival function (also known reliability function) and hazard function (also known failure rate function) of the EP distribution are given, respectively, by −1 s(x; ) = 1 − F (x; ) = 1 − e exp(−x) 1 − e , −−x+ exp(−x) 1 − e f (x; ) e . h(x; ) = = s(x; ) 1 − e− 1 − e exp(−x) Both functions have simple forms in contrast to those of some DFR distributions such as Gamma with DFR. The hazard function is decreasing because the DFR property follows from the results of Barlow et al. (1963) on mixtures. The EP hazard functions for selected parameter values are displayed in Fig. 3. −1 and h(∞; ) = , respectively. They are both finite The initial and long-term hazards are h(0; ) = e e − 1 in contrast to those of Weibull distribution with h(0; ) = ∞ and h(∞; ) = 0.
4500
C. Ku¸s / Computational Statistics & Data Analysis 51 (2007) 4497 – 4509
.
.
.
.
.
.
.
Fig. 3. Hazard functions of the EP distribution for = 1 and = 0.5, 1.0, 1.5, identified, respectively, by the y-axis intercepts appearing in increasing order of magnitude.
The mean residual lifetime of the EP distribution is given by m (x0 ; ) = E (X − x0 |X x0 ; ) =
e−x0 F2,2 [1, 1], [2, 2], e−x0 , e exp(−x0 ) − 1
where Fp,q (n, d, ) is the generalized hypergeometric function. Because the EP is a DFR distribution m (x0 ; ) increases monotonically up to its limit at infinity. Clearly, m (x0 ; ) → E(X; ) as x0 → 0 and m (x0 ; ) → −1 as x0 → ∞. 4. Estimation of the parameters 4.1. Estimation by maximum likelihood The log-likelihood function based on the observed sample size of n, yobs = (xi ; i = 1, 2, . . . , n), from the EP distribution is given by (; yobs ) = n log() −
n
xi − n +
i=1
n
e−xi − n log 1 − e−
i=1
and subsequently the associated gradients are found to be n −1 j (; yobs ) −xi , = e + n −1 − 1 − e− j
(4)
i=1
j (; yobs ) = n−1 − xi − xi e−xi . j n
n
i=1
i=1
The equation j (; yobs ) /j = 0 could be solved exactly for , namely
n −1 1
−1 Xi exp −Xi = −X n
(5)
(6)
i=1
conditional upon the value of , where and are the maximum likelihood estimates (MLEs) for the parameters and , respectively. Eq. (6) might be used for building a profile likelihood. In the following, Theorems 1 and 2 express when the other parameter is given (or known), what the conditions are there to obtain the existence and uniqueness of the MLE.
C. Ku¸s / Computational Statistics & Data Analysis 51 (2007) 4497 – 4509
4501
in (4), where is the true value of the Theorem 1. Let h (; , yobs ) denote the function on the RHS of the expression , is unique for n/2 < ni=1 e−xi . parameter. Then, the root of h (; , yobs ) = 0, Theorem 2. Let g (; , yobs ) denote the function on the RHS of the expression in (5) and x = n−1 ni=1 xi , where 2 , lies in the interval is the true value of the parameter. Then, for a given ∈ 0, e the root of g (; , yobs ) = 0, −1 −1 x( + 1) , x and is unique. 4.2. An EM algorithm The MLEs of and must be derived numerically. Newton–Raphson algorithm is one of the standard methods to determine the MLEs of the parameters. To employ the algorithm, second derivatives of the log-likelihood are required for all iteration. EM algorithm is a very powerful tool in handling the incomplete data problem (Dempster et al., 1977; McLachlan and Krishnan, 1997). It is an iterative method by repeatedly replacing the missing data with estimated values and updating the parameter estimates. It is especially useful if the complete data set is easy to analyze. As pointed out by Little and Rubin (1983), the EM algorithm will converge reliably but rather slowly (as compared to the Newton–Raphson method) when the amount of information in the missing data is relatively large. Recently, EM algorithm has been used by several authors such as Adamidis and Loukas (1998), Adamidis (1999), Ng et al. (2002), Karlis (2003) and Adamidis et al. (2005). To start the algorithm, hypothetical complete-data distribution is defined with density function f (x, z; ) = −1 , x > 0, z = 1, 2, . . . , , > 0. Thus, it is straightforward to verify that the Eze−zx e− z −1 (z + 1) 1 − e− step of an EM cycle requires the computation of the conditional expectation of Z|X; (h) , where (h) = (h) , (h) is the current estimate of . Using that P (z|x; ) = e−xz− exp(−x)+x z−1 −1 (z), z ∈ N, this is found to be E(Z|X; ) = 1 + e−x . The EM cycle is completed with M-step, which is complete data maximum likelihood over , with the missing Z’s replaced by their conditional expectations E(Z|X; ) (Adamidis and Loukas, 1998). Thus, an EM iteration, taking (h) into (h+1) , is given by
(h+1)
=n
n
−1
xi 1 +
(h) −(h) xi
e
,
i=1
(h+1)
=n
−1
1−e
−(h+1)
n
1 + (h) e−
(h)
xi
.
(7)
i=1
It can be seen that only a one-dimensional search such as Newton–Raphson is required for M-step of an EM cycle. Ng et al. (2002) used the same method for estimation of the parameters of the Weibull distribution based on progressively type-II right censored sample. 4.3. Asymptotic variances and covariance of the MLEs Applying the usual large sample approximation, the MLE of can be treated as being approximately bivariate normal with mean and variance–covariance matrix, which is the inverse of the expected information matrix J () = E(I ; ), where I = I (; yobs ) is the observed information matrix with elements Iij = −j2 /ji jj with i, j = 1, 2 and the expectation is to be taken with respect to the distribution of X. Differentiating (4) and (5), the elements of the symmetric, second-order observed information matrix are found to be −2 , I11 = n 1 + e2 − 2 e − 2e 1 − e− I12 =
n
i=1
xi e
−xi
,
I22 = n
−2
−
n
i=1
xi2 e−xi .
4502
C. Ku¸s / Computational Statistics & Data Analysis 51 (2007) 4497 – 4509
The elements of the expected information matrix, J (), are calculated by taking the expectations of Iij , i, j = 1, 2 with respect to the distribution of X, i.e. following expectations are required: E Xe−X =
e− F2,2 ([2, 2], [3, 3], ), 4 1 − e−
E X2 e−X =
e− F3,3 ([2, 2, 2], [3, 3, 3], ). 4 1 − e− 2
Thus J () are derived in the form −2 , J11 = n 1 + e2 − 2 e − 2e 1 − e− J12 =
ne− F2,2 ([2, 2], [3, 3], ), 4 1 − e−
J22 = n−2 −
n2 e− F3,3 ([2, 2, 2], [3, 3, 3], ). 42 1 − e−
The inverse of J (), evaluated at provides the asymptotic variance–covariance matrix of the MLEs. Alternative estimates can be obtained from the inverse of the observed information matrix since it is a consistent estimator of J −1 ().
4.4. Simulation study Simulation results are obtained based on the assumption of Theorem 1, to guarantee the existence and uniqueness of the solution in each iteration of (h+1) in EM algorithm. No restriction has been imposed on the maximum number of iterations and convergence is assumed when the absolute differences between successive estimates are less than 10−5 . For each value of , the parameter estimates have been found by employing (7) with different initial values. In order to assess the accuracy of the approximation of the variances and covariance of the MLEs determined from the information a simulation matrix, study (based on 10 000 simulations) has been carried out. The simulated values of Var , Var and Cov , as well as the approximate values determined by averaging the corresponding values obtained from the expected and observed information matrix are presented in Table 1. From Table 1, it is observed that the approximate values determined from expected and observed information matrix are quite close to the simulated values for large values of n. Furthermore, it is noted that the approximation becomes quite accurate as n increases. In addition, variances and covariances of the MLEs obtained from the observed information matrix is quite close to the that of expected information matrix for large values of n. Also simulations have been performed to investigate the convergence of the proposed EM scheme. Ten thousand samples of size 100 and 500 each of which are randomly sampled from the EP for each of five values of are generated. The results from simulated data sets are reported in Table 2, which gives the averages of the 10 000 MLEs av and averages numbers of iterations to convergence av(t) together with standard errors se and se(t). Table 2 indicates the following results: convergence has been achieved in all cases, even when the starting values are poor and this emphasizes the numerical stability of the EM algorithm. The values of av and se suggest the EM estimates performed consistently. Standard errors of the MLEs decrease when sample size increases. MLEs of the parameters and have bias. However, bias of the MLEs are reduced for large sample size in five cases of parameters in Table 2.
C. Ku¸s / Computational Statistics & Data Analysis 51 (2007) 4497 – 4509
4503
Table 1 Variances and covariance of the MLEs
n
Simulated Var
Var
, Cov
From expected information , Var Var Cov
From observed information , Var Var Cov
50
(0.5, 1) (1, 0.5) (0.5, 2) (2, 0.5) (2, 2)
1.5599 1.8175 1.6120 2.2499 2.2971
0.0662 0.0242 0.2688 0.0502 0.7984
−0.2564 −0.1715 −0.5273 −0.2817 −1.1291
3.9604 5.2131 4.0185 8.5655 8.4945
0.1507 0.0545 0.5990 0.1178 1.8565
−0.6135 −0.4319 −1.2233 −0.8435 −3.3073
3.1078 3.5839 3.1327 4.5887 4.4832
0.1291 0.0422 0.5208 0.0744 1.1730
−0.4909 −0.3052 −0.9892 −0.4571 −1.7868
100
(0.5, 1) (1, 0.5) (0.5, 2) (2, 0.5) (2, 2)
1.0886 1.4367 1.1664 2.0511 2.1500
0.0429 0.0173 0.1747 0.0385 0.6070
−0.1781 −0.1352 −0.3744 −0.2477 −1.0097
1.3737 2.0329 1.4191 4.2078 4.2756
0.0661 0.0245 0.2655 0.0566 0.8995
−0.2389 −0.1811 −0.4825 −0.4132 −1.6525
1.3121 1.8245 1.2934 2.7439 2.8582
0.0651 0.0224 0.2579 0.0404 0.6671
−0.2344 −0.1605 −0.4625 −0.2685 −1.1210
500
(0.5, 1) (1, 0.5) (0.5, 2) (2, 0.5) (2, 2)
0.1126 0.2821 0.1119 0.9741 0.9941
0.0098 0.0049 0.0391 0.0139 0.2241
−0.0289 −0.0330 −0.0574 −0.1080 −0.4383
0.1333 0.2063 0.1328 0.7133 0.7222
0.0118 0.0045 0.0471 0.0116 0.1851
−0.0355 −0.0273 −0.0708 −0.0811 −0.3263
0.1325 0.2145 0.1323 0.7444 0.7630
0.0118 0.0045 0.0469 0.0114 0.1825
−0.0354 −0.0280 −0.0706 −0.0082 −0.3292
1000
(0.5, 1) (1, 0.5) (0.5, 2) (2, 0.5) (2, 2)
0.0570 0.0949 0.0575 0.5030 0.4862
0.0054 0.0022 0.0215 0.0074 0.1176
−0.0155 −0.0135 −0.0311 −0.0573 −0.2239
0.0631 0.0919 0.0631 0.3081 0.3055
0.0058 0.0022 0.0231 0.0059 0.0941
−0.0171 −0.0131 −0.0342 −0.0393 −0.1570
0.0633 0.0918 0.0636 0.3437 0.3528
0.0058 0.0022 0.0232 0.0061 0.0977
−0.0171 −0.0131 −0.0344 −0.0422 −0.1700
10 000
(0.5, 1) (1, 0.5) (0.5, 2) (2, 0.5) (2, 2)
0.0061 0.0086 0.0062 0.0236 0.0233
0.0006 0.0002 0.0023 0.0006 0.0092
−0.0017 −0.0012 −0.0034 −0.0035 −0.0142
0.0062 0.0087 0.0062 0.0232 0.0237
0.0006 0.0002 0.0023 0.0006 0.0093
−0.0017 −0.0013 −0.0034 −0.0036 −0.0143
0.0062 0.0087 0.0062 0.0237 0.0237
0.0006 0.0002 0.0023 0.0006 0.0093
−0.0017 −0.0013 −0.0034 −0.0036 −0.0143
5. Illustrative examples The fit of the EP distribution of real data is examined by graphical methods using MLEs. It is also compared with the EG, Weibull and Gamma models with respective densities −2 f1 (x; 1 ) = 1 1 − 1 e−1 x 1 − 1 e−1 x , x > 0, 1 > 0, 1 ∈ (0, 1),
2
f2 (x; 2 ) = 2 2 2 x 2 −1 e−(2 x ) ,
x > 0, 2 , 2 > 0,
−1 f3 (x; 3 ) = 3 3 x 3 −1 e−(3 x ) 3 , x > 0, 3 , 3 > 0, where j = j , j , j = 1, 2, 3. The first set consists of the number of successive failures for the air conditioning system of each member in a fleet of 13 Boeing 720 jet airplanes. The pooled data, yielding a total of 213 observations, were first analyzed by Proschan (1963) and further discussed in Dahiya and Gurland (1972), Gleser (1989) and Adamidis and Loukas (1998). In the second set, the data are 109 observations on the period between successive coal-mining disasters and can be found in Cox and Lewis (1978, p. 4) and further discussed in Adamidis and Loukas (1998). Now, we focus on the analysis of earthquakes in the last century in North Anatolia fault zone between 39.00◦ –42.00◦ North latitude and 30.00◦ –40.00◦ East longitude. In Table 3, the dates of the successive earthquakes with magnitudes greater than or equal to 6 Mw (moment magnitude), which are recorded with their exact locations, magnitudes and depths between the years 1900 and 2000. Also, their exact locations are shown on a partial map in Fig. 4. The third
4504
C. Ku¸s / Computational Statistics & Data Analysis 51 (2007) 4497 – 4509
Table 2 The means and standard errors of the EM estimators for = (, ) and iterations to convergence with initial value (0) = (0) , (0) from 10 000 samples of size 100 and 500 generated from the EP se av(t) se(t) av n (0) 100
500
(0.5, 1) (1, 0.5) (0.5, 2) (2, 0.5) (2, 2)
(0.5, 1) (1, 0.5) (0.5, 2) (2, 0.5) (2, 2)
(0.906, 0.930) (1.271, 0.488) (0.930, 1.854) (2.056, 0.544) (2.083, 2.162)
(0.010, 0.002) (0.012, 0.001) (0.011, 0.004) (0.014, 0.002) (0.015, 0.008)
626.980 646.143 622.885 732.388 740.541
6.087 7.785 6.332 8.135 8.463
(0.5, 1) (1, 0.5) (0.5, 2) (2, 0.5) (2, 2)
(0.1, 0.1) (0.1, 0.1) (0.1, 0.1) (0.1, 0.1) (0.1, 0.1)
(0.906, 0.930) (1.282, 0.490) (0.916, 1.860) (2.057, 0.545) (2.032, 2.190)
(0.010, 0.002) (0.012, 0.001) (0.011, 0.004) (0.014, 0.002) (0.014, 0.008)
645.543 691.886 653.821 784.597 786.981
6.412 7.560 6.175 8.094 8.779
(0.5, 1) (1, 0.5) (0.5, 2) (2, 0.5) (2, 2)
(0.5, 1) (1, 0.5) (0.5, 2) (2, 0.5) (2, 2)
(0.546, 0.994) (1.018, 0.503) (0.540, 1.988) (2.124, 0.504) (2.137, 2.007)
(0.003, 0.001) (0.005, 0.001) (0.003, 0.002) (0.010, 0.001) (0.010, 0.005)
497.997 443.625 501.344 837.308 857.118
2.720 3.554 2.712 10.862 11.224
(0.5, 1) (1, 0.5) (0.5, 2) (2, 0.5) (2, 2)
(0.1, 0.1) (0.1, 0.1) (0.1, 0.1) (0.1, 0.1) (0.1, 0.1)
(0.539, 0.993) (1.016, 0.503) (0.539, 1.987) (2.137, 0.503) (2.146, 2.003)
(0.003, 0.001) (0.005, 0.001) (0.003, 0.002) (0.010, 0.001) (0.014, 0.007)
543.856 543.028 538.514 992.273 1009.841
1.977 3.176 2.201 11.617 17.176
Table 3 Earthquakes in North Anatolia fault zones Date
Longitude
Latitude
Magnitude (Mw)
Depth (km)
04.12.1905 09.02.1909 25.06.1910 24.01.1916 18.05.1929 19.04.1938 26.12.1939 30.07.1940 20.12.1940 08.11.1941 11.12.1942 20.12.1942 20.06.1943 26.11.1943 01.02.1944 26.10.1945 13.08.1951 07.09.1953 20.02.1956 26.05.1957 22.07.1967 03.09.1968 13.03.1992 08.03.1997 12.11.1999
39 38 34 36.83 37.9 33.79 39.51 35.25 39.2 39.5 34.83 36.8 30.51 33.72 32.69 33.29 32.87 33.01 30.49 31 30.69 32.39 39.63 35.44 31.21
39 40 41 40.27 40.2 39.44 39.8 39.64 39.11 39.74 40.76 40.7 40.85 41.05 41.41 41.54 40.88 41.09 39.89 40.67 40.67 41.81 39.72 40.78 40.74
6.8 6.3 6.2 7.1 6.1 6.6 7.9 6.2 6 6 6.1 7 6.5 7.2 7.2 6 6.9 6.4 6.4 7.1 7.2 6.5 6.1 6 7.2
0 60 0 10 10 10 20 50 0 0 40 16 10 10 10 50 10 40 40 10 33 5 23 5 25
C. Ku¸s / Computational Statistics & Data Analysis 51 (2007) 4497 – 4509
30° 42°
4505
32°
34°
36°
38°
40°
32°
34°
36°
38°
40°
756 896 461 3709
409 8592 1821 979
40°
30°
Locations of successive earthquakes Centrums Fig. 4. The successive earthquakes occurred in North Anatolia fault zone.
Table 4 Time intervals of the successive earthquakes (To be read down the columns.) 1163 501 2039 4863
3258 616 217 143
323 398 9 182
159 67 633 2117
data set given in Table 4 includes the time intervals (in days) of the successive earthquakes mentioned above. The data and Fig. 4 are taken from University of Bosphoros, Kandilli Observatory and Earthquake Research Institute-National Earthquake Monitoring Center (KOERI-NEMC, web address: http://www.koeri.boun.edu.tr). The PP-plots are given in Figs. 5–7. As thefigures show, in all cases, four models are very similar and they provide good fits to the data sets. The MLEs = , and j = j , j , j = 1, 2, 3, Kolmogorov–Smirnov (K–S) statistics with the associated p-values under the EP, EG, Weibull and Gamma models are presented in Table 5. K–S test statistics takes the smallest values and the largest p-values for all the data sets under EP model with regard to the other models. The proposed model offers an attractive alternative to these three well-established models such as Weibull, Gamma and EG as in the analyzed real data sets. From the first data, we have J11 = 16.2796,
J22 = 2853969.0620 and by inverting the J , we have the covariance matrix of , J
−1
=
J12 = 6417.8254,
0.5412 −1.22 × 10−3
−1.22 × 10−3 3.09 × 10−6
.
Hence, approximate 0.95 confidence interval for the parameters and are given by (1.2333, 1.4309) and (7.27 × 10−3 , 7.75 × 10−3 ), respectively. From the second data, we have J11 = 5.7949,
J12 = 12009.5231,
J22 = 253888796.13
4506
C. Ku¸s / Computational Statistics & Data Analysis 51 (2007) 4497 – 4509
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
. .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Fig. 5. PP-plots for the first data set.
and by inverting the J , we have the covariance matrix of , J
−1
=
8.7642 −4.15 × 10−3
−4.15 × 10−3 2 × 10−6
.
0.95 confidence interval for the parameters and are given by (2.6590, 3.7706) and Hence, approximate 1.4 × 10−3 , 1.94 × 10−3 , respectively. From the third data, we have J11 = 1.4583,
J12 = 13312.5888, J22 = 125643909.8 and by inverting the J , we have the covariance matrix of , J
−1
=
20.9354 −2.22 × 10−3
−2.22 × 10−3 2.43 × 10−7
.
0.95 confidence interval for the parameters and are given by (0.8071, 4.4683) and Hence, approximate 1.59 × 10−4 , 5.53 × 10−4 , respectively. Acknowledgements The author would like to thank the Editor-in-Chief, the associate editor and referees for carefully reading the paper and for their great help in improving the article. The author also wishes to express his gratitude to Ms. ˙Ilknur Ku¸s for her support. This research has been funded in part by both S.Ü. BAP Office and TUBITAK.
C. Ku¸s / Computational Statistics & Data Analysis 51 (2007) 4497 – 4509
.
4507
. .
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
. .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Fig. 6. PP-plots for the second data set.
Appendix Proof of Theorem 1. Since lim→0 h (; , yobs )= ni=1 e−xi −n/2 then h (; , yobs ) > 0 when ni=1 e−xi > n/2. On the other hand, since lim→∞ h (; , yobs ) = ni=1 e−xi − n then h (; , yobs ) < 0 and, therefore, there is at least one root of h (; , yobs ) = 0. For the proof of theorem, we need to show that the function h is strictly decreasing in . To prove that the h is strictly decreasing in , the first derivative of h is considered and it is given by h (; , yobs ) = −n
exp(−2) − 2 exp(−) − 2 exp(−) + 1 2 (1 − exp(−))2
.
If z() = exp(−2) − 2 exp(−) − 2 exp(−) + 1 > 0, then h (; , yobs ) < 0 and, hence, h is strictly decreasing in . Since z() = exp(−2) − 2 exp(−) − 2 exp(−) + 1 = 4 exp(−)[sinh(/2) + /2][sinh(/2) − /2] and sinh(/2) − /2 > 0, then z() > 0. This completes the proof.
Proof of Theorem 2. Let w (; , yobs ) = ni=1 xi e−xi , then it is clear that w is strictly decreasing in and lim→∞ w (; , yobs ) = 0. It follows that g (; , yobs ) < n−1 −
n
i=1
xi − lim w (; , yobs ) = n−1 − →∞
n
i=1
xi .
4508
C. Ku¸s / Computational Statistics & Data Analysis 51 (2007) 4497 – 4509
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
. .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
Fig. 7. PP-plots for the third data set.
Table 5 Parameter estimates, Kolmogorov–Smirnov statistics and p-values obtained from the fit of each of the four distribution for the data set 1,2 and 3 Data set 1 (n = 213)
Distribution EP EG Weibull Gamma
2 (n = 109)
EP EG Weibull Gamma
3 (n = 24)
EP EG Weibull Gamma
Estimates = 1.3321, 7.51 × 10−3 1 = 8.01 × 10−3 , 0.4276 2 = 1.12 × 10−2 , 0.9240 3 = 9.9 × 10−3 , 0.9220
= 1 = 2.65 × 10−3 , 0.5602 2 = 4.5 × 10−3 , 0.8650 3 = 3.44 × 10−3 , 0.8280 3.2148, 1.67 × 10−3
= 2.6377, 3.56 × 10−4 1 = 6.995 × 10−4 , 1.154 × 10−5 2 = 8.12 × 10−4 , 0.7854 3 = 4.98 × 10−4 , 0.7117
K–S
p-value
0.0470
0.7351
0.0494
0.6759
0.0509
0.6393
0.0634
0.3586
0.0625
0.7876
0.0761
0.5524
0.0773
0.5325
0.0852
0.4076
0.0972
0.9772
0.1839
0.3914
0.1004
0.9690
0.1239
0.8551
C. Ku¸s / Computational Statistics & Data Analysis 51 (2007) 4497 – 4509
and, hence, g (; , yobs ) < 0 when > x −1 . On the other hand, lim→0 w (; , yobs ) = g (; , yobs ) > n−1 −
n
xi − lim w (; , yobs ) = n−1 − ( + 1) →0
i=1
n
4509
n
i=1 xi
so that
xi .
i=1
−1 and, therefore, there is at least one root of g (; , yobs ) = 0 in the < x( + 1) Hence, g(; , yobs ) > 0 when −1 −1 ; this proof is analogous to that of Theorem 4.1 in Adamidis and Loukas (1998). interval x( + 1) , x To prove that the root is unique, we consider the first derivative of g, g (; , yobs ), given by
g (; , yobs ) = −n
−2
+
n
xi2 e−xi .
i=1
The function
g
is not always monotonic for all > 0. However, if ∗ is a solution of g (; , yobs ) = 0, then
n∗−1 = ∗
n
xi2 e−xi
i=1
and n
∗ ∗ g ∗ ; , yobs = xi e− xi ∗ xi − e xi / − 1 . i=1
For simplicity, denote ∗ xi − e
∗
xi / − 1
by m(p) = p − ep / − 1 with p = ∗ xi . Then it is not difficult to see that
m(p) log() − 2,
> 0. Therefore, m(p) < 0 for ∈ 0, e2 . Hence, g ∗ ; , yobs < 0 for all ∗ , xi > 0, i = 1, 2, . . . , n and ∈ 0, e2 . 2 This nimplies that g is negative at its stationary points for ∈ 0, e . Taking into account the lim→∞ g (; , yobs ) = − i=1 xi , the uniqueness of the root in the specified interval is proven. References Adamidis, K., 1999. An EM algorithm for estimating negative binomial parameters. Austral. New Zealand J. Statist. 41 (2), 213–221. Adamidis, K., Loukas, S., 1998. A life time distribution with decreasing failure rate. Statist. Probab. Lett. 39, 35–42. Adamidis, K., Dimitrakopoulou, T., Loukas, S., 2005. On an extension of the exponential–geometric distribution. Statist. Probab. Lett. 73, 259–269. Barlow, R.E., Marshall, A.W., 1964. Bounds for distributions with monotone hazard rate I and II. Ann. Math. Statist. 35, 1234–1274. Barlow, R.E., Marshall, A.W., 1965. Tables of bounds for distributions with monotone hazard rate. J. Amer. Statist. Assoc. 60, 872–890. Barlow, R.E., Marshall, A.W., Proschan, F., 1963. Properties of probability distributions with monotone hazard rate. Ann. Math. Statist. 34, 375–389. Cox, D.R., Lewis, P.A.W., 1978. The Statistical Analysis of Series of Events. Chapman & Hall, London. Cozzolino, J.M., 1968. Probabilistic models of decreasing failure rate processes. Naval Res. Logist. Quart. 15, 361–374. Dahiya, R.C., Gurland, J., 1972. Goodness of fit tests for the gamma and exponential distributions. Technometrics 14, 791–801. Dempster, A.P., Laird, N.M., Rubin, D.B., 1977. Maximum likelihood from incomplete data via the EM algorithm (with discussion). J. Roy. Statist. Soc. Ser. B 39, 1–38. Gleser, L.J., 1989. The gamma distribution as a mixture of exponential distributions. Amer. Statist. 43, 115–117. Gurland, J., Sethuraman, J., 1994. Reversal of increasing failure rates when pooling failure data. Technometrics 36, 416–418. Karlis, D., 2003. An EM algorithm for multivariate Poisson distribution and related models. J. Appl. Statist. 30 1, 63–77. Little, R.J.A., Rubin, D.B., 1983. Incomplete data. In: Kotz, S., Johnson, N.L. (Eds.), Encyclopedia of Statistical Sciences, vol. 4. Wiley, New York. Lomax, K.S., 1954. Business failures: another example of the analysis of failure data. J. Amer. Statist. Assoc. 49, 847–852. Marshall, A.W., Proschan, F., 1965. Maximum likelihood estimates for distributions with monotone failure rate. Ann. Math. Statist. 36, 69–77. McLachlan, G.J., Krishnan, T., 1997. The EM Algorithm and Extension. Wiley, New York. McNolty, F., Doyle, J., Hansen, E., 1980. Properties of the mixed exponential failure process. Technometrics 22, 555–565. Nassar, M.M., 1988. Two properties of mixtures of exponential distributions. IEEE Trans. Raliab. 37 (4), 383–385. Ng, H.K.T., Chan, P.S., Balakrishnan, N., 2002. Estimation of parameters from progressively censored data using EM algorithm. Comput. Statist. Data Anal. 39, 371–386. Proschan, F., 1963. Theoretical explanation of observed decreasing failure rate. Technometrics 5, 375–383. Saunders, S.C., Myhre, J.M., 1983. Maximum likelihood estimation for two-parameter decreasing hazard rate distributions using censored data. J. Amer. Statist. Assoc. 78, 664–673.
