Statistical inference of constant-stress partially accelerated life tests under type II generalized hybrid censored data from Burr III distribution

Notation

n

total number of experiment items, where n = n ₁ + n ₂

T ₁ & T ₂

termination time of the experiment, predetermined by the experimenter, where T ₁ < T ₂

y _r:n

the time of the rth failure observed from n items under GTI-HCS

y _s:n

the time of the sth failure observed from n items under GTII-HCS

r & s

the failure number obtained at time T ₁ and T ₂ under GTI-HCS and GTII-HCS, respectively

y _1i

the items allocated at use condition in CSPALT, under GTII-HCS, where, i = 1, …, n ₁

y _2i

the items allocated at accelerated condition, in CSPALT, under GTII-HCS, where i = 1, …, n ₂

r ₁

the failed items observed from n ₁ before the censoring time h _1, in case of used condition under GTII-HCS

r ₂

the failed items observed from n ₂ units before the censoring time h ₂, in case of accelerated condition, under GTII-HCS

d ₁₁ & d ₂₁

the failure number obtained at time T ₁₁ and T ₂₁, respectively, at used condition in CSPALT under GTII-HCS

d ₁₂ & d ₂₂

the failure number obtained at time T ₁₂ and T ₂₂, respectively, at accelerated condition in CSPALT under GTII-HCS

D ₁ & D ₂

the number of total failures in the experiment up to time h ₁ & h ₂, respectively, in CSPALT under GTII-HCS

T ₁₁ & T ₁₂

time points predetermined by the experimenter in used condition of CSPALT under GTII-HCS

T ₂₁ & T ₂₂

time points predetermined by the experimenter in accelerated condition of CSPALT under GTII-HCS

Abbreviations

Abias

absolute Bias

AD

Anderson-Darling

AIC

Akaike information criterion

Asy-CIs

asymptotic confidence intervals

BCIs

Bayesian credible intervals

BEs

Bayesian estimates

BIC

Bayesian information criterion

BIIID

Burr III distribution

CAIC

corrected AIC

CDF

cumulative distribution function

CIs

confidence intervals

CP

coverage probability

CSPALT

constant-stress partial accelerated life testing

CVM

Cramer-von Mises

G1 & G2

first group and second group, respectively

GTI−HCS

generalized type I – hybrid censoring sampling

GTII−HCS

generalized type II – hybrid censoring sampling

HCS

hybrid censored sampling

HF

hazard rate function

HQIC

Hannan–Quinn information criterion

K-S

Kolmogorov–Smirnov

KSD

K-S distance

KSPV

K-S p-value

LACI

length of asymptotic confidence intervals

LCCI

length of Bayesian credible intervals

LXF

linear exponential loss function

MCMC

Markov chain Monte Carlo

MH

Metropolis-Hastings

MLE

maximum likelihood estimate

PALT

partial accelerated life test

PP

probability-probability

QQ

quantile-quantile

StEr

standard error

TTT

total time on test

V-CM

variance–covariance matrix

1 Introduction

In investigations of reliability and life testing, censored data are crucial when it becomes impossible to obtain all the data needed for the experiment due to factors like expense or schedule. Type-I and Type-II censoring are two of the most common types of censoring. A hybrid censored sampling (HCS) approach, which combines Type-I and Type-II censoring schemes, was proposed by Epstein [1]. Suppose that Y _1:n < Y _2:n <…< Y _n:n denote the ordered failure times of the experimental units. In Type-I HCS (TI−HCS), the life-testing experiment ends at time T ₁ = min (Y _r:n , T), where Y _r:n represents the time of the rth failure out of n items, and T is the maximum time point for the test, r ∈ ( 1 , 2 , … , n ) , and T ∈ ( 0 , ∞ ) . The test continues until a predefined number of units fail, at which point the experiment is terminated. Since the maximum amount of time allowed for the test under this system is T, there is a chance that relatively few failures in this HCS strategy will happen before time T ₁. As a result of this, Childs et al. [2] suggested a new HCS that ensures a fixed number of failures, known as Type-II HCS (TII−HCS). In this situation, the life test ends at time T ₂ = max(Y _r:n , T), where r ∈ ( 1 , 2 , … , n ) , and T ∈ ( 0 , ∞ ) . Chandrasekar et al. [3] enhanced these censored sampling strategies by establishing two extensions of this kind, known as generalized TI−HCS (GTI−HCS) and generalized TII−HCS (GTII−HCS). In the context of GTI−HCS, s , r ∈ ( 1 , 2 , … , n ) and T ∈ ( 0 , ∞ ) are fixed so that s < r. The experiment ends at Y _s:n if the sth failure happens after time T. Thus, a minimal number s of failures is ensured in GTI−HCS. In GTII−HCS, one fixes r ∈ ( 1 , 2 , … , n ) and T 1 , T 2 ∈ ( 0 , ∞ ) so that T ₁ < T ₂. The experiment ends at time T ₁ if the rth failure happens before that time T ₁. If the rth failure happens between T ₁ and T ₂, the experiment is ended at Y _r:n . Lastly, the experiment ends at time T ₂ if the rth failure happens after time T ₂.

Rapid technological advancement, consumer demand for consistent products, and competitive markets have all placed pressure on producers to supply consistent, quality products. It is extremely difficult to observe the failure time for complicated, high-reliability devices during life testing, such as lasers, optical fibers, semiconductors, metal fatigue, electric cables, and insulating materials. Due to this, accelerated life testing (ALT) or partial ALT (PALT) is recommended in the industrial sector to gather the essential failure data quickly enough to determine the failure’s link to external stress factors. Test items in ALTs are only tested at accelerated circumstances, or pressures higher than normal, in order to cause an early failure. A physically acceptable statistical model then extrapolates the data gathered under such accelerated settings to predict the lifespan distribution under typical use conditions. This test might save a significant amount of money, time, labor, and resources. There are several ways to administer stress, such as step stress, progressive stress, and continuous stress (see [4]).

The PALT is the most appropriate test to use, especially when test items are administered under both normal and higher-than-normal stress levels [5,6,7,8]. It is also the most crucial technique to apply when calculating the acceleration factor and projecting the accelerated test results to the situation. PALT integrates both standard life and ALT. As a result, PALT makes sense for calculating the acceleration factor (β > 1), which is defined as “the ratio of the hazard rate at the accelerated condition to that at normal conditions.” PALT may be broadly classified into two categories: step-stress PALT and constant-stress PALT (CSPALT). In step-stress PALT, after a certain number of failures or at a predetermined period, the test conditions for the remaining objects in the experiment change from normal use to higher stress. In CSPALT, every group of test units is subjected to different usage and accelerated conditions, which is the primary focus of this work. CSPALTs were examined by several authors, such as [9,10] for the Pareto distribution, [11] inverse Weibull distribution, [12,13] Gompartz distribution, [14] Weibull distribution, [15] Kumaraswamy distribution, [16] Lomax distribution, [17] weighted Lomax distribution, and [18] linear exponential distribution. For more studies, refer to [19,20,21,22,23].

To fit different failure life time data, Burr [24] created 12 cumulative distribution functions (CDFs). A lot of attention has recently been paid to the Burr type-III distribution (BIIID) among others since it can accept various hazard lifetime data. It can also accurately approximate a wide range of popular distributions for fitting lifetime data, including the gamma, Weibull, and log-normal distributions. Many scientific domains have employed the BIIID. Examples of these include modeling forestry-related events [25], fracture roughness [26], actuarial literature [27], meteorological literature [28], reliability theory [29,30], and operational risk [31], among others. The probability density function (PDF) of the BIIID is given by:

where ϑ , λ > 0 , are the shape parameters. The CDF of the BIIID is

The hazard rate function (HF) of the BIIID is given by:

Numerous scholars have examined inferences for the parameters of the BIIID. A Bayesian estimator for the BIIID based on double censoring was examined by Abd-Elfattah and Alharbey [32]. Kim and Kim [33] discussed the Bayesian and classical estimation methodologies of the BIIID parameters based on dual generalized order statistics. The BIIID’s estimate and prediction issues were addressed using TII censored and progressive Type II hybrid censored samples by Altindag et al. [34] and Gamchi et al. [35], respectively,. The unified hybrid censored sample was used by Panahi [36] to develop the statistical inference of a BIIID. Hassan et al. [37] addressed the estimation of the lifetime performance index for BIIID based on progressive censoring scheme. Dutta and Kayal [38] provided estimation and prediction for the BIIID using unified progressive hybrid censoring. Bayesian and non-Bayesian inferences of the BIIID under joint progressive TII censoring were examined by Hassen et al. [39].

So far, there have been no articles specifically focused on estimating parameters of the BIIID in CSPALT with GTII-HCS data. This article aims to address this gap, given the flexibility and versatility of the GTII-HCS, along with the widespread application of BIIID in modeling lifetime data. Consequently, this paper’s primary goals are as follows:

1. Using the GTII-HCS, the classical (maximum likelihood) estimation method will be utilized to derive the point estimators and asymptotic confidence intervals (Asy-CIs) for the parameters and acceleration factor in CSPALT model.

2. The Bayes procedure will be applied to obtain Bayesian estimators and the Bayesian credible intervals (BCIs) of the BIIID parameters. By assuming independent gamma priors for the parameters, we can compute Bayes estimates under different loss functions.

3. The Markov chain Monte Carlo (MCMC) technique will be employed, implementing Gibbs sampling within the Metropolis Hastings (MH) framework, to generate samples from the posterior distributions.

4. A simulation study will be conducted to assess the efficiency of the estimators based on several precision measures.

5. To examine two real datasets to show the applicability of the proposed estimators.

The format of the paper is as follows: Section 2 presents the model description along with the assumptions made. In Section 3, we discuss the maximum likelihood (ML) and Asy-CI estimators for the model parameters and acceleration factor. Section 4 focuses on the Bayesian and BCI estimators for the model parameters and the acceleration factor. A simulation study to evaluate the performance of these estimates for the specified model is conducted in Section 5. The analysis of actual data is detailed in Section 6, while Section 7 concludes with some key findings.

2 Model description and assumption

This section gives a description of the model and its assumptions.

3 ML estimators

This section provides the ML technique to generate the point estimators of the unknown parameters λ , ϑ , and acceleration factor ρ based on the GTII−HCS with CSPALT. Furthermore, the Asy-CIs of the model parameters and acceleration factor are determined.

Suppose that y 1 i = ( y 11 < … < y 1 n 1 ) , i = 1 , … , n 1 be n ₁ GTII−HCS drawn from BIIID with PDF (1) and CDF (2) under normal usage conditions. Also, assume that y 2 i ≡ ( y 21 < … < y 2 n 2 ) , i = 1 , 2 , … , n 2 be n ₂ GTII−HCS drawn from BIIID with PDF(5) and CDF(4) under accelerated conditions. Hence, the likelihood function in CSPALT for the BIIID under the GTII-HCS is obtained by combining the likelihood functions L ₁ and L ₂ as follows:

where d ℓ , ℓ = 1 , 2 is the number of total failures in the experiment up to time h ℓ and their values are given by:

Then, inserting Eqs. (1), (2), (4), (5) in (6) gives:

where ζ ≡ ( λ , ϑ , ρ ) . The natural logarithms of likelihood function is

where, υ 1 ( y 2 i , λ , ϑ ) = [ 1 − ( 1 + y 2 i − λ ) − ϑ ] , υ 2 ( h 1 , λ , ϑ ) = [ 1 − ( 1 + h 1 − λ ) − ϑ ] , and υ 3 ( h 2 , λ , ϑ ) = [ 1 − ( 1 + h 2 − λ ) − ϑ ] .

The ML estimators λ ˆ , ϑ ˆ , and ρ ˆ , of λ , ϑ , and ρ are produced by simultaneously solving the following normal equations:

where

One may determine the ML estimator of acceleration factor ρ as a function of the parameters λ and ϑ by using (9) and for fixed value, as seen below.

Then inserting Eq. (10) in (7) and (8) and solving them numerically, we get the ML estimators λ ˆ and ϑ ˆ . It is important to note that these equations cannot be solved analytically. Some numerical methods, such as the Newton-Raphson technique, can be used to get the necessary estimate. The ML estimator ρ ˆ of ρ is then produced by inserting λ ˆ and ϑ ˆ in Eq. (10).

It might be more useful to find a range of values that, with a given probability, include the unknown parameters rather than point estimates for them. Interval estimates are the name given to these ranges. The Asy-CIs of the unknown parameters ζ = ( λ , ϑ , ρ ) T are constructed here by employing the asymptotic characteristics of the ML estimators. Using the large sample theory, we can determine that the asymptotic distribution of the ML estimators ζ ˆ = ( λ ˆ , ϑ ˆ , ρ ˆ ) T is a normal distribution with a variance–covariance matrix (VC-M) of I − 1 ( ζ ) and a mean of ζ . To estimate I − 1 ( ζ ) , where the asymptotic VC-M I − 1 ( ζ ˆ ) is employed, which may be generated by inverting the observed Fisher information matrix. The asymptotic VC-M in this instance is as follows:

In Appendix 1, the second partial derivatives are provided. Asymptotic normality of the ML estimators allows for the construction of the 100(1 − δ )% Asy-CIs of λ , ϑ , and ρ , respectively.

where z δ / 2 is the upper δ / 2 th percentile point of the standard normal distribution.

4 Bayesian estimation

In this part, the Bayes estimators for the unknown parameters λ , ϑ , and the acceleration factor ρ of the BIIID are presented. It is assumed that the independent parameter ρ has a truncated gamma distribution ρ − 1 ∼ gamma ( a 3 , b 3 ) , and that the random variables λ and ϑ have independent gamma prior distributions, λ ∼ Gamma ( a 1 , b 1 ) , and ϑ ∼ Gamma ( a 2 , b 2 ) . The joint prior of parameters and acceleration factor is as follows:

where, a 1 , a 2 , a 3 , b 1 , b 2 , and b ₃ are the hyper-parameters. Interestingly, the gamma prior density for the acceleration factor ρ was first used by DeGroot and Goel [40]. In the Bayesian analysis, independent gamma priors have been used because of the gamma distribution’s remarkable flexibility (see, for instance, Dey et al. [41] and Kundu and Howlader [42]). The likelihood function (6) and the joint prior distribution (11) are used to get the joint posterior distribution of λ , ϑ , and ρ using Bayes theorem.

where,

According to the Bayes method, selecting a loss function that matches each potential estimator will help determine which estimator is the best. In this case, estimates are obtained for two distinct categories of loss functions: symmetric and asymmetric loss functions. The first kind will be represented by the squared error loss function (SLF), and the second type will be illustrated by the use of the linear exponential loss function (LXF). The SLF is inappropriate when either an overestimation or an underestimating takes place. As an alternate choice to estimate the parameters in this situation, LXF can be used. When there is more substantial overestimation than underestimation, and vice versa, the LXF is helpful. Bayes estimates of the function U ( ζ ) = U ( λ , ϑ , ρ ) , under SLF and LXF, are given, respectively, by:

and

where the parameter q represents the sign that indicates the direction of asymmetry. The integral given by Eq. (12) cannot generally be obtained in a closed form. Here, the MCMC method is used to generate samples from the posterior distributions, and the Bayesian estimators for each parameter and acceleration factor are then computed. Making a decision from several MCMC plans that are offered might be difficult. Notable subclasses of MCMC methods include Gibbs sampling and the more general Metropolis inside Gibbs samplers. With the MCMC methodology, we can always achieve an appropriate parameter interval estimate since we generate the probability intervals based on the empirical posterior distribution, which gives it an edge over the ML method. This is commonly unavailable while using an ML estimator. The MCMC samples may actually be used to completely quantify the posterior uncertainty about the parameters λ , ϑ , and ρ . This also applies to any function that depends on the inputs by employing a kernel estimate of the posterior distribution. The following are λ , ϑ , and ρ ’s conditional posterior densities: The conditional posterior densities of λ , ϑ , and ρ are as follows:

and

The following MH-within-Gibbs sampling steps can be used to obtain samples of λ , ϑ , and ρ .

Step 1: Set the initial values λ ( 0 ) = λ ˆ , ϑ ( 0 ) = ϑ ˆ and ρ ( 0 ) = ρ ˆ .

Step 2: Set I = 1.

Step 3: Generate λ ∗ from a normal distribution with a mean ( λ ˆ ) and variance ( V λ ˆ ) , represented as N ( λ ˆ , V λ ˆ ) . Similarly, generate ϑ ∗ from the normal distribution with a mean ( ϑ ˆ ) and variance ( V θ ˆ ) , represented as N ( ϑ ˆ , V ϑ ˆ ) . Also, generate ρ ∗ from the normal distribution with a mean ( ρ ˆ ) and variance ( V ρ ˆ ) , represented as N ( ρ ˆ , V ρ ˆ ) .

Step 4: Obtain ℏ λ = min 1 , π ( λ ⁎ ∣ ϑ ( I − 1 ) , ρ ( I − 1 ) , y ̲ ) π ( λ ( I − 1 ) ∣ ϑ ( I − 1 ) , ρ ( I − 1 ) , y ̲ ) , ℏ ϑ = min 1 , π ( ϑ ⁎ ∣ λ ( I − 1 ) , ρ ( I − 1 ) , y ̲ ) π ( ϑ ( I − 1 ) ∣ λ ( I − 1 ) , ρ ( I − 1 ) , y ̲ ) , and ℏ ρ = min 1 , π ( ρ ⁎ ∣ λ ( I − 1 ) , ϑ ( I − 1 ) , y ̲ ) π ( ρ ( I − 1 ) ∣ λ ( I − 1 ) , ϑ ( I − 1 ) , y ̲ ) .

Step 5: Generate samples U j ; j = 1 , 2 , 3 from the uniform U(0, 1) distribution.

Step 6: If U 1 ≤ ℏ λ , U 2 ≤ ℏ ϑ , and U 3 ≤ ℏ ρ , then set λ ( I ) = λ ⁎ , ϑ ( I ) = ϑ ⁎ , ρ ( I ) = ρ ∗ ; otherwise λ ( I ) = λ ( I − 1 ) , ϑ ( I ) = ϑ ( I − 1 ) , and ρ ( I ) = ρ ( I − 1 ) , respectively.

Step 7: Set I = I + 1.

Step 8: Repeat steps 3–7 B times and obtain λ ( I ) , ϑ ( I ) , and ρ ( I ) , for I = 1, 2, …, B.

5 Simulation study

This section explores the model’s performance through simulations. First, we present an illustrative example using a simulated dataset. Subsequently, we conduct a more extensive Monte Carlo simulation study. All results are accompanied by detailed discussions to explain their significance. It is important to note that finding exact solutions for the model parameters (mentioned earlier) is mathematically challenging due to the system of nonlinear equations involved. In our case, with three parameters, we were dealing with three complex equations. Therefore, we employed numerical methods to estimate the desired parameters. Since we were able to derive the first and second derivatives of the objective functions (functions used for parameter estimation), we employed Newton’s method to find the ML estimates (MLEs). Convergence issues can arise in optimization processes if researchers don’t pick good starting values. To address this, we randomly generated initial values close to the actual model parameters. For brevity, we won’t delve into the details of Newton’s method, which is known for its efficient convergence properties.

This section focuses on a simulation study designed to compare the performance of different parameter estimation methods and their corresponding confidence intervals in Tables 1–6. We’ll be evaluating two types of estimates: MLEs and Bayesian estimates (BEs). The comparison will involve two key metrics:

• Mean squared error (MSE) and absolute bias (Abias): These metrics will be used to assess the accuracy of the estimates themselves. Lower MSE and Abias indicate estimates that are closer to the true parameter values on average.

• Average length and coverage probability (CP): We’ll compare the average lengths of two types of CIs: Asy-CIs (LACI) and BCIs (LCCI). Additionally, we’ll examine their CP. Ideally, CIs should have an average length that captures the true parameter value within a certain confidence level (e.g., 95%) most of the time.