Bayesian Estimation and Prediction for Flexible Weibull Model under Type-II Censoring Scheme

Abstract

We have developed the Bayesian estimation procedure for flexible Weibull distribution under Type-II censoring scheme assuming Jeffrey's scale invariant (noninformative) and Gamma (informative) priors for the model parameters. The interval estimation for the model parameters has been performed through normal approximation, bootstrap, and highest posterior density (HPD) procedures. Further, we have also derived the predictive posteriors and the corresponding predictive survival functions for the future observations based on Type-II censored data from the flexible Weibull distribution. Since the predictive posteriors are not in the closed form, we proposed to use the Monte Carlo Markov chain (MCMC) methods to approximate the posteriors of interest. The performance of the Bayes estimators has also been compared with the classical estimators of the model parameters through the Monte Carlo simulation study. A real data set representing the time between failures of secondary reactor pumps has been analysed for illustration purpose.

1. Introduction

In reliability/survival analysis, generally, life test experiments are performed to check the life expectancy of the manufactured product or items/units before products produced in the market. But in practice, the experimenters are not able to observe the failure times of all the units placed on a life test due to time and cost constraints or due to some other uncertain reasons. Data obtained from such experiments are called censored sample. Keeping time and cost constraints in mind, many types of censoring schemes have been discussed in the statistical literature named as Type-I censoring, Type-II censoring and progressive censoring schemes, and so forth. In this paper, Type-II censoring scheme is considered. In Type-II censoring scheme, the life test is terminated as soon as a prespecified number (say, 𝑟) of units have failed. Therefore, out of 𝑛 units put on test, only first 𝑟 failures will be observed. The data obtained from such a restrained life test will be referred to as a Type-II censored sample.

Prediction of the lifetimes of future items based on censored data is very interesting and valuable topic for researchers, engineers and reliability practitioners. In predictive inference, In predictive Inference, we can infer about the lifetimes of the future items using observed data. The future prediction problem can be classified into two types: (1) one-sample prediction problem. (2) two-sample prediction problem. In one-sample prediction problem, the variable to be predicted comes from the same sequence of variables observed and dependent of the informative sample. In the second type, the variable to be predicted comes from another independent future sample. Reference [1] has developed the Bayesian procedure to the prediction problems of future observations and use the concept of Bayesian predictive posterior distribution. Many authors have focussed on the problem of Bayesian prediction of future observations based on various types of censored data from different lifetime models (see [2–9], and references cited therein).

The flexible Weibull distribution is a new two-parameter generalization of the Weibull model which has been proposed by [10]. The density function of the flexible Weibull distribution is given by 𝑓(𝑥)=𝛽𝛼+𝑥²𝛽exp𝛼𝑥−𝑥exp−𝑒^{𝛼𝑥−𝛽/𝑥},𝛼,𝛽>0. (1)

The corresponding CDF is given by 𝐹(𝑥)=1−exp−𝑒^{𝛼𝑥−𝛽/𝑥}. (2)

They have shown that this distribution is able to model various ageing classes of lifetime distributions including IFR, IFRA, and MBT (modified bathtub). They have also checked the goodness-of-fit of this distribution to the failure time between secondary reactor pumps data in the comparison of various extensions of Weibull model and found that this model gives the better fit. Therefore, this distribution can be considered as an alternative lifetime model of the various well-known generalizations of the Weibull distribution. Some statistical properties and the classical estimation procedure for the flexible Weibull distribution have been discussed by [10]. It is to be mentioned here that this distribution has not been considered under Bayesian setup in the earlier literature.

It is well-known that the squared error loss function is the most widely used loss function in Bayesian analysis. This loss is well justified in the classical paradigm on the ground of minimum variance unbiased estimation procedure. In most of the cases, Bayesian estimation procedures has been developed under the same loss function. For Bayesian estimation, we also need to assume a prior distribution for the model parameters involved in the analysis. In this paper, Bayesian analysis have been preformed under the squared error loss function (SELF) assuming both Jeffreys scale invariant and Gamma priors.

A major difficulty to the implementation of Bayesian procedure is that of obtaining the posterior distribution. The process often requires the integration which is very difficult to calculate especially when dealing with complex and high-dimensional models. In such a situation, Monte Carlo Markov chain (MCMC) methods, namely, Gibbs sampling [11] and Metropolis-Hastings (MH) algorithms [12, 13], are very useful to simulate the deviates from the posterior density and produce the good approximate results.

The rest of the paper is organized as fallows. In Section 2, we have discussed the point estimation procedures for the parameters of the considered model under classical set-up. The confidence/bootstrap intervals have been constructed in Section 3. In Section 4, we have developed the Bayesian estimation procedure under the assumption that model parameters have the gamma prior density function. We have also derived the one-sample and two-sample predictive densities and corresponding survival functions of the future observables in Sections 5 and 6, respectively. The predictive bounds of future observations under one-sample and two-sample predictions have also been discussed in respective sections. For comparing the performance of the classical and Bayesian estimation procedures, the Monte Carlo simulation study has been presented in Section 7. To check the applicability of the proposed methodologies, a real data set has been analysed in Section 8. Finally, the conclusions have been given in Section 9.

2. Classical Estimation

Let 𝑥₁,𝑥₂,…,𝑥_𝑛 be the IID random sample from (1) and let 𝑥₍₁₎,𝑥₍₂₎,…,𝑥_(𝑟),(𝑟≤𝑛) be the ordered sample obtained under Type-II censoring scheme. Then, the likelihood function for such type of censored sample can be defined as 𝐿=𝑛!(𝑛−𝑟)!^𝑟_𝑖=1𝑓𝑥_(𝑖)1−𝐹𝑥_(𝑟)^𝑛−𝑟. (3)

Substituting (1) and (2) in (3), we have 𝐿=𝑛!(𝑛−𝑟)!^𝑟_𝑖=1𝛽𝛼+𝑥²_(𝑖)exp^𝑟_𝑖=1𝛼𝑥_(𝑖)−𝛽𝑥_(𝑖)×exp−^𝑟_𝑖=1exp𝛼𝑥_(𝑖)−𝛽𝑥_(𝑖)−(𝑛−𝑟)exp𝛼𝑥_(𝑟)−𝛽𝑥_(𝑟). (4)

The log-likelihood function is given by log𝐿=ln𝑛!()!+𝑛−𝑟^𝑟_𝑖=1𝛽ln𝛼+𝑥²_(𝑖)+^𝑟_𝑖=1𝛼𝑥_(𝑖)−𝛽𝑥_(𝑖)−^𝑟_𝑖=1exp𝛼𝑥_(𝑖)−𝛽𝑥_(𝑖)−(𝑛−𝑟)exp𝛼𝑥_(𝑟)−𝛽𝑥_(𝑟). (5)

The MLEs 𝛼 and 𝛽 of 𝛼 and 𝛽 can be obtained as the simultaneous solution of the following two nonlinear equations: ^𝑟_𝑖=11𝛼+𝛽/𝑥²_(𝑖)+^𝑟_𝑖=1𝑥_(𝑖)−^𝑟_𝑖=1𝑥_(𝑖)𝑒^{(𝛼𝑥(𝑖)−𝛽/𝑥(𝑖))}−(𝑛−𝑟)𝑥_(𝑟)𝑒^{(𝛼𝑥(𝑟)−𝛽/𝑥(𝑟))}=0,_𝑟𝑖=11/𝑥²_(𝑖)𝛼+𝛽/𝑥²_(𝑖)−^𝑟_𝑖=11𝑥_(𝑖)+^𝑟_𝑖=1𝑒^{(𝛼𝑥(𝑖)−𝛽/𝑥(𝑖))}𝑥_(𝑖)+()𝑒𝑛−𝑟^{(𝛼𝑥(𝑟)−𝛽/𝑥(𝑟))}𝑥_(𝑟)=0 (6)(7)

It can be seen that the above equations cannot be solved explicitly and one needs iterative method to solve them. Here, we proposed the use of the fixed point iteration method, which can be routinely applied as follows:

The equation (6) can be rewrite as 𝛼=𝜙₁=∑𝛼,𝛽^𝑟𝑖=11/1+𝛽/𝛼𝑥²_(𝑖)ℎ₁.𝛼,𝛽 (8)

Similarly, from (7), we have 𝛽=𝜙₂=∑𝛼,𝛽^𝑟𝑖=11/𝑥²_(𝑖)/𝛼/𝛽+1/𝑥²_(𝑖)ℎ₂,𝛼,𝛽 (9)

where ℎ₁𝛼,𝛽=^𝑟_𝑖=1𝑥_(𝑖)𝑒^{(𝛼𝑥(𝑖)−𝛽/𝑥(𝑖))}+(𝑛−𝑟)𝑥_(𝑟)𝑒^{(𝛼𝑥(𝑟)−𝛽/𝑥(𝑟))}−^𝑟_𝑖=1𝑥_(𝑖),ℎ₂𝛼,𝛽=^𝑟_𝑖=11𝑥_(𝑖)−^𝑟_𝑖=1𝑒^{(𝛼𝑥(𝑖)−𝛽/𝑥(𝑖))}𝑥_(𝑖)−()𝑒𝑛−𝑟^{(𝛼𝑥(𝑟)−𝛽/𝑥(𝑟))}𝑥_(𝑟). (10)

The following steps followed to obtain the solution of the normal equations (6) and (7).

Step 1. Start with initial starting points, say 𝛼⁰ and 𝛽⁰.

Step 2. By using 𝛼⁰ and 𝛽⁰, obtain 𝛼=𝜙₁(𝛼⁰,𝛽⁰) from (8).

Step 3. Then, obtain 𝛽=𝜙₂(𝛼⁰,𝛽⁰) from (9).

Step 4. If |𝛼⁰−𝛼|≤𝜖 and |𝛽⁰−𝛽|≤𝜖, where 𝜖 is some preassigned tolerance limit, then (𝛼,𝛽) will be the desired solution of (8) and (9).

Step 5. If |𝛼⁰−𝛼|>𝜖 and |𝛽⁰−𝛽|>𝜖, then set 𝛼⁰=𝛼 and 𝛽⁰=𝛽 and repeat Steps 2–5, until tolerance limit is achieved.

3. Confidence Intervals

3.1. Asymptotic Confidence Intervals

The exact distribution of MLEs cannot be obtained explicitly. Therefore, the asymptotic properties of MLEs can be used to construct the confidence intervals for the parameters. Under some regularity conditions, the MLEs (𝛼,𝛽) are approximately bivariate normal with mean (𝛼,𝛽) and variance matrix 𝐼⁻¹(𝛼,𝛽), where 𝐼(𝛼,𝛽) is the observed Fishers information matrix and is defined as 𝐼𝛽=−𝛼,𝜕log𝐿𝜕𝛼²−𝜕log𝐿−𝜕𝛽𝜕𝛼𝜕log𝐿−𝜕𝛼𝜕𝛽𝜕log𝐿𝜕𝛽²_(𝛼=𝛼^,𝛽=𝛽₎, (11)

where 𝜕log𝐿𝜕𝛼²=−^𝑟_𝑖=11𝛼+𝛽/𝑥²_(𝑖)²−^𝑟_𝑖=1𝑥²_(𝑖)𝑒^{(𝛼𝑥(𝑖)−𝛽/𝑥(𝑖))}−(𝑛−𝑟)𝑥²_(𝑟)𝑒^{(𝛼𝑥(𝑟)−𝛽/𝑥(𝑟))},𝜕log𝐿𝜕𝛽²=−^𝑟_𝑖=11/𝑥²_(𝑖)²𝛼+𝛽/𝑥²_(𝑖)²+^𝑟_𝑖=1𝑒^{(𝛼𝑥(𝑖)−𝛽/𝑥(𝑖))}𝑥²_(𝑖)+()𝑒𝑛−𝑟^{(𝛼𝑥(𝑟)−𝛽/𝑥(𝑟))}𝑥²_(𝑟),𝜕log𝐿=𝜕𝛽𝜕𝛼𝜕log𝐿𝜕𝛼𝜕𝛽=−^𝑟_𝑖=11/𝑥²_(𝑖)𝛼+𝛽/𝑥²_(𝑖)²+^𝑟_𝑖=1𝑒^{(𝛼𝑥(𝑖)−𝛽/𝑥(𝑖))}+(𝑛−𝑟)𝑒^{(𝛼𝑥(𝑟)−𝛽/𝑥(𝑟))}. (12)

The diagonal elements of 𝐼⁻¹(𝛼,𝛽) provide the asymptotic variances for the parameters 𝛼, 𝛽 and respectively. A two-sided 100(1−𝛾)% normal approximation confidence interval of 𝛼 can be obtained as 𝛼∓𝑍_𝛾/2var(𝛼). (13)

Similarly, a two-sided 100(1−𝛾)% normal approximation confidence interval of 𝛽 can be obtained as 𝛽∓𝑍_𝛾/2var𝛽, (14)

where 𝑍 is the standard normal variate (SNV).

3.2. Bootstrap Confidence Intervals

In this subsection, we have discussed another method for obtaining the confidence intervals proposed by [14]. They have developed the computer-based technique that can be routinely applied without any heavy theoretical consideration. The bootstrap method is very useful when an assumption regarding the normality is invalid. For bootstrap procedure, the computational algorithm is given as follows:

Step 1. Generate sample {𝑥₁,𝑥₂,…,𝑥_𝑛} of size 𝑛 form (1) by using inversion method. Then estimated distribution function is given by 𝐹(𝑥,Θ), where Θ={𝛼,𝛽;𝛼,𝛽>0}.

Step 2. Generate a bootstrap sample {𝑥^∗₁,𝑥^∗₂,…,𝑥^∗_𝑛} of size 𝑛 from 𝐹(𝑥,Θ). Obtain bootstrap estimates Θ^∗={𝛼^∗,𝛽^∗} of Θ={𝛼,𝛽} using bootstrap sample.

Step 3. Repeat Step 2, 𝐵-times. Obtain the bootstrap estimates {𝛼^∗₁,𝛼^∗₂,…,𝛼^∗_𝐵} and {𝛽^∗₁,𝛽^∗₂,…,𝛽^∗_𝐵}.

Step 4. Let {𝜃^∗₍₁₎,𝜃^∗₍₂₎,…,𝜃^∗_(𝐵)} be the ordered values of a sequence {𝜃^∗₁,𝜃^∗₂,…,𝜃^∗_𝐵} of a variable 𝜃. Then, the empirical distribution function (EDF) of {𝜃^∗₁,𝜃^∗₂,…,𝜃^∗_𝐵} is given by 𝐺(𝑡)={Numberof(𝜃^∗<𝑡)/𝐵}. The 100(1−𝛾)% boot-p confidence intervals for 𝜃 can be obtained by the following formula: (𝜃^∗_𝐵𝛾/2,𝜃^∗_{𝐵(1−𝛾/2)}). By using the above definition, the two-sided 100(1−𝛾)% boot-p confidence intervals for 𝛼 and 𝛽 are given by (𝛼^∗_𝐵𝛾/2,𝛼^∗_{𝐵(1−𝛾/2)}) and (𝛽^∗_𝐵𝛾/2,𝛽^∗_{𝐵(1−𝛾/2)}), respectively.

Step 5. The bootstrap measure of symmetry can be defined as 𝜃Shape=^{∗𝐵1−𝛾/2}−𝜃^∗𝜃^∗−𝜃^∗_𝐵𝛾/2. (15)

Using the above formula, we can also easily obtain the measure of symmetry (Shape) for 𝛼 and 𝛽. For standard normal approximate confidence intervals, the shape is always equals to one.

4. Bayes Estimation

In Bayesian scenario, we need to assume the prior distribution of the unknown model parameters to take into account uncertainty of the parameters. The prior densities for 𝛼 and 𝛽 are given as 𝑔₁(𝛼)∝𝛼^𝑏−1𝑒^−𝛼𝑎𝑔,𝑎,𝑏,𝛼>0,₂𝛽∝𝛽^𝑑−1𝑒^−𝛽𝑐,𝑐,𝑑,𝛽>0. (16)

Further, it is assumed that the parameters 𝛼 and 𝛽 are independent. Therefore, the joint prior of 𝛼 and 𝛽 is given by 𝑔𝛼,𝛽∝𝑔₁(𝛼)×𝑔₂𝛽, (17)

where 𝑎,𝑏,𝑐, and 𝑑 are the hyperparameters. Then, the joint posterior PDF of 𝛼 and 𝛽 can be readily defined as 𝜋=𝐿𝑥𝑔𝛼,𝛽∣𝑥∣𝛼,𝛽𝛼,𝛽∫^∞₀∫^∞₀𝐿𝑥∣𝛼,𝛽𝑔𝛼,𝛽𝑑𝛼𝑑𝛽. (18)

If 𝜆(𝛼,𝛽) is the function of 𝛼 and 𝛽, then the Bayes estimates of 𝜆(𝛼,𝛽) are given by 𝜆=∫𝛼,𝛽∣𝑥^∞0∫^∞₀𝜆𝛼,𝛽𝐿𝑥∣𝛼,𝛽𝑔𝛼,𝛽𝑑𝛼𝑑𝛽∫^∞₀∫^∞₀𝐿𝑥∣𝛼,𝛽𝑔𝛼,𝛽𝑑𝛼𝑑𝛽. (19)

The above expression cannot be obtained in nice closed form. The evaluation of the posterior mean of the parameters will be complicated and it will be the ratio of two intractable integrals. In such situations, Monte Carlo Markov chain (MCMC) method, namely, Gibbs sampling techniques can be effectively used. For implementing the Gibbs algorithm, the full conditional posterior densities of 𝛼 and 𝛽 are given by 𝜋₁𝛼∣𝛽,𝑥=𝑒^{−𝛼(∑𝑟𝑖=1𝑥(𝑖)𝑟+𝑎)}_𝑖=1𝛽𝛼+𝑥²_(𝑖)×exp−^𝑟_𝑖=1exp𝛼𝑥_(𝑖)−𝛽𝑥_(𝑖)×𝛼^𝑏−1exp−(𝑛−𝑟)exp𝛼𝑥_(𝑟)−𝛽𝑥_(𝑟),𝜋₂𝛽∣𝛼,𝑥=𝑒^{−𝛽(∑𝑟𝑖=1(1/𝑥(𝑖)𝑟+𝑐))}_𝑖=1𝛽𝛼+𝑥²_(𝑖)×exp−^𝑟_𝑖=1exp𝛼𝑥_(𝑖)−𝛽𝑥_(𝑖)×𝛽^𝑑−1exp−(𝑛−𝑟)exp𝛼𝑥_(𝑟)−𝛽𝑥_(𝑟). (20)

The simulation algorithm consists of the following steps.

Step 1. Start with 𝑗=1 and the initial values of {𝛼⁽⁰⁾,𝛽⁽⁰⁾}.

Step 2. Using the initial values {𝛼⁽⁰⁾,𝛽⁽⁰⁾}, generate candidate points {𝛼_𝑐^(𝑗),𝛽_𝑐^(𝑗)} from proposal densities {𝑞₁(𝛼^(𝑗),𝛼^(𝑗)),𝑞₂(𝛽^(𝑗),𝛽^(𝑗−1))}, where, 𝑞(𝛼^(𝑗),𝛼^(𝑗−1))=𝑞(𝛼^(𝑗)→𝛼^(𝑗−1)) is the probability of returning a value of 𝛼^(𝑗) given a previous value of 𝛼^(𝑗−1).

Step 3. Generate 𝑈 uniform variate on range 0-1; that is, 𝑢∼𝑈(0,1).

Step 4. Calculate the ratios at the candidate point 𝛼_𝑐^(𝑗)and previous point 𝛼^(𝑗)𝑅₁=𝜋₁𝛼_𝑐^(𝑗)∣𝛽^(𝑗−1),𝑥𝑞₁𝛼_𝑐^(𝑗),𝛼^(𝑗−1)𝜋₁𝛼^(𝑗−1)∣𝛽^(𝑗−1),𝑥𝑞₁𝛼^(𝑗−1),𝛼_𝑐^(𝑗). (21)

Step 5. If 𝑢≤min(1,𝑅₁), accept the candidate point with probability min(1,𝑅₁), that is, 𝛼^(𝑗)=𝛼_𝑐^(𝑗). Otherwise set 𝛼^(𝑗)=𝛼^(𝑗−1).

Step 6. Similarly from Step 4, the ratio 𝑅₂=𝜋₂𝛽_𝑐^(𝑗)∣𝛼^(𝑗),𝑥𝑞₂𝛽_𝑐^(𝑗),𝛽^(𝑗−1)𝜋₂𝛽^(𝑗−1)∣𝛼^(𝑗),𝑥𝑞₂𝛽^(𝑗−1),𝛽_𝑐^(𝑗). (22)

Step 7. If 𝑢≤min(1,𝑅₂), accept the candidate point with probability min(1,𝑅₂), that is, 𝛽^(𝑗)=𝛽_𝑐^(𝑗). Otherwise set 𝛽^(𝑗)=𝛽^(𝑗−1).

Step 8. Repeat Steps 2–7 for all 𝑗=1,2,…,𝑀 and obtain (𝛼₁,𝛽₁),(𝛼₂,𝛽₂),…,(𝛼_𝑀,𝛽_𝑀).

Note that if the candidate point is independent on previous point, that is, if 𝑞(𝑥→𝑦)=𝑞(𝑥), then the M-H algorithm is called independence M-H sampler. The acceptance function becomes 𝜋min1,¹𝛼_𝑐^(𝑗)∣𝛽^(𝑗−1),𝑥𝑞₁𝛼_𝑐^(𝑗)𝜋₁𝛼^(𝑗−1)∣𝛽^(𝑗−1),𝑥𝑞₁𝛼^(𝑗−1). (23)

The Bayes estimates under SELF of the parameters can be obtained as the mean of the generated samples from the posterior densities by using the algorithm discussed previously. The formulae are given by 𝛼=𝐸_𝜋≈1𝛼∣𝑥𝑀−𝑀₀^𝑀_𝑘=𝑀0+1𝛼_𝑘,𝜆=𝐸_𝜋≈1𝜆∣𝑥𝑀−𝑀₀^𝑀_𝑘=𝑀0+1𝜆_𝑘, (24)(25)

where 𝑀₀ is the burn-in-period of Markov Chain. The HPD credible intervals for 𝛼 and 𝜆 can be constructed by using the algorithm given in [15]. Let {(𝛼_(𝑖),𝛽_(𝑖));𝑖=1,2,…,𝑀} be the corresponding ordered MCMC sample of {(𝛼_𝑖,𝛽_𝑖),𝑖=1,2,…,𝑀}. Then construct all the 100(1−𝜓)% credible intervals of 𝛼 and 𝜆 as 𝛼_[1],𝛼_{[𝑀(1−𝜓)]},…,𝛼_[𝑀𝜓],𝛼_[𝑀],𝛽_[1],𝛽_{[𝑀(1−𝜓)]},…,𝛽_[𝑀𝜓],𝛽_[𝑀]. (26)

Here, [𝑥] denotes the largest integer less than or equal to 𝑥. Then, the HPD credible interval is that interval which has the shortest length.

5. One-Sample Prediction

In a Type-II censoring scheme, the life test consists only of few observed items (say 𝑟) out of all units/items (say 𝑛) under study due to time and cost constraints. In practice, the experimenter may be interested to know the life time of the (𝑛−𝑟) removed surviving units on the basis of informative units. In such a situation, one-sample prediction technique may be helpful to give an idea about the expected life of the removed units. Let 𝑥₍₁₎,𝑥₍₂₎,…,𝑥_(𝑟);(𝑟≤𝑛) the observed censored sample and 𝑦₍₁₎,𝑦₍₂₎,…,𝑦_{(𝑛−𝑟)} be the unobserved future ordered sample from the same population. Let 𝑌_(𝑠)=𝑋_(𝑟+𝑠) represent the failure lifetimes of the remaining surviving units. From [16], the conditional PDF of 𝑌_(𝑠) given 𝑋=𝑥 can be obtained as 𝑓𝑦_(𝑠)∣𝑥_(𝑟)=(𝑛−𝑟)!()!()!×𝑦𝑠−1𝑛−𝑟−𝑠1−𝐹_(𝑠)^{𝑛−𝑟−𝑠}𝐹𝑦_(𝑠)−𝐹𝑥_(𝑟)^𝑠−1𝑓𝑦_(𝑠)1−𝐹𝑥_(𝑟)^𝑛−𝑟. (27)

Putting (1) and (2) in (27), we get 𝑓𝑦_(𝑠)∣𝑥_(𝑟)=(𝑛−𝑟)!()!()!𝛽𝑠−1𝑛−𝑟−𝑠𝛼+𝑦²_(𝑠)×^𝑠−1∞_𝑗=0^𝑖_𝑖=0𝑙=0−1^(𝑖+𝑗)𝑗𝑖𝑙𝑖!𝑠−1𝑛−𝑟−𝑠+𝑗+1^𝑖×exp(𝑖+𝑙−1)𝛼𝑦_(𝑠)−𝛽𝑦_(𝑠)+𝑙𝛼𝑥_(𝑟)−𝛽𝑥_(𝑟). (28)

Then, the predictive posterior density of future observables under Type-II censoring scheme is given by 𝑓₁𝑦_(𝑠)∣𝑥=^∞₀^∞₀𝑓𝑦_(𝑠)∣𝛼,𝛽,𝑥𝜋𝛼,𝛽∣𝑥𝑑𝛼𝑑𝛽. (29)

Equation (29) cannot be evaluated analytically. Therefore, to obtain the consistent estimator for 𝑓₁(𝑦_(𝑠)∣𝑥), MCMC sample obtained through Gibbs algorithm is used. The consistent estimate of 𝑓₁(𝑦_(𝑠)∣𝑥) is obtained as 𝑓^∗₁𝑦_(𝑠)=1∣𝑥𝑀−𝑀₀^𝑀−𝑀0_𝑖=1𝑓𝑦_(𝑠)∣𝛼_𝑖,𝛽_𝑖,𝑥. (30)

To obtain the estimate of future sample, we used the M-H algorithm, to draw the sample from (29). Similarly, form (24), we can estimate the future observations under SELF as the mean of simulated sample drawn from (29). The survival function of future sample can be simply defined as 𝑆_𝑦(𝑠)(𝑇)=1−^𝑇_{𝑦(𝑠)=𝑥(𝑟)}𝑓₁𝑦_(𝑠)∣𝑥𝑑𝑦_(𝑠)=1−^𝑇_{𝑦(𝑠)=𝑥(𝑟)}^∞₀^∞₀𝑓𝑦_(𝑠)∣𝛼,𝛽,𝑥×𝜋𝛼,𝛽∣𝑥𝑑𝛼𝑑𝛽𝑑𝑦_(𝑠). (31)

We can also obtain the two-sided 100(1−𝛼)% prediction intervals (𝐿_𝑠,𝑈_𝑠) for 𝑦_(𝑠) by solving the following two equations: 𝑃𝑦_(𝑠)>𝑈_𝑠=𝛼∣𝑥2,𝑃𝑦_(𝑠)>𝐿_𝑠𝛼∣𝑥=1−2. (32)

Confidence intervals can be obtained by using any suitable iterative procedure as the above equations cannot be solved directly.

6. Two-Sample Prediction

In some situations, only lifetime model is given, and no priori information is available then assumes 𝑝(𝑦_(𝑘)∣𝛼,𝛽,𝑥)=𝑝(𝑦_(𝑘)∣𝛼,𝛽). This leads to the two-sample prediction problems. That is, the experimenters are interested in the 𝑘th failure time in a future sample of size 𝑁 following the same life time distribution. From [16], the PDF of 𝑘th order statistics is given by 𝑝𝑦_(𝑘)∣𝛼,𝛽=𝑁!()!()!×𝐹𝑦𝑘−1𝑁−𝑘_{(𝑘)𝑘−1}1−𝐹𝑦_(𝑘)^𝑁−𝑘𝑓𝑦_(𝑘). (33)

Putting (1) and (2) in (29), we get 𝑝𝑦_(𝑘)∣𝛼,𝛽=𝑁!()!()!𝛽𝑘−1𝑁−𝑘𝛼+𝑦²_(𝑘)×^𝑘−1∞_𝑗=0𝑖=0−1^𝑗𝑗𝑖!𝑠−1𝑁−𝑘+𝑗+1^𝑖×exp(𝑖+1)𝛼𝑦_(𝑘)−𝛽𝑦_(𝑘). (34)

The predictive posterior density of future observables under Type-II censoring scheme is given by 𝑝₁𝑦_(𝑘)∣𝑥=^∞₀^∞₀𝑝𝑦_(𝑘)∣𝛼,𝛽,𝑥𝜋𝛼,𝛽∣𝑥𝑑𝛼𝑑𝛽. (35)

Equation (35) cannot be evaluated analytically. Therefore, to obtain the consistent estimator for 𝑝₁(𝑦_(𝑘)∣𝑥), MCMC sample obtained through Gibbs algorithm is used. The consistent estimate of 𝑝₁(𝑦_(𝑘)∣𝑥) is given by 𝑝^∗₁𝑦_(𝑘)=1∣𝑥𝑀−𝑀₀^𝑀−𝑀0_𝑖=1𝑝𝑦_(𝑘)∣𝛼_𝑖,𝛽_𝑖. (36)

To obtain the estimate of future sample, again M-H algorithm is used to draw the sample from (35). Similarly, form (24), we can estimate the future observations under SELF as the mean of simulated sample drawn from (35). The survival function of future sample can be simply defined as 𝑆_𝑦(𝑘)(𝑇)=1−^𝑇_𝑦(𝑘)=0𝑝₁𝑦_(𝑘)∣𝑥𝑑𝑦_(𝑘)=1−^𝑇_𝑦(𝑘)=0^∞₀^∞₀𝑝𝑦_(𝑘)∣𝛼,𝛽,𝑥×𝜋𝛼,𝛽∣𝑥𝑑𝛼𝑑𝛽𝑑𝑦_(𝑘). (37)

We can also obtain the two sided 100(1−𝛼)% prediction intervals (𝐿_𝑘,𝑈_𝑘) for 𝑦_(𝑘) by solving the following two nonlinear equations: 𝑃𝑦_(𝑘)>𝑈_𝑘=𝛼∣𝑥2,𝑃𝑦_(𝑘)>𝐿_𝑘𝛼∣𝑥=1−2. (38)

Confidence intervals can be obtained by using any suitable iterative procedure as the above equations cannot be solved directly.

7. Simulation Study

This section consists of the simulation results to compare the performance of the classical and Bayesian estimation procedures under different Type-II censoring schemes and parameter combinations. The comparison between the MLEs and Bayes estimators of the model parameters made in terms of their mean square errors (MSEs). We have also compared the average lengths of the asymptotic confidence intervals, bootstrap intervals and HPD credible intervals. For this purpose, we generate the sample of sizes 𝑛=20 small, 30 medium, and 50 large from (1) for fixed values of 𝛼=2 and 𝛽=2. We have considered the different Type-II censoring schemes for each sample size so that the sample contains 100%, 80%, and 60% of the available information.

The choice of the hyperparameters is the main concerning issue in the Bayesian analysis. Reference [17] argues that when information is not in compact form it is better to perform the Bayesian analysis under the assumption of non-informative prior. If we take 𝑎=𝑏=𝑐=𝑑=0, then posterior becomes as obtained under Jeffery’s scale invariant prior. For the choice of hyper parameters under the subjectivism, we have taken prior means equals to the true values of the parameters with varying variances. The prior variance indicates the confidence of our prior guess. A large prior variance shows less confidence in prior guess and resulting prior distribution is relatively flat. On other hand, small prior variance indicates greater confidence in prior guess. In this study, we have taken prior variance equals to 1 (small) and 8 (large), and we call these Gamma-1 and Gamma-2 priors, respectively.

For obtaining the Bayes estimates, we generate posterior deviates for the parameters 𝛼 and 𝛽 using algorithm discussed in Section 4. First thousand MCMC iterations (Burn-in period) have been discarded from the generated sequence. We have also checked the convergence of the sequences of 𝛼 and 𝛽 for their stationary distributions through different starting values. It was observed that all the Markov chains reached to the stationary condition very quickly.

For the unknown model parameters, we have computed MLEs and Bayes estimates under informative and non-informative priors along with their asymptotic confidence/bootstrap/HPD intervals. We repeat the process 1000 times, and the average estimates with the corresponding mean square errors (MSEs) of the estimators, and average confidence/bootstrap/HPD intervals are recorded. The simulation results are summarized in Tables 1, 2, 3, 4, 5, 6, 7, and 8. All necessary computational algorithms are coded in R-environment [18] and codes are available with the authors.

Table 1 

Average estimates and MSEs (in the brackets) of the estimators of 𝛼 (in the first row of each cell) and 𝛽 (in the second row of each cell) under different Type-II censoring schemes for fixed values of 𝛼=2 and 𝛽=2.

Table 2 

Average estimates and MSEs (in the brackets) of the estimators of 𝛼 (in the first row of each cell) and 𝛽 (in the second row of each cell) for complete sample 𝑛=20 with varying parameters.

Table 3 

Average estimates and MSEs (in the brackets) of the estimators of 𝛼 (in the first row of each cell) and 𝛽 (in the second row of each cell) for fixed 𝑛=20, 𝑟=16 with varying parameters.

Table 4 

Average estimates and MSEs (in the brackets) of the estimators of 𝛼 (in the first row of each cell) and 𝛽 (in the second row of each cell) for fixed 𝑛=20, 𝑟=12 with varying parameters.

Table 5 

Average 95% confidence/HPD/bootstrap intervals (in brackets) along with their width and shape (in square brackets) for bootstrap method under different Type-II censoring schemes for fixed 𝛼=2 and 𝛽=2.

Table 6 

Average 95% confidence/HPD intervals (in brackets) and bootstrap method along with their width and shape (in square brackets) for bootstrap method for complete sample 𝑛=20 with varying parameters.

Table 7 

Average 95% confidence/HPD intervals (in brackets) and bootstrap method along with their width and shape (in square brackets) for bootstrap method for complete sample 𝑛=20, 𝑟=16 with varying parameters.

Table 8 

Average 95% confidence/HPD intervals (in brackets) and bootstrap method along with their width and shape (in square brackets) for bootstrap method for complete sample 𝑛=20, 𝑟=12 with varying parameters.

On the basis of the results summarised in Tables 1–8, some conclusions can be drawn which are stated as follows: (i)The MSE of all the estimators decreases as sample increases (i.e., as 𝑛 and 𝑟 increases) for fixed values of 𝛼 and 𝛽. (ii)The MSE of all the estimators increases with increasing the value of the parameters for any fixed values of 𝑛 and 𝑟. (iii)The MSE of the maximum likelihood and Bayes estimator of 𝛼 increases with increasing 𝛼 for given values of 𝛽, 𝑛 and 𝑟. (iv)The MSE of the maximum likelihood and Bayes estimator of 𝛽 increases with increasing 𝛽 for given values of 𝛼, 𝑛 and 𝑟. (v)The Bayes estimators have the smaller risks than the classical estimators for estimating the parameters in all the considered cases. Although the Bayes estimates obtained under Gamma-1 prior more efficient than those obtained under Jeffery’s and Gamma-2 priors. This indicates that the Bayesian procedure with accurate prior information provides more precise estimates. (vi)The width of the HPD credible intervals is smaller than the width of the asymptotic confidence/bootstrap intervals. (vii)In all the cases, the bootstrap procedure provides larger width of the confidence intervals for the parameters and the wide range of the confidence intervals helps to cover the asymmetry. (viii)The shape is greater than one in all the considered cases which indicates that the distribution of the maximum likelihood estimators is positively skewed and becomes more skewed with decreasing sample size 𝑛.

8. Real Data Analysis

In this section, we analysed the data set of time between failures of secondary reactor pumps. This data set has been originally discussed in [19]. The chance of the failure of the secondary reactor pump is of the increasing nature in early stage of the experiment and after that it decreases. It has been checked by [10] that flexible Weibull distribution is well fitted model to this data set. The times between failures of 23 secondary reactor pumps are as follows: 2.160, 0.150, 4.082, 0.746, 0.358, 0.199, 0.402, 0.101, 0.605, 0.954, 1.359, 0.273, 0.491, 3.465, 0.070, 6.560, 1.060, 0.062, 4.992, 0.614, 5.320, 0.347, and 1.921.

For analysing this data set under Type-II censoring scheme, we generate two artificial Type-II censored samples from this real data set by considering two different values of 𝑟(=18,15). In the real applications, we have nothing in our hand other than few observations following any distribution function. Let us sketch the likelihood profile with respect to the parameters. The log-likelihood function is plotted over the whole parameter space in Figure 1. The contour plots of the likelihood function are also plotted in Figure 2. The maximum likelihood estimates, Bayes estimates, and corresponding confidence/bootstrap/HPD intervals are presented in Table 9 for different values of 𝑛 and 𝑟. The simulation runs and the histogram plot of simulated 𝛼 and 𝛽 are plotted in Figure 3.

Table 9 

Estimates and confidence/HPD/bootstrap intervals for real data set.

Figure 1 

Likelihood profile with respect to parameters for real data set.

Figure 2 

Likelihood contour plot with respect to parameters for real data set.

(a) Trace plot for Simulated 𝛼

(b) Posterior density of 𝛼

(c) Trace plot for Simulated 𝛽

(d) Posterior density of 𝛽

Figure 3 

MCMC runs and posterior density plot of 𝛼 and 𝛽 for complete real data set.

The summary of one-sample predictive densities of future samples is presented in Table 10 for different values of 𝑠. Two-sample predictive density functions for different values of 𝑟 and 𝑘 are summarised in Table 11. From Tables 10 and 11, it can be observed that the standard error (SE) of future observables increases as the values of 𝑠 and 𝑘 increase for one-sample and two-sample prediction, respectively. Two-sample predictive density of the first future ordered sample is plotted in Figure 4. The density of 𝑦₍₃₎th future sample in case of one-sample prediction is plotted in Figure 5.

Table 10 

The summary of the one sample predictive densities for real data set when 𝑟=20.

Table 11 

The summary of the two sample predictive densities for real data set.

Figure 4 

Two-sample predictive density function when 𝑘=1 for complete real data set.

Figure 5 

One-sample predictive density function when 𝑆=3 and 𝑟=20 for real data set.

The expressions (31) and (37) do not seem to be possible to compute analytically. Therefore, we prosed to use the Monte Carlo technique to solve these two equations. To compute the integral part 𝐼₁=^𝑇_{𝑦=𝑥(𝑟)}^∞₀^∞₀𝑓𝜋𝐼𝑦∣𝛼,𝛽,𝑥𝛼,𝛽∣𝑥𝑑𝛼𝑑𝛽𝑑𝑦,₂=^𝑇_𝑦=0^∞₀^∞₀𝑝𝑦∣𝛼,𝛽,𝑥𝜋𝛼,𝛽∣𝑥𝑑𝛼𝑑𝛽𝑑𝑦 (39)

appearing in the expressions, we follow the following steps.

Step 1. The approximate value of 𝐼₁ can be obtained as 𝐼₁=^𝑇_{𝑦=𝑥(𝑟)}𝐸_𝜋𝑓≈𝑦∣𝛼,𝛽,𝑥𝑑𝑦𝑇−𝑥_(𝑟)1𝑀^𝑀_𝑖=1𝐸_𝜋𝑓𝑦_𝑖,𝑥∣𝛼,𝛽,𝑥𝑦∼𝑈_(𝑟),𝑇. (40)

Similarly, 𝐼₂1≈𝑇𝑀^𝑀_𝑖=1𝐸_𝜋𝑝𝑦_𝑖∣𝛼,𝛽,𝑥,𝑦∼𝑈(0,𝑇). (41)

𝐸_𝜋 denotes the expectation of the function with respect to the joint posterior pdf of 𝛼 and 𝜆.

Step 2. Simulate 𝛼 and 𝜆 from (19) by using the algorithm discussed in Section 2.

Step 3. Generate 𝑈₁ from uniform density of size 𝑀 on range (𝑥_(𝑟),𝑇). That is, 𝑈₁∼𝑈(𝑥_(𝑟),𝑇).

Step 4. Generate 𝑈₂ from uniform density of size 𝑀 on range (0,𝑇). That is, 𝑈₂∼𝑈(0,𝑇).

Step 5. By using 𝑈₁ and 𝑈₂, then calculate the approximate values of integrals as 𝐼₁≈𝑇−𝑥_(𝑟)1𝑀^𝑀_𝑖=1𝐸_𝜋𝑓𝑈_𝑖∶1,𝐼∣𝛼,𝛽,𝑥₂1≈𝑇𝑀^𝑀_𝑖=1𝐸_𝜋𝑝𝑈_𝑖∶2∣𝛼,𝛽,𝑥. (42)

By using the above-discussed algorithm, the survival functions (31) and (37) of future samples are plotted in Figures 6 and 7, respectively.

Figure 6 

One-sample predictive survival function for real data set when 𝑟=20.

Figure 7 

Two-sample predictive survival function for real data set.

9. Conclusion

In this paper, we have considered the problem of estimation and prediction for flexible Weibull distribution in the presence of Type-II censored sample. We have found that Bayesian procedure provides the precise estimates of the unknown parameters of flexible Weibull model with smaller mean square error. The width of HPD intervals is smaller than that of the asymptotic and the bootstrap confidence intervals. Prediction has been applied in medicine, engineering, business, and other areas, and the Bayesian approach using MCMC methods can be effectively used to solve the prediction problems. The methodology developed in this paper will be very useful to the researchers, engineers, and statisticians where such type of life test is needed and especially where the flexible Weibull distribution is used.

Acknowledgment

The authors would like to thank the editor and revivers for the valuable suggestions to improve this manuscript. The third author (Vikas Kumar Sharma) thanks to the University grant commission (UGC), New Delhi for financial assistance.