Bayesian Inference for the Kumaraswamy Distribution under Generalized Progressive Hybrid Censoring

Incomplete data are unavoidable for survival analysis as well as life testing, so more and more researchers are beginning to study censoring data. This paper discusses and considers the estimation of unknown parameters featured by the Kumaraswamy distribution on the condition of generalized progressive hybrid censoring scheme. Estimation of reliability is also considered in this paper. To begin with, the maximum likelihood estimators are derived. In addition, Bayesian estimators under not only symmetric but also asymmetric loss functions, like general entropy, squared error as well as linex loss function, are also offered. Since the Bayesian estimates fail to be of explicit computation, Lindley approximation, as well as the Tierney and Kadane method, is employed to obtain the Bayesian estimates. A simulation research is conducted for the comparison of the effectiveness of the proposed estimators. A real-life example is employed for illustration.


Kumaraswamy Distribution
Given that the classical probability distribution functions like beta, normal, log-normal, Student-t (Contreras-Reyes et al. [1]) and other empirical distributions cannot fit hydrological data quite well, Kumaraswamy [2] came up with a new two-parameter distribution to be specifically applicable for hydrological problems. The cumulative distribution function (cdf) of the Kumaraswamy distribution is where both α and β represent the positive shape parameters of the distribution, which is denoted by K(α,β) in this paper. The following presents the corresponding probability density function (pdf): In addition, an in-depth observation is made on the reliability function of the Kumaraswamy distribution, which is shown below: An increasing number of statisticians have started to study it since it has many flexible shape properties which are identical with the beta distribution. In light of different values of parameters, What is more, the Kumaraswamy distribution fits well with some natural phenomena, such as daily rainfall, water flows and other pertinent fields, see Fletcher and Ponnambalam [4], Sundar and Subbiah [5], Ponnambalam et al. [6] as well as Seifi et al. [7], especially for the outcomes of which possess upper and lower bounds, like the heights of people, test scores, air temperatures, economic data, etc. Estimation for the Kumaraswamy distribution has gradually attracted the attention of scholars in recent years. Lemonte [8] obtained modified maximum likelihood estimators which are unbiased in the second stage. Based on the bias correction estimations, they studied a bias-corrected method using parametric bootstrap. Then, statisticians attempted to study the Kumaraswamy distribution under different censoring schemes and started to combine classical and Bayesian methods. Ghosh and Nadarajah [9] discussed Bayesian estimation using two loss functions under three types of censoring schemes: left censoring, single type-II censoring and double type-II censoring with one parameter known. Sultana et al. [10] discussed and estimated parameters of the Kumaraswamy distribution with hybrid censoring scheme, and recently Sultana et al. [11] combined hybrid with progressive type I censoring schemes to explore the parameter estimation problems for the same distribution.

Generalized Progressive Hybrid Censoring Scheme
Life testing experiments are widely used in engineering, biology, machinery and other fields of study, which can be summarized as mathematical and probabilistic models of survival analysis. In reality, several restrictions, such as time and cost, prevent us from observing the failure time of all units. It is common to cease in the middle of process before all the observations fail. Such limitations result in censored data. Among all censoring cases, the two most typical schemes are Type-I as well as Type-II censoring schemes. According to previous literature, plenty of authors have discussed this aspect and one may consult Meeker and Escobar [12] which includes methods of handling Type-I as well as Type-II censored data.
An accidental pause or an unavoidable loss of the experiment units is likely to happen before the final termination. However, the constraint in those two censoring schemes is that removal cannot occur to the units in the duration of the experiment. To resolve this inflexibility, Cohen [13] first introduced progressive censoring scheme. A progressive Type-II censoring sample will be given below. Assuming that n independent units of a common lifetime distribution denoted by X 1 , X 2 , . . . , X n are put in the experiment at t = 0, when the first failure occurs, among n − 1 survivals, we take R 1 units out of the experiment at random. Similarly, when the second failure happens, we randomly remove R 2 units among the n − 2 − R 1 survivals. We conduct the repeated procedure till the observation of the mth failure. On the occasion of the m-th failure, removal occurs to the R m = n − m − R 1 − R 2 − . . . − R m−1 survivals. During the experiment, the progressive censoring scheme (R 1 , R 2 , . . . , R m ), which is considered to be progressive type-II censored scheme, is prefixed satisfying ∑ m i=1 R i + m = n. m ordered failure times are written in light of X 1 < X 2 < . . . < X m .
According to weakness harbored by the progressive type-II censored scheme, if the experimental units are highly reliable, this experiment will last quite long. Hence, the progressively hybrid censoring scheme was introduced by Kundu and Joarder [14]. For the censoring scheme, the implementation of n independent identical distributed units is used. The experimenter will cease the operation at min {T, X m }. Here, the time T as well as 1 ≤ m ≤ n is determined ahead of time. In the context of the progressive type-II censored scheme, the span of the experiment will not take a longer duration than T.
However, given that when the prefixed termination time T may be small, the observation we obtained would be insufficient. Therefore, a new mode of censoring scheme-generalized progressive hybrid censoring scheme-is proposed by Cho and Sun [15], which enables us to obtain a predetermined series of failures. How to get generalized progressive hybrid censored data is described below by a graphic illustration in Figure 3.

I 0
Removed Numbers: End t (a) when T < X k < X m , the experiment stops at k-th failure time and obtains X 1 , . . . , X k .

II 0
Removed Numbers: the experiment stops at time T and obtains X 1 , . . . , X J .

III 0
Removed Numbers: (c) when X k < X m < T, the experiment stops at m-th failure time and obtains X 1 , . . . , X m . Assume that our research group possesses n independent units of a common lifetime distribution. The corresponding lifetime is denoted by X 1 , X 2 , . . . , X n . The integers k and m (k < m) have been under predetermination between zero and n as well as R 1 , R 2 , . . . , R m which can satisfy the equation ∑ m i=1 R i + m = n function as preplanned integers. On the arrival of the first failure X 1 , we randomly remove R 1 units. When the second failure X 2 happens, we take random removal of R 2 units out of the n − 2 − R 1 survivals. The process is repeated and terminated at T * = max{min{T, X m }, X k } with the rest of the survival units under the removal. It greatly modified the previous schemes so that we can choose to continue the experiment when the sample is insufficient at the prefixed cut-off time T. On the condition of the generalized progressive hybrid censoring scheme, researchers would like to obtain m failures, while they can also adopt k failures which are regarded as the bare minimum. We denote the generalized progressive hybrid censoring scheme as (R 1 , R 2 , . . . , R m ). Let J be the observed failure times before arriving at the predetermined time T. The generalized progressive censoring scheme can be classified into these cases as below: Case I : X 1 , . . . , X J , . . . , X k , for T < X k < X m , Case II : X 1 , . . . , X k , . . . , X J , for X k < T < X m , Case III : X 1 , . . . , X k , . . . , X m , for X k < X m < T.
Case III represents the progressive Type-II censoring scheme and the mixture of Case II and Case III is the progressive hybrid censoring scheme. Hence, evidently, Case I is the modification of this scheme. Let R 1 = 0, R 2 = 0, . . . , R m = 0, we can obtain complete sample. We also assume R 1 = . . . = R m−1 = 0, R m = n − m, in order to get the type-II censored sample.
According to the generalized progressive hybrid censoring scheme, the likelihood equations of the three cases are where As far as we know, estimations for the Kumaraswamy distribution under generalized progressive hybrid censoring scheme have not been done yet in previous literature. Hence, this paper will discuss the estimations for parameters and reliability of the Kumaraswamy distribution on the basis of this model.
The rest of the article consists of the following parts. In the next section, the maximum likelihood estimators of the two unknown parameters as well as reliability function of the model will be derived theoretically. For all unknown quantities, Bayes estimators are achieved in Section 3 under three diverse loss functions employing Lindley's approximation. Then in Section 4, a simulation experiment will be carried out in light of the conclusion in Sections 2 and 3. Data analysis is demonstrated in Section 5. Finally, Section 6 concludes the thesis.

Maximum Likelihood Estimation
In dealing with reliability problems and survival analysis, an effective and classical approach widely employed by statisticians is maximum likelihood estimation (MLE). By employing this method, two unknown parameters will be derived. As a result, MLE of R(t) will also be acquired. Plugging the pdf and cdf of the Kumaraswamy distribution, i.e., (2) and (1), into the likelihood Formula (4), the likelihood functions of α and β after neglecting the constants are expressed as Case I : Case II : Case III : Disregarding the constant, the log-likelihood functions are Case I: Case II: Case III: In order to simplify the above expressions, we combine Case I, II and III and obtain the log-likelihood function as where for Case I, D = k, E(α) = 0; for Case II, D = J, E(α) = R * J+1 log(1 − T α ); for Case III, D = m, E(α) = 0.
We take the partial derivatives of the above function (5) for α and β respectively and get a set of likelihood equations as follows: and ∂l ∂α where E (1) (α) = 0 for Case I and III, By solving the roots of the above set of equations, the MLEs of two parameters are able to be attained theoretically. From the Equation (6), we gain the maximum likelihood estimate of β aŝ Placing the estimation value of β into the Equation (7), we can get where Obviously, it is hard to simplify and gain solutions of closed forms, since both (8) and (9) are nonlinear. In this situation, updating the estimates seems to be an effective method to gain approximate solution of α. This iterative algorithm has been proposed by Kundu [16]. Here, we give a brief description. Begin with an initial assumption of α, noted by α (0) , then gain α (1) = g(α (0) ) and repeat this iteration and we gain α(n + 1) = g(α (n) ). When the precision meet the tolerance limit which is set beforehand |α (n+1) − α (n) | < ε, we stop the iterative process. Once we get the MLE of α, denoted byα, the MLE of β is deduced asβ =β(α). As long as the MLEs are obtained, we can substitute these two estimates into (3) to gain the MLE of R(t) as:

Bayesian Estimation
Bayesian estimation which considers prior information as well as sample information is a fresh but efficient approach in comparison with MLE. This more comprehensive estimation method is usually more precise than the maximum likelihood estimation. In this part, Bayesian estimation of the model parameters α, β and reliability function will be discussed.

Symmetric and Asymetric Loss Functions
In statistics, we usually estimate a parameter by minimizing the loss function. There are lots of diverse symmetric or asymmetric loss functions. Here we will interpret the three typical loss functions which are taken into consideration. In all of the following cases, d(η) stands for the true value of the unknown parameter andd(η) is the corresponding estimate of d(η). The symmetric one refers to the squared error loss (SEL) function. It is the most prevalent one and can be easily proved to be right based on minimum variance-unbiased estimation. The definition and corresponding Bayes estimator are However, due to symmetry, the overestimation of SEL function has equal weight as underestimation of the same magnitude, which gives rise to the emergence of a large number of asymmetric functions. The Linex loss (LL) function [17], an extensively adopted asymmetric loss function, is another loss function discussed in this paper. The definition and corresponding Bayes estimator are The parameter p represents the deviation direction, and the degree of deviation is reflected by its magnitude. When p < 0, the underestimation is greater than the overestimation and the opposite is the case when p > 0. When parameter p converges towards zero, the linex loss function can be converted to SEL loss function.
In addition, an asymmetric loss function-the general entropy loss (EL) function is also considered whose definition and corresponding Bayes estimator are Here, the positive error is greater than the negative one when q > 0, and the opposite is the case when q > 0.

Prior and Posterior Distributions
Since a natural conjugate bivariate prior distribution for α and β does not exist, we employ the same assumption as Kundu and Pradhan [18], supposing that α and β are subjected to gamma distributions independently for the reason that gamma distribution can be adapted to various shapes depending on parameter values. Now, the joint prior distribution, ignoring the constant coefficient, is in the form of Here, the positive hyperparameters a, b, c, d embody the prior knowledge and information about α and β. On this basis, the joint posterior distribution is where X are observations X 1 , X 2 , . . .. The conditional posterior distribution is defined by Under the SEL function, the Bayes estimators can be gained aŝ where ψ stands for α, β or reliability function. Under the LL function, the Bayes estimators are obtained aŝ Under the EL function, the Bayes estimators are written aŝ By observing the above Bayes estimation expressions, we find that they are all the ratios of two integrals and the explicit expressions are hard to acquire. Therefore, proper methods to approximate the above integrals need to be employed. Thus, we introduce Lindley's approximation, as well as the Tierney and Kadane method, to get the closed estimators.

Lindley's Approximation
Lindley [19] deduced this general term formula by developing asymptotic expansions for the ratio of integrals. To apply this method, we denote µ = (µ 1 , µ 2 ) and consider the function g(µ 1 , µ 2 ) in the form of ratio of intergrals given by where g(µ 1 , µ 2 ) is a certain function of µ 1 and µ 2 , l(µ 1 , µ 2 | X) is the log-likelihood function and κ(µ 1 , µ 2 ) = log π(µ 1 , µ 2 ). Employing the Lindley method, the function g(µ 1 , µ 2 ) can be written aŝ where For our estimation problem, µ = (α, β) and now let us deduce it in detail. For convenience, we denote Also, for our problem we have for Case I and III, for Case II, for Case I and III, for Case II.
• Under the SEL function, -when g(α, β) = α, we observe that Using the above Equation (14), the Bayes estimator of α can be obtained aŝ -When g(α, β) = β, we can derive that Likewise, the Bayes estimator of β can be gained aŝ The Bayes estimator of R(t) can be computed bŷ • Under the LL function, -when g(α, β) = e −pα , we observe that The Bayes estimator of α can be obtained aŝ -When g(α, β) = e −pβ , we can derive that Similarly, the Bayes estimator of β is in the form of The Bayes estimator of R(t) are acquired aŝ • Under the EL function, -when g(α, β) = α −q , we observe that The Bayes estimator of α can be obtained aŝ -When g(α, β) = β −q , we can derive that Analogously, the Bayes estimator of β is in the form of The Bayes estimate of R(t) can be achieved bŷ
For the two-parameter Kumaraswamy distribution, we have Subsequently, the following set of equations need to be solved to get (α δ ,β δ ): Then, we obtain |Ω|, which is given by Recall that expressions |Ω * g | and δ * g in Equation (16) rely on function g.
• Under the SEL function, -g(α, β) = α for α and corresponding function δ * comes to be Later, we solve the set of following equations and gain (α δ * ,β δ * ). Then, |Ω * α | is computed as Using the above expression in Equation (16), the desired Bayes estimator of α can be written aŝ .
The Bayes estimator of β based on the SEL function can be attained in similiar way. -Now, we consider the relibility function R(t). Let g(α, β) = (1 − t α ) β , then we have Hence, we calculate (α δ * ,β δ * ) by figuring out the following set of equations: Subsequently, we deduce |Ω * Rt | as follows: After that, the Bayes estimator of reliability function turns out to bê . •
Likewise, based on the linex loss function, the Bayes estimator of β will be realized. -Now, we consider the relibility function R(t). Let g(α, β) = e −p(1−t α ) β , then we have Hence, we calculated (α δ * ,β δ * ) by figuring out the following set of equations Subsequently, we deduce |Ω * Rt | as follows: After that, the Bayes estimator of R(t) iŝ . •
Obviously, under the EL function, we can also obtain the Bayes estimator of β. -Now, we consider the relibility function R(t). Let g(α, β) = (1 − t α ) −qβ , then we have Hence, we calculate (α δ * ,β δ * ) by figuring out the following set of equations: Subsequently, we deduce |Ω * Rt | as follows: After that, the Bayes estimator of R(t) results to bê .

Simulation Study
Within the simulation experiment, generating the generalized progressive hybrid censoring sample is our first step. Before continuing further, we give the way to generate progressive Type-II censoring sample in accordance with Balakrishnan and Sandhu [22]. He presented a simple but effective simulational algorithm, which enables one to collect a series of progressive Type-II censored sample out of any continuous distribution. By adapting this method, we give the algorithm of generating the generalized progressive hybrid censoring sample as below.
Step-4: At last we set X i = F −1 (Y i ), for i = 1, . . . , m, where F −1 (Y i ) stands for the inverse culmulative density function of any distribution considered. X 1 , X 2 , . . . , X m are the desired progressive Type-II censored sample out of the distribution F(). In this article, F() is the Kumaraswamy distribution.
Without loss of generality, it is found to use the Kumaraswamy distribution which takes the value of α = 3, β = 2 and T = 0.9. The process is replicated 1000 times in each case. The results have been obtained under diverse (n, m, k). The censoring schemes employed in this simulation are We compute both the MLE and Bayes estimates. Bayes estimates are respectively on the condition of the non-informative prior distribution, which means that four hyperparameters a, b, c, d adopt values of 0, and informative prior distribution, where a = 1.5, c = 1, and b = d = 0.5. Bayes estimates are obtained under loss function subject to SEL, LL and EL functions. As for linex loss function, relative estimates are gained with p = 0.2, 0.8. General entropy loss function is considered with q = 0.2, 0.8. Finally, we use mean squared errors (MSEs), as well as absolute biases (ABs), to evaluate the accuracy of the estimations. They are where σ is the true value,σ stands for the corresponding estimate, and S represents the simulation times. The simulation results, Tables A1-A3, are shown in the Appendix A. In every two rows of Tables A1-A3, the upper one is MSEs and the lower one is ABs. We put the corresponding MSEs and ABs of MLE in the fifth column. Besides, all other columns are composed of eight values. The first two values represent the MSEs and ABs in light of Lindley approximation with noninformative prior distribution. The second two values represent the MSEs and ABs using Lindley approximation with informative prior distribution. The third two values denote the MSEs and ABs using TK method with noninformative prior distribution. The fourth two values correspond to the MSEs and ABs using TK method with informative prior distribution.
In Table A1, it can be seen that with the sample size n increasing, MSEs decrease, since as n increases, more additional information is gathered. For a given sample size, MSEs also decline with the generalized progressive hybrid censored sample R i decreasing, namely m increasing. Further, with n and m kept constant, the MSEs decline as the accepted bare minimum of failures k increases. In general, as for MSEs, Bayes estimates are superior to the MLEs, and there is no significant difference among the three schemes. Particularly, respective Bayes estimates of α under SEL make a better performance compared to the corresponding LL and EL. Besides, the Bayes estimates with informative prior provide better performance compared with the respective MSEs with noninformative prior. Estimation under LL with p = 0.2 gives a better option than p = 0.8, while q = 0.2 tends to produce lower MSEs than q = 0.8 under EL. Overall, with proper prior information, the Bayes estimates based on the SEL function and TK method do better than other estimates for α.
In Table A2, it can be seen that with the sample size n increasing, MSEs decrease, since as n increases, more additional information is gathered. For a given sample size, MSEs also decline with the generalized progressive hybrid censored sample R i decreasing, namely m increasing. Further, with n and m kept constant, the MSEs decline as the accepted bare minimum of failures k increases. In general, as for MSEs, Bayes estimates are superior to the MLEs, and there is no significant difference among the three schemes. Particularly, respective Bayes estimates of β under LL perform better than the corresponding SEL, as well as EL. Besides, the Bayes estimates with informative prior provide better performance compared with the respective MSEs with noninformative prior. Estimation under LL with p = 0.2 gives a better option than p = 0.8, while q = 0.2 tends to produce lower MSEs than q = 0.8 under EL. As for β, the Bayes estimates with proper prior based on the lines loss function and Lindley method perform better than other estimates.
In Table A3, it can be seen that with the sample size n increasing, MSEs decrease, since as n increases, more additional information is gathered. For a given sample size, MSEs also decline with the generalized progressive hybrid censored sample R i decreasing, namely m increasing. Further, with n and m kept constant, the MSEs decline as the accepted bare minimum of failures k increases. The maximum likelihood estimates outperform noninformative Bayes estimates, and there is no significant difference among the three schemes. Particularly, respective Bayes estimates of R(t) under LL perform worse than the corresponding SEL, as well as EL. Besides, the Bayes estimates with informative prior provide better performance compared with the respective MSEs with noninformative prior. In the case of LL, the choice p = 0.8 is preferred, while in the case of EL, the option q = 0.2 produces better results.

Data Analysis
This part illustrates a real-life dataset, aiming at analyzing the validity of the presented estimation methods. The following dataset shows the monthly water capacity of the Shasta reservoir in the time range of August and December from 1975 to 2016. These data have been used before by some statisticians, such as Kohansal [23]. Because K(α, β) is defined for 0 < x < 1, all data are divided by the total capacity of Shasta reservoir 4,552,000 acre-foot. The transformed data are tabulated in Table 1 and presented as a histogram in Figure 4.  One may doubt whether the considered dataset comes from the Kumaraswamy distribution or not. To verify the reasonableness, we fit the Kumaraswamy distribution to the dataset, competing with two other distributions-exponentiated exponential distribution (EED) and Lomax distribution (LD). The corresponding cdf and pdf are given below.
We employ the Akaike's information criterion (AIC), defined by −2ln(L) + 2k, the associated second order criterion (AICc), defined by AIC + 2k(k+1) n−k−1 and Bayesian information criterion (BIC), defined by −2ln(L) + kln(n), where L is maximized value of the likelihood function, k is the number of the parameters, n is the number of observations, as well as Kolmogorov-Smirnov (K-S) statistics with its p-value. The results are shown in Table 2. Besides, the empirical cumulative distribution functions and the quantile-quantile plots are given respectively in Figures 5 and 6.  Considering that the Kumaraswamy distribution has lower AIC, AICc, BIC and K-S statistics and higher p-value, it is reasonable to say we fail to reject the hypothesis that the data come from Kumaraswamy.
The following data is generalized progressive hybrid censored sample artificially by using m = 21 and R 1 = R 2 = . . . = R 21 = 1. The ordered progressive Type-II censored dataset is given in Table 3.
In this illustration, take T = 0.75 and k = 19 for Case I, T = 0.75 and k = 14 for Case II, as well as T = 0.9 and k = 14 for Case III. Table 4 presents the estimation of α, β, as well as the reliability function of the generalized progressive hybrid censored sample.

Conclusions
In this paper, we estimate two parameters and the reliability of the Kumaraswamy distribution when the dataset is sampled in a generalized progressive hybrid censoring scheme. The maximum likelihood estimators are inferred. In addition, the Bayes estimators are achieved by the Lindley method, as well as the TK method, based on general entropy, squared error and linex loss function, which are all conducted in light of noninformative and informative priors. The performance of the estimates is contrasted according to absolute bias values and mean squared error. The Bayes estimates of proper informative prior are revealed to work better than corresponding noninformative prior estimates. Besides, Bayes estimates produce lower MSEs for two parameters, while for reliability estimation MLEs have lower values instead. A real-life example is also intensively investigated. In the future, more in-depth researches are worth discussing, such as the use of bayesian estimation for the Kumaraswamy distribution in non-linear regression models (Contreras-Reyes et al.