Estimation and Prediction for Gompertz Distribution under General Progressive Censoring

In this article, we discuss the estimation of the parameters for Gompertz distribution and prediction using general progressive Type-II censoring. Based on the Expectation–Maximization algorithm, we calculate the maximum likelihood estimates. Bayesian estimates are considered under different loss functions, which are symmetrical, asymmetrical and balanced, respectively. An approximate method—Tierney and Kadane—is used to derive the estimates. Besides, the Metropolis-Hasting (MH) algorithm is applied to get the Bayesian estimates as well. According to Fisher information matrix, we acquire asymptotic confidence intervals. Bootstrap intervals are also established. Furthermore, we build the highest posterior density intervals through the sample generated by the MH algorithm. Then, Bayesian predictive intervals and estimates for future samples are provided. Finally, for evaluating the quality of the approaches, a numerical simulation study is implemented. In addition, we analyze two real datasets.


Introduction
Gompertz distribution has wide applications in describing human mortality, establishing actuarial tables and other fields. Historically, it was originally introduced by Gompertz (see Reference [1]). The probability density function (PDF) and cumulative distribution function (CDF) of the Gompertz distribution are defined as f (x|α, β) = αβe βx e −α(e βx −1) , 0 < x < +∞, (1) and F(x|α, β) = 1 − e −α(e βx −1) , 0 < x < +∞, (2) where the unknown parameters α and β are positive. The Gompertz distribution possesses a unimodal PDF; in addition to this, it also has an increasing hazard function. Many researchers have contributed to the properties of the Gompertz distribution. In recent years, Reference [2] studied the relations between other distributions and Gompertz distribution, for instance, the Type I extreme value and Weibull distributions. Reference [3] obtained the weighted and unweighted least squares estimations under censored and complete samples. Reference [4] calculated the maximum likelihood estimates (MLEs), and completed the establishment for exact confidence interval and joint confidence region base on progressive Type-II censoring. Reference [5] studied the statistical inferences for Gompertz distribution under generalized progressively hybrid censoring. They derived the MLEs by Newton's iteration method and used Markov chain Monte Carlo method to obtain Bayes estimates under generalized entropy and other loss functions. Bayesian predictions based on this censoring scheme were provided by oneand two-sample predictive approaches. Finally, they compared the proposed methods by simulation. Reference [6] obtained the MLEs and Bayesian estimates for the parameters where Q = ( n r )(n − r) , and x = (x r+1 , . . . , x m ) denotes an observed value of X = (X r+1 , . . . , X m ). In addition, ( n r ) is the binomial coefficient, that is n! (n−r)!r! . Substituting (1) and (2), (3) can be written in the following form:

Point Estimation with EM Algrithm
A classical method for obtaining MLE is the Newton-Raphson method, which requires the second-order partial derivatives of the log-likelihood function and the derivatives are usually complicated in the case of censoring. Therefore, it is necessary to seek other methods. Following Reference [16], we use the EM algorithm to derive the MLEs. This algorithm is powerful for handling incomplete data problems because only the pseudo-log likelihood function of complete data needs to be maximized. It is an iterative method by using current estimation of the parameters to expect the log-likelihood function filled with censored data which is called E-step and maximize it to get the next estimation which is called M-step.
We use Z = (Z r , Z r+1 , . . . , Z m ) to represent the censored sample, where Z r is a 1 × r vector Z r = (Z r1 , . . . , Z rr ) of the first r unobserved failures, and Z i , i = r + 1, . . . , m denotes a 1 × R i vector Z i = (Z i1 , . . . , Z iR i ) of the censored data after X i failed. The observed and complete sample are denoted by X = (X r+1 , . . . , X m ) and K respectively, then K = (X, Z). Let z 0 = (z r1 , . . . , z rr ), z i = (z i1 , . . . , z iR i )(i = r + 1, . . . , m) and x = (x r+1 , . . . , x m ) represent the corresponding observations. Under the complete data, we can express the log-likelihood function by L c (α, β, K) = n ln(αβ) + β To conduct the E-step smoothly, first we compute the expectation of (5), the pseudo-log likelihood function is then expressed as and •
The observed information can be extracted by means of a program proposed by Reference [17] when using EM method in handling incomplete sample problems to derive MLEs. Let θ be the unknown parameter (α, β). I K (θ) and I X (θ) denote complete information and observed information respectively. Furthermore, missing information is represented as I K|X (θ). The main concept of this program can be described as the principle of missing information where I K (θ) can be derived by I W|X is the expected information of the distribution of Z = (Z r , Z r+1 , . . . , Z m ) given X = (X r+1 , . . . , X m ). According to Reference [18], given the observed sample of general progressive Type-II censoring, we obtain the distribution of Z as then the I K|X is Let ]. Following Reference [18], we know that given the observed sample (X r+1 , . . . , X m ) = (x r+1 , . . . , x m ), the components of Z r are independent of each other and have the PDF f Z 0 (z rk |x r+1 , θ), k = 1, . . . , r. Similarly, the components of Z i , i = r + 1, . . . , m are independent of each other and have the PDF f Z i (z ik |x i , θ), k = 1, . . . , R i . Therefore, the I K|X can be restated as where and Now we can figure out the elements of the above matrices as follows: and Further, using asymptotic normality of MLEθ = (α,β),θ ∼ N(θ, I −1 (θ)), the 100(1 − γ)% ACIs for the two unknown parameters are obtained as α − η γ 2 Var(α),α + η γ 2 Var(α) and β − η γ 2 Var(β),β + η γ 2 Var(β) , (24) where η γ 2 represents the upper γ 2 -th quantile for standard normal distribution, Var(α) and Var(β) denote the principal diagonal elements of I −1 (θ) respectively.

Bootstrap Confidence Interval
As is widely known, asymptotic confidence interval on the basis of MLE requires a large sample to support its accuracy but, in many practical cases, the sample size tends to not be enough. Reference [19] proposed the bootstrap method to construct the confidence interval (CI), which is more suitable for small sample. In this part, the parametric bootstrap method is employed in the establishment of percentile bootstrap (bootstrap-p) and bootstrap-t CIs for a parameter λ (here λ is α or β). Interested readers may refer to References [20,21] for more information about bootstrap and Reference [22] for the algorithm of generating the general progressive type-II censored sample.

Bayesian Estimation
Bayesian statistics are different from traditional statistics in that they allow the incorporation of subjective prior knowledge about life parameters into the inferential procedure in reliability analysis. Therefore, for the same quality of inferences, Bayesian methods tend to require fewer sample data than traditional statistical methods do. This makes it extremely important in expensive life tests.
We investigate the Bayesian estimates in this section. Suppose that α and β independently have gamma prior distributions with the parameters (a, b) and (c, d). Afterwards, we can obtain their joint prior distribution, that is where the positive constants a, b, c and d are hyperparameters. Let x = (x r+1 , . . . , x m ) be an observed value of X = (X r+1 , . . . , X m ). Based on the joint prior distribution and likelihood function, the joint posterior function is It is clear that (26) is analytically tricky. Furthermore, the Bayesian estimation of a function with α and β is also intractable because it is related to a ratio of two integrals. For solving the corresponding ratio of two integrals, some approximate approaches have been presented in the literature. Among them, the TK method was proposed by Reference [23] to obtain the approximate posterior expectations. Besides, the MH algorithm is a simulation method with wide applications in sampling from posterior density function. In this article, we use the TK method and the MH algorithm to derive approximate explicit forms for the Bayesian estimates.

Loss Functions
In Bayesian statistics, the selection of loss function is a fundamental step. There are many symmetric loss functions, among which squared error loss (SEL) function is wellknown for its good mathematical properties. Let δ be a Bayesian estimate of θ. The form of SEL function is then under SEL function the Bayesian estimate for θ can be computed by δ S = E (θ|X) (θ|X). However, in many practical situations, overestimation and underestimation result in different losses, and the consequence is likely to be quite serious if one uses symmetric loss function indiscriminately. In these cases, asymmetrical loss functions are considered to be more suitable. In the literature, many different asymmetric loss functions were used. Among them, LINEX is dominant, and this loss function can be expressed as Without loss of generality, here, we take ζ = 1. Thus Bayesian estimate for θ is given by the expression [24] proposed a balanced loss function that has a more generalized form where δ 0 is a known estimate of θ such as the MLE, and ρ is a loss function selected arbitrarily. By choosing ρ as SEL function given by (27), the (29) is transformed into balanced squared error loss (BSEL) function. According to BSEL we can give the Bayesian estimate by δ BS = (1 − σ)E (θ|X) (θ|X) + σδ 0 . It can be clearly seen that the balanced loss function is more general since it includes special cases of the MLE, symmetric and asymmetric loss functions. For instance, by setting δ 0 to be the MLE of parameter, under BSEL function, the Bayesian estimate is exactly equal to MLE if σ = 1, and it is simplified into the Bayesian estimate under SEL function when σ = 0. Similarly, if ρ is chosen as LINEX loss function given by (28), L 3 (·) is called BLINEX function. When σ = 1 and σ = 0, the Bayesian estimates under the BLINEX loss function correspondingly reduce to MLE and the case of LINEX loss function.
In this article, we derive Bayesian estimates under SEL, BSEL functions and LINEX loss functions, respectively. Next, the TK method is suggested to deal with the ratio of the integrals problem on posterior expectation estimation.
We set: Maximizing ϕ(α, β) and ϕ * u (α, β) individually, we derive (α 1 ,β 1 ) and (α u ,β u ). Then, the approximate posterior expectation of u(α, β) by applying the TK method is obtained aŝ where | ∑ | and | ∑ * u | represent the corresponding determinants of negative inverse Hessian matrix of ϕ(α, β) and ϕ * u (α, β). Next, ignoring the constant term, we note that Now, we compute the partial derivatives of ϕ: and ∂ϕ ∂β Similarly, the second derivatives can be derived as and Through (36)-(38), | ∑ | is obtained as and As a result, | ∑ * | is Finally in the above calculation processes, setting u(α, β) as α and β respectively, the estimates on the basis of SEL function are given bŷ Further, we are able to calculate the Bayesian estimates based on BSEL function using the equation Similarly, by treating u 1 (α, β) as e −hα and u 2 (α, β) as e −hβ , the estimates for unknown parameters under LINEX loss function are given bŷ

MH Algorithm
We derive Bayesian estimates for α and β by the MH algorithm (see Reference [20]). First, we suppose that the bivariate normal distribution is the proposal distribution for the parameter θ = (α, β), then the MH algorithm can generate samples from the bivariate normal distribution and finally get convergent samples from the posterior distribution. Using the samples, we first compute the Bayesian estimates under different loss functions, thereafter, establish HPD intervals. The MH algorithm can be summarized as: (1) Begin with an initial value θ 0 = (α 0 , β 0 ), set n = 1.
Proceeding similarly, the desired estimates under BSEL function can be obtained easily. Further the Bayesian estimates under LINEX can be computed as Now, we can establish the 100(1 − γ)% HPD interval (see Reference [25]) for the unknown parameter α. Sort the remaining D − D 0 values in ascending order to be α (1) , α (2) , . . . , α (D−D 0 ) . The 100(1 − γ)% HPD interval of α is given as where w * is selected when the following equation is satisfied: . Likewise, the HPD interval of β can be obtained .

Bayesian Prediction
Now we obtain the prediction estimates for the future sample on the basis of available sample and obtain predictive intervals. Bayesian prediction for future sample is a fundamental subject in many fields such as medical, agricultural and engineering experiments. Interested readers may refer to Reference [11].
Suppose that the existing X = (X r+1 , . . . , X m ) is a group of general progressive censored data observed from a population with Gompertz distribution. Let Y 1 ≤ Y 2 ≤ · · · ≤ Y W denote the ordered failures time for a future sample with size W, which is also obtained from Gompertz distribution. We aim to obtain their predictive estimation (two-sample prediction). Suppose that Y v (1 ≤ v ≤ W) represents the v-th failure time of the future sample. Then, for given α and β, we can obtain the density function of Y v as Consequently, the posterior predictive density function for Y v is derived as It is infeasible to compute (57) analytically. By using the MH algorithm, we can obtain its approximate solution as follows: Further, the survival function is computed as The posterior predictive survival function for Y v can be derived by Then, we construct the 100(1 − γ)% Bayesian predictive interval (L 0 , U 0 ) of Y v by finding the solution of the equations Further, it is convenient to derive the predictive estimate of the future v-th ordered lifetime, which is given bŷ By using the MH algorithm described in the previous section, the prediction estimate of Y k is derived asŷ

Simulation and Data Analysis
For evaluating the quality of the approaches, a numeric simulation study is carried out. In addition, we also analyze two real data sets for further illustration.

Simulation Study
For the sake of simulation, first we generate general progressive censored sample with the algorithm discussed by Reference [22]. The procedures are as below: (1) Generate B m from Beta(n − r, r + 1).
(2) Generate independent U r+k from Uniform(0, 1), k = 1, 2, . . . , m − r − 1. Then we get the desired general progressive censored data X i , i = r + 1, . . . , m drawn from Gompertz distribution. In our experiment, the true values for (α, β) are selected to be (0.3, 1.2). The MLEs for the two parameters are calculated by means of EM algorithm. In the aspect of Bayesian estimation and prediction, (0.2, 7.8, 0.1, 3.7) are chosen to be the values of hyperparameters (a, b, c, d) respectively. Moreover, Bayesian estimates are obtained by TK and MH methods under different loss functions. Comparison between the results is made according to mean-square error (MSE).
Furthermore, different intervals have also been constructed, including ACIs on the basis of Fisher information matrix, parametric bootstrap intervals, and HPD intervals based on the sample generated from the MH algorithm. Table 2 presents their average length (AL) and coverage probabilities (CPs). The tabulated values indicate that the AL of the HPD intervals is the shortest among those obtained from other interval estimates. Besides, we also find that the ACIs have better performance according to CPs. In general, bootstrap-t and bootstrap-p intervals behave similarly, and their CPs tend to be below the 95% confidence level. Table 3 lists the results of point prediction and 95% interval prediction. We give the prediction results of y 3 , y 7 and y 10 in a future sample with size 10. Furthermore, we discover that the interval length becomes wider as v increases.

Data Analysis
Dataset 1: First we analyze a real dataset about the breaking stress of carbon fibers (in Gba) (n = 66) (see Reference [26]). It is listed as follows: 3 In order to analyze these data, we calculate the MLEs for the two parameters, and then for the Gompertz distribution we conduct a goodness-of-fit test with some practical guidelines like the Kolmogorov-Smirnov (K-S) statistics, and other criteria, for example, the Akaike Information Criterion (AIC), as well as the Bayesian Information Criterion (BIC). For comparison, some other life distributions have also been tested for goodness-of-fit, including Generalized Exponential (GE), Inverse Weibull and Exponential distributions. Their PDFs have the following forms, respectively: (1) The PDF of GE distribution: (2) The PDF of Inverse Weibull distribution: (3) The PDF of Exponential distribution:  The results of the tests are presented in Table 4 together with the MLEs. Note that the distribution is more suitable to fit the data when K-S, AIC and BIC values are smaller.
Comparing the values, we can conclude that the Gompertz distribution is more appropriate.
To illustrate the proposed methods, three groups of general progressive censored data have been randomly drawn from the parent sample as follows:  With the EM algorithm, we calculate the MLEs, and the corresponding Bayesian estimates are also derived by TK and MH methods. For the sake of no prior information, all the hyperparameters are set close to zero values. We list the average MLEs and Bayesian estimates in Tables 5 and 6. In Table 7, the 90% interval estimates are tabulated, which are ACIs, parametric bootstrap and HPD intervals. Finally in Table 8, point prediction and 95% interval prediction of y 1 and y 6 in a future sample with size 6 are presented. In order to analyze these data, Reference [28] assumed that the number of tumor-free days obeys the Gompertz distribution. To illustrate the methods discussed, here we also suppose that the distribution for these data is Gompertz with (α, β). Let m − r = 20, we set up two censoring schemes, respectively s 1 = (0 * 4 , 5, 0 * 7 , 2, 0 * 7 ), r = 3 and s 2 = (0 * 2 , 4, 0 * 6 , 5, 0 * 10 ), r = 1. Then we have obtained the sample under s 1 :   63, 66, 66, 66, 68, 70, 70, 77, 77, 84, 91, 91, 94, 98, 101, 105, 108, 112, 115 In Tables 9 and 10, we calculate the average MLEs, and average Bayesian estimates are derived by the TK method and the MH algorithm, respectively. The interval estimates are presented in Table 11, including ACIs, parametric bootstrap and HPD intervals. Finally Table 12 presents the results of the point prediction and the 95% interval prediction of y 1 and y 5 with W = 5.

Conclusions
In summary, we discuss the classical and Bayes inferences for the Gompertz distribution using the general progressive censoring. First, the MLEs are acquired by the Expectation-Maximization algorithm. Then, according to the asymptotic normality of MLEs and the principle of missing information, we provide the asymptotic confidence intervals. Moreover, we derive parametric percentile bootstrap and bootstrap-t intervals. In Bayesian statistics, three loss functions are considered, which are symmetrical, asymmetrical and balanced, respectively. Since the posterior expectation is intractable to obtain in explicit form, the TK method is employed to calculate approximate Bayesian estimates. Besides, the Metropolis-Hasting algorithm is applied to get the Bayesian estimates and establish HPD intervals. Furthermore, we derive the prediction estimates of future samples. Finally, a numerical simulation is executed to appraise the quality of the approaches, and two real data sets are also analyzed. The results indicate that these approaches have good performance. In addition, the methods in this article can be extended to other distributions.

Conflicts of Interest:
The authors declare no conflict of interest.