On Progressive-Stress ALT under Generalized Progressive Hybrid Censoring Scheme for Quasi Xgamma Distribution

1  Introduction

Advancements in science and technology have led to the development of highly durable and complex products, such as silicone seals, computers, missiles, and LEDs. Reliability is essential for maintaining product quality, prompting manufacturers to invest heavily in testing and design. However, for highly reliable products, their extended lifespans often result in few or no failures during standard testing under normal conditions. Traditional life-testing methods struggle to handle such scenarios.

In traditional life testing and reliability experiments, obtaining sufficient failure data can be challenging, especially for highly reliable products with long lifespans. Under normal conditions, tests conducted within limited timeframes often result in very few failures. To address this, accelerated life testing (ALT) is commonly used. ALT involves subjecting products to elevated stress levels, such as increased humidity, temperature, pressure, voltage, or vibration, to induce failures more quickly. The resulting failure data are then analyzed to estimate the product’s life characteristics under normal operating conditions.

ALT methods are broadly categorized into three main types based on how the stress is applied over time: constant-stress, step-stress, and progressive-stress models [1]. These categories reflect different strategies to accelerate failure in highly reliable products in order to estimate their lifetime distribution within practical timeframes and costs.

•   Constant-stress ALT: Applies a fixed, elevated stress level throughout the entire test period.

•   Step-stress ALT: Increases the stress in discrete steps at predetermined time intervals or after specific events.

•   Progressive-stress ALT: Increases the stress continuously (e.g., linearly) over time, which better mimics gradual degradation processes.

Fig. 1 illustrates the conceptual differences between these ALT strategies.

Figure 1: ALT types

Constant-stress ALT applies a fixed high stress until failure or a preset time, and has been widely used in reliability analysis. Several studies have proposed models and estimation techniques under this scheme, including works by Kumar et al. [2], Hakamipour [3], Balakrishnan et al. [4], Abd El-Raheem et al. [5], Sief et al. [6], El-Sherpieny et al. [7], and Abd El-Raheem et al. [8]. Step-stress ALT gradually increases stress at set intervals or after specific failures and is well established in reliability studies. Key contributions include works by Gouno et al. [9], Balakrishnan et al. [10], Xu et al. [11], Mohie El-Din [12], Xu et al. [13], Alotaibi et al. [14,15], and Hassan and Abdelghaffar [16], who developed various inference methods under different censoring and lifetime models.

Progressive-stress ALT elevates stress over time, unlike constant- and step-stress methods. The ramp-stress test, which increases tension linearly, is a common example. This method closely resembles real-world deterioration processes and is popular in dependability studies. For instance, Yin and Sheng [17] examined maximum likelihood estimation for exponential models under progressive stress. Bai and Cha [18] studied optimal test design and failure rate modeling under increasing load. Wang and Fei [19] studied statistical inference for Weibull distributions in modified failure rate models under increasing stress. Moreover, AL-Hussaini et al. [20] introduced Bayesian prediction ranges for the half-logistic distribution using Type-II censored data, demonstrating its versatility and effectiveness.

Censored data are commonly used in studies on reliability and life testing. Due to factors like preserving working experimental units for future use, reducing the overall time for the test, and financial limits, investigators have to gather data using censored samples. Time-censoring (Type-I) and failure-censoring (Type-II) strategies are the two widely used censoring strategies in life-testing and reliability studies (see, for additional details, the work by Bain and Engelhardt [21]). These methods are not flexible enough to allow units to be removed from the experiment at any point other than the terminal point.

The progressive Type-II censoring scheme (PTII) is a widely utilized method in reliability and survival analysis, often preferred over the standard Type-II censoring scheme due to its adaptability in practical applications such as industrial and medical studies. Unlike the conventional approach, PTII permits the removal of still-operational units during the experiment, rather than at its conclusion [22]. In PTII, suppose m out of n units are intended to fail during a life test. A predetermined sequence (R1,R2,…,Rm) dictates the number of units removed after each failure. When the first failure occurs at time X1:m:n, R1 remaining units are arbitrarily removed. Similarly, after the second failure at X2:m:n, R2 additional units are removed. This process continues until the m-th failure at Xm:m:n, where all Rm remaining units are discarded, concluding the test. While this method has been extensively studied, it has a notable drawback: if the units are highly reliable, the experiment duration may become excessively long.

To address this limitation, Kundu and Joarder [23] introduced the progressive hybrid censoring scheme (PHCS), where the test ends at min(T,Xm:m:n), with T being a predetermined maximum time. However, PHCS struggles when only a few failures occur before T, leading to limited or imprecise parameter estimation. To mitigate these challenges, Cho et al. [24] proposed the generalized progressive hybrid censoring scheme (GPHCS). This scheme ensures a minimum number of failures are observed, improving parameter estimation accuracy. Under GPHCS, both the number of failures and their corresponding lifetimes are predetermined, balancing test duration and cost. The experiment concludes at T∗=max{Xk:m:n,min{Xm:m:n,T}}, where k, m, and (R1,…,Rm) are predefined, see [23]. The observations fall into one of three categories:

•   Case I: X1:m:n,X2:m:n,…,Xk:m:n, if T<Xk:m:n<Xm:m:n,

•   Case II: X1:m:n,…,Xk:m:n,…,XD:m:n, if Xk:m:n<T<Xm:m:n,

•   Case III: X1:m:n,…,Xk:m:n,…,Xm:m:n, if Xk:m:n<Xm:m:n<T,

a schematic representation of the generalized progressive hybrid censoring scheme is presented in Fig. 2.

Figure 2: Diagram of GPHCS strategy

More paper used GPHCS for model a new dialog of censored sample as: Koley and Kundu [25] introduced GPHCS in the presence of competing risks models. Maswadah [26] improved maximum likelihood estimation of the shape-scale family based on GPHCS. Salem et al. [27] introduced joint Type-II of GPHCS. Abdelwahab et al. [28] discussed classical and Bayesian inference for lifetime distribution based on GPHCS. Mohie El-Din et al. [29] obtained predictions of lifetimes under GPHCS. Mahto et al. [30] developed a partially observed competing risk model under GPHSC. Zhu obtained reliability inference for the multi-component stress–strength model under GPHCS. Wang et al. [31] discussed the dependence of competing risks with partially observed failure causes from bivariate distribution under GPHCS. Çetinkaya [32] discussed inference of P (X > Y) under GPHCS. Shi et al. [33] introduced reliability analysis of the J-out-of-N system under GPHCS.

Research on progressive-stress ALT continues to expand, incorporating diverse statistical models and censoring schemes. Rong-hua and Heliang [19] applied it to Weibull distributions using tampered failure rate models. Abdel-Hamid and Al-Hussaini [34] combined progressive stress with finite mixture distributions, Bayesian analysis, and progressive censoring for various lifetime distributions. Abdel-Hamid and Al-Hussaini [35] also studied inference for Weibull distributions under PTII censoring. Recent works utilizing progressive-stress ALT with censoring schemes include studies by Mohie et al. [36], Zhuang et al. [37], Mahto et al. [38], Abushal and Abdel-Hamid [39], Ismail [40], Hussam et al. [41], and Alotaibi et al. [42].

The GPHC scheme has proven to be a valuable tool for analyzing censored data in ALTs, particularly in models like competing risks and step-stress models. However, a significant gap exists in the application of GPHC to Progressive Stress ALT models, which involve progressively increasing stress levels during the test. Most of the existing studies, such as those by Wang et al. [43] and Pandey et al. [44], focus on constant or step-stress models, without addressing the complexities introduced by progressive stress scenarios.

This paper extends the progressive-stress ALT framework by incorporating the QXG distribution—a flexible lifetime model suitable for diverse failure behaviors—under the GPHCS with binomial removal. To the best of our knowledge, this integration has not been addressed in the literature. The study contributes by:

•   Developing classical and Bayesian estimation procedures for the QXG parameters under PSALT with GPHCS.

•   Providing predictive survival probabilities under different censoring schemes.

•   Conducting a thorough simulation to evaluate the performance of the estimators.

•   Addressing a critical gap in the literature and paves the way for more realistic modeling of accelerated life tests.

•   Demonstrating the methodology with real-world datasets.

This combination of modeling flexibility, advanced censoring, and estimation robustness aims to provide a more realistic and practically applicable framework for analyzing accelerated life test data.

The remainder of this paper is organized as follows: Section 2 presents the fundamental assumptions and details of the Quasi Xgamma distribution under the progressive-stress ALT model. Sections 3 and 4 describe the estimation methods, including both maximum likelihood and Bayesian approaches, respectively. Section 5 provides a comprehensive simulation study to assess the performance of the proposed estimators. Section 6 demonstrates the methodology using real-world accelerated life test datasets. Section 7 concludes the paper with final remarks and potential directions for future research.

2  Test Assumptions

2.1 Quasi Xgamma Distribution

The quasi xgamma distribution, introduced as a generalization of the xgamma distribution, adds flexibility in modeling lifetime data. Developed by Sen and Chandra [45], this distribution includes an additional parameter, α, which, when combined with the primary scale parameter θ, allows for more accurate modeling of real-world data, especially where classical models fall short. The probability density function (PDF) of the quasi xgamma (QXG) distribution is given by:

f(x;α,θ)=θ(1+α)(α+(θx)22)e−θx;α>0,θ>0,x≥0.(1)

The cumulative distribution function (CDF) of the QXG distribution is given by:

F(x;α,θ)=1−(1+α+θx+(θx)22)(1+α)e−θx;x≥0.(2)

The Quasi Xgamma (QXG) distribution was chosen for its flexibility in modeling lifetime data, especially under progressive-stress ALT. It includes a shape parameter α that adjusts the distribution’s form, allowing it to generalize gamma and xgamma distributions as special cases [45]. Structurally, QXG is a mixture of exponential and gamma distributions, enhancing its ability to capture varied failure behaviors. Prior applications in survival analysis, such as bladder cancer data, have demonstrated its superior fit and predictive accuracy compared to classical models [45]. These features make QXG particularly suitable for reliability studies involving complex censoring and stress structures.

This model is quite versatile in nature and resembles probabilistic behavior in the gamma distribution and exponential distribution; see Sen and Chandra [45]. Recently, papers used QXG distribution as a baseline to obtain a new good distribution as: Sen et al. [46] discussed the QXG-geometric distribution with applications in medicine. Sen et al. [47] introduced QXG-Poisson distribution. Ahsan-ul-Haq et al. [48] discussed analysis, estimation, and practical implementations of the discrete power QXG distribution. Hassan et al. [49] obtained a new generalized QXG distribution applicable to survival times. Wani and Shafi [50] introduced the generalized Lindley-QXG distribution.

Assumption: Progressive-stress ALT

1.   Lifetimes of units follow QXG(α,θ).

2.   We have stress function S(t)=νt, ν>0.

3.   The S(t) depends upon t. Also, λ is a function of time t.

4.   We also have λ(t)=1θ1[S(t)]θ2, where θ1,θ2>0 and are to be estimated.

5.   Linear cumulative exposure model (LCEM) is considered for modeling the effect of stress change, for more details, see Nelson [1].

Under consideration of LCEM, the CDF under progressive-stress Si(t) is expressed as

Gi(t)=Fi(Δt),i=1,2,…,s,

where

Δt=∫0tdwλi(w)=θ1νiθ2tθ2 + 1θ2+1.

Then by setting the value of the scale parameter θ=1, we have the CDF for the i-th progressive stress given as:

Gi(t;Ω)=1−(1+α+θ1νiθ2tθ2 + 1θ2 + 1+12(θ1νiθ2tθ2 + 1θ2+1)2)(1+α)e−θ1νiθ2tθ2 + 1θ2 + 1;t>0,α,θ1,θ2>0,(3)

where Ω is vector of parameter as (α,θ1,θ2), and the corresponding PDF is given by:

gi(t;Ω)=θ1νiθ2tθ21+α(α+12(θ1νiθ2tθ2 + 1θ2+1)2)e−θ1νiθ2tθ2 + 1θ2 + 1.(4)

Algorithm 1 has been used To generate sample from generalized progressive hybrid censoring scheme.

3  Estimation Methods

In this section, we utilize the maximum likelihood estimation method based on PSALT to estimate the model parameters α, θ1, and θ2. A concise explanation of GPHCS under PSALT is provided. Consider s stress levels, where ni items are subjected to testing at each level, with i=1,2,…,s. Additionally, the observed items, denoted as mi, ki, and Ti, for i=1,2,…,s, follow predefined censoring schemes ℛi1,ℛi2,…,ℛimi at failure times ti1:mi:ni,ti2:mi:ni,…,timi:mi:ni, for i=1,2,…,s.

During the first failure at the ith stress level, ℛi1 surviving units are removed from the experiment. Similarly, at the second failure in the ith stress level, ℛi2 surviving units are withdrawn from the remaining (ni−2−ℛi1) units. Finally, at the mith failure, all the remaining units, calculated as ℛimi=ni−mi−∑j=1mi−1ℛij, are withdrawn, bringing the experiment to an end. Refer to Fig. 2 for illustration.

Let the observed GPHCS data under the progressive-stress level Si(t), for i=1,2,…,s be represented by tij. Here, i=1,2,…,s and j=1,2,…,mi. The corresponding likelihood function with vector parameter Ω is then expressed as follows:

L(Ω)=∏i=1s[Ci∏j=1Digi(tij;Ω)(1−Gi(tij;Ω))Rij(1−Gi(Ti;Ω))δRT∗],(5)

where Ci is constant and does not depend on parameters Ω, the RT∗ can be define as ni−mi−∑j=1DiRij.

Using Eqs. (3) and (4) in (5), the likelihood function of the QXG distribution under GPHCS within the framework of the PSALT model is derived as:

L(Ω)=∏i=1s[Ci∏j=1Diθ1νiθ2tijθ2(1+α)Rij + 1(α+12(θ1νiθ2tijθ2 + 1θ2+1)2)(1+α+θ1νiθ2tijθ2 + 1θ2+1+12(θ1νiθ2tijθ2 + 1θ2+1)2)Rij×e−Rijθ1νiθ2tijθ2 + 1θ2 + 1((1+α+θ1νiθ2Tθ2 + 1θ2 + 1+12(θ1νiθ2Tθ2 + 1θ2 + 1)2)(1+α)e−θ1νiθ2Tθ2+1θ2+1)δRT∗],(6)

hence, the log-likelihood function, excluding the constant term Ci, is expressed as:

ℓ(Ω)=∑i=1s∑j=1Di(ln⁡θ1+θ2ln⁡νi−(Rij+1)ln⁡(α+1))+θ2∑i=1s∑j=1Diln⁡tij+∑i=1s∑j=1Diln⁡(α+12(θ1νiθ2tijθ2 + 1θ2+1)2)+∑i=1s∑j=1DiRijln⁡(1+α+θ1νiθ2tijθ2 + 1θ2+1+12(θ1νiθ2tijθ2 + 1θ2+1)2)−∑i=1s∑j=1DiRijθ1νiθ2tijθ2+1θ2+1+δRT∗[ln⁡(1+α+θ1νiθ2Tθ2 + 1θ2+1+12(θ1νiθ2Tθ2 + 1θ2+1)2)−ln⁡(α+1)−θ1νiθ2Tθ2 + 1θ2+1].(7)

3.1 Newton–Raphson Algorithm

The likelihood equations are obtained by taking the first partial derivatives of the log-likelihood function in (7) with respect to α, θ1, and θ2, expressed as:

∂ℓ(Ω)∂α=∑i=1s∑j=1Di−(Rij+1)α+1+∑i=1s∑j=1Di1α+12(θ1νiθ2tijθ2 + 1θ2+1)2+∑i=1s∑j=1DiRij1+α+θ1νiθ2tijθ2 + 1θ2+1+12(θ1νiθ2tijθ2 + 1θ2+1)2+δRT∗(1+α+θ1νiθ2Tθ2 + 1θ2+1+12(θ1νiθ2Tθ2 + 1θ2+1)2)−δRT∗α+1,(8)

∂ℓ(Ω)∂θ1=∑i=1sDiθ1+∑i=1s∑j=1Diθ1(νiθ2tijθ2 + 1θ2+1)2α+12(θ1νiθ2tijθ2 + 1θ2+1)2+∑i=1s∑j=1DiRijνiθ2tijθ2 + 1θ2+1+θ1(νiθ2tijθ2 + 1θ2+1)21+α+θ1νiθ2tijθ2+1θ2 + 1+12(θ1νiθ2tijθ2 + 1θ2+1)2−∑i=1s∑j=1DiRijνiθ2tijθ2+1θ2+1+δRT∗[νiθ2Tθ2 + 1θ2+1+θ1(νiθ2Tθ2 + 1θ2+1)21+α+θ1νiθ2Tθ2 + 1θ2+1+12(θ1νiθ2Tθ2 + 1θ2+1)2−νiθ2Tθ2 + 1θ2+1],(9)

and

∂ℓ(Ω)∂θ2=∑i=1sDiln⁡νi+∑i=1s∑j=1Diln⁡tij+∑i=1s∑j=1Di(θ1νiθ2tijθ2 + 1θ2+1)2A(θ2,tij)α+12(θ1νiθ2tijθ2 + 1θ2+1)2+∑i=1s∑j=1DiRijθ1νiθ2tijθ2 + 1θ2+1A(θ2,tij)+(θ1νiθ2tijθ2 + 1θ2+1)2A(θ2,tij)1+α+θ1νiθ2tijθ2 + 1θ2+1+12(θ1νiθ2tijθ2 + 1θ2+1)2−∑i=1s∑j=1DiRijθ1νiθ2tijθ2 + 1θ2+1A(θ2,tij)+δRT∗θ1νiθ2Tθ2 + 1θ2+1A(θ2,T)+(θ1νiθ2Tθ2 + 1θ2+1)2A(θ2,T)1+α+θ1νiθ2Tθ2 + 1θ2+1+12(θ1νiθ2Tθ2 + 1θ2+1)2−δRT∗θ1νiθ2Tθ2 + 1θ2+1A(θ2,T),(10)

where A(θ2,t)=ln⁡νi+ln⁡tij−1θ2+1.

The maximum likelihood estimates of the unknown parameters α, θ1, and θ2 are determined by solving the three Eqs. (8)–(10) by setting them to zero. Since these equations are nonlinear, a numerical method is employed to compute the estimates such as the Newton–Raphson (NR) algorithm are employed. The NR algorithm begins with an initial estimate (α(0),θ1(0),θ2(0)) and iteratively improves the estimates through the following update formula:

Ω(k)=Ω(k−1)−[f′(Ω(k−1))]−1f(Ω(k−1)).

Here, Ω denotes the vector of unknown parameters, f(Ω) corresponds to the vector of first-order derivatives, and f′(Ω) represents the Hessian matrix of second-order derivatives:

Ω=[αθ1θ2],f(Ω)=[∂ℓ(Ω)∂α∂ℓ(Ω)∂θ1∂ℓ(Ω)∂θ2],f′(Ω)=[∂2ℓ(Ω)∂α2∂2ℓ(Ω)∂α∂θ1∂2ℓ(Ω)∂α∂θ2∂2ℓ(Ω)∂θ1∂α∂2ℓ(Ω)∂θ12∂2ℓ(Ω)∂θ1∂θ2∂2ℓ(Ω)∂θ2∂α∂2ℓ(Ω)∂θ2∂θ1∂2ℓ(Ω)∂θ22].

3.2 Asymptotic Confidence Interval

In this subsection, approximate confidence intervals for the parameters are derived using the asymptotic distribution of the MLEs of Ω=(α,θ1,θ2). The asymptotic distribution of the MLEs for (α,θ1,θ2) is expressed as:

((α^−α),(θ^1−θ1),(θ^2−θ2))∼N(0,f′(Ω)−1),

where f′(Ω)ij−1 represents the variance–covariance matrix of the parameters Ω. The elements of the 3×3 matrix f′(Ω)−1 denoted as f′(Ω)ij for i,j=1,2,3 for (α,θ1,θ2), can be approximated by f′(Ω^)ij. Accordingly, an approximate 100(1−p)% confidence interval for the parameter Ω is given by:

Ω^i±Z1−p/2f′(Ω)ii−1,i=1,2,3,

where Ω^1≡α^, Ω^2≡θ^1, and Ω^3≡θ^2.

The following section explores the incorporation of prior knowledge into the analysis to achieve improved results.

4  Bayesian Estimation

In this section, Bayes estimates for the model parameters α, θ1, and θ2 are derived under symmetric and asymmetric loss functions, specifically the squared error loss function and the entropy loss function, respectively. When all the model parameters are unknown, a joint conjugate prior may not exist. Therefore, piecewise independent priors are considered. It is assumed that the parameters α, θ1, and θ2 follow informative gamma priors with hyperparameters (pi,qi). Consequently, the priors for α, θ1, and θ2 are defined as follows:

π1(α)∝αp1−1e−αq1,π2(θ1)∝θ1p2−1e−θ1q2,π3(θ2)∝θ2p3−1e−θ2q3,

where α,θ1θ2>0,pi,qi>0;i=1,2,3.

The joint prior density function of α, θ1, and θ2, under the assumption that these three parameters are independent, is expressed as:

π(Ω)=αp1−1e−αq1θ1p2−1e−θ1q2θ2p3−1e−θ2q3.(11)

The mean and variance of the MLE for each of the jth samples are computed and equated to the mean and variance of the gamma prior distribution to determine the hyperparameters. The estimators α^, θ^1, and θ^2 for α, θ1, and θ2, respectively. For more information about this crtirea see [52]. Solving the mean and variance of the MLE and mean and variance of the gamma prior distribution yields the estimators for the hyper-parameters:

pi^=(1I∑j=1IΩi^j)/(1I−1∑j=1I(Ωi^j−1I∑j=1IΩi^j)2)

and

qi^=(1I∑j=1IΩi^j)2/(1I−1∑j=1I(Ωi^j−1I∑j=1IΩi^j)2).

For more information about gamma prior, hyperparameters, and posterior distribution, see [53–56].

By combining Eqs. (11) and (6) and applying Bayes’ theorem, the joint posterior distribution is obtained as:

Π(Ω|t)∝∏i=1s[∏j=1Diνiθ2tijθ2(1+α)Rij + 1(α+12(θ1νiθ2tijθ2 + 1θ2+1)2)(1+α+θ1νiθ2tijθ2 + 1θ2+1+12(θ1νiθ2tijθ2 + 1θ2+1)2)Rij×αp1−1e−αq1θ1p2+∑i=1sDi−1e−θ1q2θ2p3−1e−θ2q3e−Rijθ1νiθ2tijθ2 + 1θ2+1×((1+α+θ1νiθ2Tθ2 + 1θ2+1+12(θ1νiθ2Tθ2 + 1θ2+1)2)(1+α)e−θ1νiθ2Tθ2 + 1θ2 + 1)δRT∗].(12)

Deriving analytical Bayes estimators for the three parameters α, θ1, and θ2 are challenging, necessitating the use of numerical techniques. The next section explores the MCMC method for obtaining the required estimates.

4.1 MCMC Method

The marginal posterior distribution of a parameter is derived by integrating the joint posterior distribution over the other parameter. Consequently, the marginal posterior probability density functions for α, θ1, and θ2 are expressed as follows:

Π(α|θ1,θ2,t)∝∏i=1s[∏j=1Di1(1+α)Rij + 1(α+12(θ1νiθ2tijθ2 + 1θ2+1)2)(1+α+θ1νiθ2tijθ2 + 1θ2+1+12(θ1νiθ2tijθ2 + 1θ2+1)2)Rij×αp1−1e−αq1((1+α+θ1νiθ2Tθ2 + 1θ2+1+12(θ1νiθ2Tθ2 + 1θ2+1)2)(1+α)e−θ1νiθ2Tθ2 + 1θ2+1)δRT∗],(13)

Π(θ1|α,θ2,t)∝∏i=1sθ1p2 + ∑i=1sDi−1[∏j=1Di(α+12(θ1νiθ2tijθ2+1θ2+1)2)(1+α+θ1νiθ2tijθ2 + 1θ2+1+12(θ1νiθ2tijθ2 + 1θ2+1)2)Rij×e−θ1(q2 + Rijνiθ2tijθ2 + 1θ2+1)((1+α+θ1νiθ2Tθ2 + 1θ2+1+12(θ1νiθ2Tθ2+1θ2+1)2)(1+α)e−θ1νiθ2Tθ2 + 1θ2 + 1)δRT∗],(14)

and

Π(θ2|α,θ1,t)∝θ2p3−1∏i=1s[∏j=1Diνiθ2tijθ2(α+12(θ1νiθ2tijθ2 + 1θ2+1)2)(1+α+θ1νiθ2tijθ2 + 1θ2+1+12(θ1νiθ2tijθ2 + 1θ2+1)2)Rij× e−θ2q3e−Rijθ1νiθ2tijθ2 + 1θ2 + 1((1+α+θ1νiθ2Tθ2 + 1θ2+1+12(θ1νiθ2Tθ2 + 1θ2+1)2)(1+α)e−θ1νiθ2Tθ2 + 1θ2 + 1)δRT∗].(15)

The Metropolis-Hastings (MH) algorithm is the preferred method for generating posterior samples to compute the desired Bayes estimates. This is because obtaining the conditional posterior distributions of the parameters α, θ1, and θ2 in the form of well-known distributions is challenging, making direct sampling from these distributions difficult. To address this, the conditional posterior distributions are approximated using well-known distributions. The following steps of the MH algorithm are then used to generate random samples from these approximated conditional distributions:

Steps of the Metropolis-Hastings algorithm:

1.    Initialize the parameter values (α,θ1,θ2) as (α0,θ10,θ20).

2.    Set the iteration index j=1.

3.    Generate proposed values for the parameters as:

α∼N(αj,[f′(Ω)]11−1),θ1∼N(θ1j,[f′(Ω)]22−1),θ2∼N(θ2j,[f′(Ω)]33−1),

where [f′(Ω)]−1 represents the variance-covariance matrix.

4.    Calculate the acceptance ratio:

P=π(αj,θ1j,θ2j∣x)π(αj−1,θ1j−1,θ2j−1∣x).

5.    Accept the proposed values (α<j>,θ1<j>,θ2<j>) with probability min(1,P).

6.    Repeat steps (3) to (5) B times to generate B samples for the parameters (α,θ1,θ2).

As described in [57], the Bayes estimators for the parameters α, θ1, and θ2 under the general entropy loss function can be defined as:

α~=(1B′∑j=1B′(α<j>)−c)−1/c,θ1~=(1B′∑j=1B′(θ1<j>)−c)−1/c,θ2~=(1B′∑j=1B′(θ2<j>)−c)−1/c.

We remove first B0 samples to discard the possible dependency on the initial guess, where B′=B−B0 and B0 is also known as burn-in sample.

The Bayes estimators for the parameters α, θ1, and θ2 under the squared error loss function can be defined as:

α~=1B′∑j=1B′(α<j>),θ1~=1B′∑j=1B′(θ1<j>),θ2~=1B′∑j=1B′(θ2<j>).

4.2 Highest Posterior Density (HPD) Credible Interval

In this subsection, the HPD credible interval (L,U) is determined for a random variable Ω~ by ensuring that it satisfies the following expression:

P(L≤Ω~≤U)=∫LUπ∗(Ω~∣t)dΩ~=1−α.

Given the difficulty of determining the interval (L,U) analytically, we utilize the posterior samples derived to calculate the desired HPD credible intervals using the approach outlined by Chen and Shao [58].

5  Simulation

To evaluate the performance of the methods presented in this paper, a Monte Carlo simulation is carried out. The MLEs and Bayes estimates are compared based on their relative absolute bias (RAB), and mean squared errors (MSEs). Additionally, the asymptotic confidence intervals (ACI) and HPD intervals are assessed in terms of their average interval lengths and coverage probabilities. The simulation study involves three different progressive-stress accelerated life tests (ALTs): the first is a simple ramp-stress ALT with two stress levels (s=2), and the second is a multiple ramp-stress ALT with four stress levels (s=4). Furthermore, two censoring schemes are considered by binomial random removal. Assume that the removal of an individual unit from the life test is independent of the other units, with each unit having the same probability p of being removed. Consequently, the number of units removed at each failure time follows a binomial distribution, expressed as:

P(Ri1=ri1)=(ni−miri1)pri1(1−p)ni−mi−ri1,0≤ri1≤ni(16)

and

P(Rij=rij∣Ri(j−1)=ri(j−1),…,Ri1=ri1)=(ni−mi−∑k=1(j−1)rikrij)prij(1−p)ni−mi−∑k=1jrik,(17)

where 0≤rij≤ni−mi−∑k=1j−1rik,j=2,3,…,mi−1, and i=1,2,...,s is the stress levels. Also supposing further that Ri is independent of tj for all stress level. Therefore, the joint likelihood function with binomial removal can be expressed as:

L(Ω,p)=L(Ω)P(R=r).(18)

Substituting (16) and (17) into (18), we get

P(R=r)=(n−m)!∏i=1m−1ri!(n−m−∑i=1m−1ri)!p∑i=1m−1ri(1−p)m−1∏i=1m−1(n−mri).

In all scenarios considered, we have assumed p=0.2, 0.9 without any loss of generality. Samples were generated for specified values of each stress level ni and mi using a binomial removal technique in conjunction with a predefined sampling scheme.

The true values of the parameters α, θ1, and θ2 are assumed as follows: Case 1: α=0.5,θ1=0.3, and θ2=0.4; Case 2: α=1.5,θ1=0.3, and θ2=0.8; Case 3: α=1.5,θ1=2, and θ2=0.8; Case 4:α=1.5,θ1=2, and θ2=2. For Bayesian estimation, it is typical to exclude a portion of the initial samples as a burn-in period to ensure that the Markov chain has converged to the correct stationary distribution. In our approach, we generate a total of M-H B=12,000 samples, with the first B0=2000 samples discarded as burn-in. The remaining samples can then be used to construct point estimates and interval estimates for α and λ. In the first ramp-stress ALT, the ramp stresses are set to v1=170 and v2=200. For the second ramp-stress ALT, the ramp stresses are v1=170, v2=200, v3=240 and v4=280. The simulation results, based on 5000 samples, are summarized in Tables 1–6.

The entire study is carried out for the first ramp-stress is Small sample: (n1=20), (n2=15), large sample: (n1=60), (n2=40); for the second ramp-stress is Small sample: (n1=20), (n2=15), (n3=12), (n4=10), large sample: (n1=40), (n2=30), (n3=20), (n4=10) and the ratio of censored sample mi=rm∗ni, and ki=rk∗ni. The resulting MLE and Bayesian estimates, along with the ACI and HPD interval lengths, are presented in Tables 3 and 4 for each of the two sets of (n,m) values. The RAB, and MSE are shown in brackets beneath the estimated values, while the coverage probabilities are displayed in brackets beneath the interval lengths. All results are based on 1000 simulation iterations.

The results shown in Tables 1–6 indicate that Bayes estimates of α,θ1,θ2, using an informative gamma prior, outperform MLE estimates in terms of RAB, and MSEs across all two censoring schemes and two different ramp stress ALTs. However, Bayes estimates of parameters α,θ1 and θ2 based on asymmetric entropy loss function generally perform better than MLE in all cases. Comparing the results for the two sample size sets (small and large), it is evident that the estimation methods point estimates for large sample size perform better across all three ramp-stress ALTs.

Regarding the confidence intervals, Bayesian credible intervals for α,θ1 and θ, under ramp-stress ALTs with two or four stress levels, exhibit shorter lengths compared to asymptotic intervals. Similarly, Bayesian intervals for α,θ1 and θ2 outperform both boot-p and boot-t intervals. In terms of coverage probabilities (CP), MLE and Bayesian intervals for α,θ1 and θ2 consistently achieve probabilities above the nominal level of 95%.

For comparison, consider the censoring scheme using a binomial removal, where items are withdrawn immediately after the first failure, vs. predicted values for items withdrawn at the r-th failure. The predicted estimates under the first censoring (p = 0.2) scheme are smaller than those under the second censoring scheme (p = 0.9). In ramp-stress ALTs with two stress levels, higher stress levels correspond to smaller predicted estimates, while lower stress levels lead to shorter lifetime values, indicating that higher stress levels result in earlier failure. Similar trends are observed in ramp-stress ALTs with four stress levels.

6  Application

To illustrate the proposed methodologies, we analyze an accelerated life test data set from [59,60].

First data set: The data presented in Chapter 5 of [59] were obtained from ramp-voltage tests conducted on miniature light bulbs designed to operate at a stress (voltage) of 2 V. In these tests, 62 miniature light bulbs were tested at a ramp rate of 2.01 V/h, while 61 bulbs were tested at a ramp rate of 2.015 V/h. The lifetime data from these ramp-voltage tests are provided by [36].

Second data includes failure times for the time-dependent dielectric breakdown of metal-oxide-semiconductor integrated circuits. The test was conducted at three elevated temperatures: 170°C, 200°C, and 250°C. For this study, we focus on the data collected at 170°C and 200°C.

To evaluate the suitability of the data for a QXG distribution, we calculate AIC, BIC, the Kolmogorov-Smirnov distance (KSD), and the corresponding p-values for each stress level. These results, along with the each data set, are presented in Tables 7 and 8, respectively. The QXG distribution is contrasted with four rival models, including the log-logistic (LogLog) [38], exponentiated Weibull (EW) [14], extended exponential (ExEx) distribution [36], and the logistic exponential (LogEx) distribution [61]. These models are evaluated using statistical measures such as AIC, BIC, KSD, and p-values.

Tables 7 and 8 identify the best-fitting model as the one with the lowest AIC, BIC, and KSD values, combined with the highest p-values, indicating a superior fit to the QXG distribution.

Fig. 3 displays the likelihood profiles for the parameters α, θ1, and θ2 of the QXG model based on PSALT using data II. Each plot shows the log-likelihood (l(α,θ1,θ2)) as a function of one parameter, with the others fixed at their optimal values. The red dots represent the MLEs, with the profiles exhibiting distinct peaks, indicating well-defined and identifiable parameter estimates.

Figure 3: Likelihood profile for parameters of QXG based on PSALT with Data I

The results presented in Tables 9 and 10 summarize the MLE of parameters for various models under the PSALT framework, evaluated for two distinct datasets: Data I and Data II. The models compared include QXG, ExEx, LogLog, EW, and LogEx, with model performance assessed using the AIC and BIC. Lower values of AIC and BIC generally indicate better model fit for QXG under the PSALT model.

For Data I (Table 9), the QXG model emerges as the best-fitting model, with the lowest AIC (841.0324) and BIC (846.6460). The LogEx model performs competitively, but its AIC (844.3204) and BIC (852.0134) remain slightly higher than QXG. The ExEx model performs the worst, as evidenced by its significantly higher AIC (885.5127) and BIC (891.1263). The LogLog and EW models provide moderate fits, but neither outperforms QXG. Notably, the parameter estimates under QXG show relatively small standard errors, indicating precise and reliable estimates.

For Data II (Table 10), the QXG model again demonstrates superior performance, achieving the lowest AIC (550.5654) and BIC (555.2314). The ExEx and LogEx models follow closely, with AIC and BIC values that are only slightly higher, suggesting these models also offer reasonable fits for Data II. The EW model, despite its added parameter (δ), performs the worst, with the highest AIC (559.7851) and BIC (566.0065). As with Data I, the parameter estimates under QXG for Data II show small standard errors, indicating precision and stability in the estimation process.

Table 11 presents progressively censored data for a specific dataset of amp-voltage tests conducted on miniature light bulbs designed to operate at a stress. The table shows different censoring schemes defined by k1, k2, and censoring probability (p). For each scheme, the table lists the observed failure times (x1) and the censored observations (R, where R>0 indicates a censored observation). This data is likely used to study the effects of different censoring schemes on statistical analyzes, such as parameter estimation and survival analysis.

Figs. 4 and 5 show the MCMC diagnostic plots for the parameters α, θ1, and θ2 of the QXG model using PSALT under GPHCS conditions (P=0.2, k1=15, k2=18), and k1=25, k2=25) with data I. Each row corresponds to one parameter and includes a trace plot (left), running mean plot (middle), and posterior density plot (right). The trace plots indicate good mixing, as the chains fluctuate around stable means without visible trends, while the running mean plots demonstrate convergence by stabilizing over iterations. The posterior density plots show well-defined unimodal distributions, highlighting the reliability of the parameter estimates. These diagnostics confirm that the MCMC sampling converged effectively, providing robust posterior samples for α, θ1, and θ2.