Bayesian Inference with binomial removal for type II heavy-tailed Weibull distribution under unified progressive-hybrid censoring: simulation and application

1 Introduction

Scientific research is very important in our lives, especially in reliability analysis. So studying the lifetimes of items and components is very important these days. There are many challenges in studying item spans, including cost and time. With applications to COVID-19 mortality, engineering, business, medical, and data sets, this study advances statistical inference for a TII-HTW distribution under UPHCS. Censoring and life-testing are essential for time-to-event analysis when comprehensive follow-up is unattainable. They permit investigators to derive meaningful findings from incomplete data prevalent in health, engineering, and the social sciences, while maintaining the efficiency and cost-effectiveness of investigations to facilitate timely decision-making. There are many types of censoring.

There is a type I censoring scheme that depends only on time, ignoring whether failures occurred. Also, there is type II censoring, which depends only on the number of failures, regardless of time, and both schemes have many drawbacks; see Almetwally et al. [1], who studied the Marshall–Olkin alpha-power extended Weibull distribution under different censoring schemes. In both schemes, there were many disadvantages, such as the time being too long, no failures occurring before this time, or the number of failures taking a long time, which is very costly and time-consuming. So many authors used a new scheme called the hybrid scheme, which may be type I or type II. Generalized hybrid (GH) censoring strategies were initially introduced by Chandrasekar et al. [2] and have proven effective. In practical scenarios, each GH censoring technique presents distinct challenges. Type I GH censoring does not guarantee a minimum number of failures (e.g., m) and may lead to a situation with few or no observed failures prior to the specified threshold time.

To address these deficiencies, Balakrishnan et al. [3] introduced a unified hybrid (UH) censoring methodology that integrates the benefits of diverse GH censoring schemes, thus providing a more equitable and effective censoring strategy. Górny and Cramer [4] discussed this point in details. There are many drawbacks in previous distributions, such as the lack of a failure count before the predetermined time and a predetermined number of failures. The unified hybrid censoring scheme (UHCS) addresses these problems and is now very useful in reliability analysis because it ensures a specified number of failures. For more reading in detail about UHCS, see Mohammad et al. [5], who examined the unified hybrid censoring samples from the power Pratibha distribution, and Hasaballah et al. [6] used a unified hybrid censoring scheme on the generalized inverted exponential distribution. Moreover, Taher and Hussam [7] recently investigated the Perks distribution under unified hybrid censoring with applications to carbon-fiber gauge data.

In related work, Prakash et al. [8] studied progressive-stress accelerated life testing under UPHCS using classical, Bayesian, and machine-learning approaches, while Prakash et al. [9] developed inference for a reduced Type-I heavy-tailed Weibull distribution under progressive Type-II censoring. These studies are directly relevant to our work.

In this paper, we apply UPHCS with binomial removals to the TII-HTW distribution, and compare three estimation methods: MLE, MPS, and Bayesian inference. To the best of our knowledge, this is the first study to investigate the TII-HTW distribution under UPHCS with binomial removals. The MPS method is also applied for the first time within this framework. The proposed methodology is illustrated using four real datasets spanning engineering, business, medical, and public health contexts. Also, many authors discussed this scheme in detail; for example, Almetwally [10], who worked on a one-parameter lifetime model. Elshahhat and Seyam [11] worked on the very flexible Weibull model under UPHCS. They estimated its parameters and applied shock data to the model. Also, Yousef et al. [12] made a scientific contribution in the UHCS.

1.1 Unified progressively hybrid censoring scheme

In lifetime experiments, censoring schemes are essential for balancing duration and cost while ensuring sufficient failure observations. Let a life-test commence at time 0 with n identical units. Define two predetermined failure counts p and q satisfying 1 ≤ p < q ≤ n, along with two temporal thresholds T₁ and T₂ such that 0 < T₁ < T₂. A progressive removal scheme R=(R1,…,Rq) is specified, where Ri denotes the number of surviving units randomly withdrawn immediately after the i-th observed failure, subject to n=q+∑i=1qRi. Let Y_*i*:*q*:*n* represent the i-th progressively censored order statistic (the i-th observed failure time), and define the survival function S(·;*θ) = 1−F(·;θ*).

The experiment terminates at a stopping time T governed by the following four distinct scenarios:

T={Yq:q:n,if Yq:q:n<T1,T1,if Yp:q:n<T1≤Yq:q:n,Yp:q:n,if T1<Yp:q:n≤T2,T2,if T2<Yp:q:n.(1)

Suppose that ν are the items that failed till time T, where ν∈{q, d₁, p, d₂}, with d_*i* denoting the quantity of failures recorded over time T_*i* (i = 1, 2). So we can consider that the failed items till ending the experiment are R*, which is expressed as

R*={n-q-∑i=1q-1Ri,T=Yq:q:n,n-d1-∑i=1d1Ri,T=T1,n-p-∑i=1p-1Ri,T=Yp:q:n,n-d2-∑i=1d2Ri,T=T2.(2)

It is important to distinguish between termination at an observed failure time and termination at a fixed time. If the experiment terminates at Y_*q*:*q*:*n* or Y_*p*:*q*:*n*, no separate scheduled removal is made at the terminating failure; instead, all surviving units at that time are included in R^*. However, if termination occurs at T₁ or T₂, the removals associated with all failures observed before the fixed termination time are included in the progressive-removal product, and R^* denotes the remaining units withdrawn at the terminal time.

The removal probability p₁ is assumed to be a known and pre-specified design parameter. Conditional on the number of units available for removal, the progressive removals are generated sequentially according to a binomial mechanism. For a censoring design with m scheduled failures, the first m−1 removals are generated sequentially, while the final scheduled removal is defined by

Rm=n-m-∑i=1m-1Ri,(3)

so that

m+∑i=1mRi=n.

The removal mechanism is assumed to be independent of the lifetime variables and noninformative with respect to the parameters of the TII-HTW distribution. Thus, the inferential procedures are conditional on the observed removal vector R = (R₁, …, R_*m*). Accordingly, p₁ does not appear in the conditional likelihood of the lifetime parameters.

As we can see, this UPHCS overcomes the disadvantages of older schemes by guaranteeing failures before time ends, ensuring the experiment will not last forever.

To enhance the work's engagement, we employed maximum-likelihood and Bayesian inferential methodologies for the model parameters, along with the corresponding reliability metrics under the BR scheme. Bayesian computations are performed using Markov chain Monte Carlo (MCMC) with gamma priors, employing the Metropolis-Hastings (MH) algorithm. Also, we estimated both the confidence intervals and Bayesian credible intervals. Additionally, we made a Monte Carlo simulation encompassing several censoring configurations and evaluated the estimation methods it using numerous accuracy metrics, thus demonstrating the robustness of the Bayesian processes.

The subsequent sections of the paper are structured as follows. Section 2 delineates the rationale and merits of the distribution. Section 3 explores the classical point and interval estimation for the model outlined in UPHCS. Section 4 formulates Bayesian estimators for the model parameters and their interval estimation. Section 5 covers the Monte Carlo analysis utilizing simulated experiments and the final remarks regarding the simulations conducted with various censoring schemes under binomial removal. Section 6 examines the empirical data sets. Section 7 delineates the conclusion and the principal results of the paper.

2 Motivation for a heavy-tailed model

Heavy-tailed distributions have become essential tools in reliability engineering and survival analysis, as real-world phenomena frequently exhibit extreme observations whose probabilities decay more slowly than predicted by traditional light-tailed models [13]. Empirical evidence from the real datasets analyzed in this study reveals marked right-skewness, substantial dispersion, and extreme observations. These characteristics are manifested by a pronounced separation between the mean and the median, as well as sample kurtosis values exceeding those typically associated with light-tailed models. Such features are indicative of heavy-tailed behavior, where extreme outcomes occur with non-negligible probability.

In medical and reliability studies, rare but severe events, such as sudden increases in mortality rates during epidemic waves or unusually large failure counts in engineering systems, play a critical role in risk assessment and inference. Standard Weibull models, due to their light-tailed nature, may underestimate the probability of extreme observations, leading to inadequate tail inference, particularly in the presence of censoring. The real data analyses conducted in this paper provide empirical support for these considerations. For both the COVID-19 mortality data and the aviation reliability data, the TII-HTW model yields consistently better goodness-of-fit performance than the classical Weibull distribution, as evidenced by lower information criteria and improved tail accommodation. These findings confirm that the proposed TII-HTW distribution, which preserves the Weibull kernel while allowing controlled tail inflation and flexible hazard-rate behavior, offers a more appropriate and reliable modeling framework for censored medical and engineering data.

2.1 Research gap and main contribution

The gap is substantial, as the TII-HTW model has not been used before under this censoring scheme. We used a model that effectively captures intricate lifespan characteristics observed in contemporary engineering and physical systems. This model is very flexible and can fit many types of data, either medical or business, and no one has used it to estimate parameters via MPS with BR. So that motivated us to cover this gap. Based on empirical data across various scientific disciplines, we compared it with some competitors using four real data sets from many fields of life.

The novelty and main contributions of this paper are that we applied the TII-HTW distribution in a reliability experiment under UPHCS with binomial removals to medical and engineering data. We also discussed the MPS method for the first time and applied it at UPHCS in the simulation design. Moreover, we calculated Bayesian estimates and compared them with classical methods to determine which performed best. Finally, we conducted real data analysis and obtained important results, comparing our model with its rivals and demonstrating its superiority based on empirical data across various scientific disciplines using four real datasets from different fields of life.

2.2 The type II heavy-tailed Weibull distribution

In fields such as radiation biology and actuarial science, accurately modeling rare but catastrophic events is critical for robust decision-making [14]. Recently, novel construction methods based on transformer functions have led to the development of flexible heavy-tailed families. Specifically, Zhao et al. [15] introduced the type I heavy-tailed (TI-HT) family, while El-Saeed et al. [14] proposed the type II heavy-tailed (TII-HT) family as an extension offering greater flexibility in tail weight and hazard rate shapes. Building on this foundation, this study specializes the TII-HT framework to the Weibull baseline, yielding a new three-parameter distribution that retains the interpretability of the Weibull model while providing greater capacity to accommodate heavy-tailed phenomena across various censoring schemes.

Let the baseline cumulative distribution function (CDF) and probability density function (PDF) of the Weibull model be given by

F(x;κ,λ)=1-e-λxκ, f(x;κ,λ)=κλxκ-1e-λxκ, x>0,(4)

where κ > 0 is the shape parameter and λ > 0 is the scale parameter.

Following the TII-HT generator [14], the CDF of the TII-HTW distribution can be formulated and written as seen below

G(x;ψ,κ,λ)=1-(e-2λxκ1-(1-ψ)(1-e-λxκ))ψ, ψ,κ,λ>0,x>0,(5)

where ψ > 0 is a tail-weight parameter that controls the heaviness of the right tail. The corresponding PDF is

g(x;ψ,κ,λ)=ψκλxκ-1[2ψ+(1-ψ)e-λxκ]e-2ψλxκ[1-(1-ψ)(1-e-λxκ)]ψ+1.(6)

The survival function S(x) and hazard rate function h(x) are

S(x)=(e-2λxκ1-(1-ψ)(1-e-λxκ))ψ,(7)

h(x)=ψκλxκ-1[2ψ+(1-ψ)e-λxκ]1-(1-ψ)(1-e-λxκ).(8)

Remark. When ψ = 1, the TII-HTW distribution reduces to the standard two-parameter Weibull distribution with scale parameter λ^* = 2λ. For ψ < 1, the tail becomes heavier than that of the Weibull, making the model suitable for data with extreme outliers.

3 Classical method for estimating parameters

We used many types of classical methods of estimation, such as the Maximum Likelihood Estimation (MLEs) under UPHCS and the MPS. We will compare its results with Bayesian estimates. Also, we will estimate its confidence intervals. See Balakrishnan et al. [3] for more reading and details.

3.1 Maximum likelihood estimation under UPHCS

Suppose x = (x₁, x₂, …, x_ν) represents the observed failure times from the TII-HTW distribution with parameter vector *Θ* = (ψ, κ, λ)^⊤. Suppose that x = (x₁, …, x_ν) denotes the observed failure times and let

a={ν-1,if T=Yq:q:nor T=Yp:q:n,ν,if T=T1or T=T2.

Accordingly, the likelihood function under the UPHCS framework takes the standard progressive censoring can be written in the unified form

L(Θ∣x)∝[∏i=1νg(xi;Θ)][∏i=1aS(xi;Θ)Ri]S(T;Θ)R*,(9)

where g(·;*Θ*) and S(·;*Θ*) are the probability density function and survival function of the TII-HTW distribution, respectively, substituting the explicit forms from Equations 6, 7. Thus, terminal removals are included only once.

The likelihood in Equation 9 is conditional on the observed removal vector R. Since p₁ is fixed by the experimental design and the removal mechanism is assumed to be independent of the lifetimes, the probability mass function of R does not involve *Θ* = (ψ, κ, λ)^⊤ and can be absorbed into the proportionality constant. If p₁ were treated as unknown, the joint likelihood would be

Ljoint(Θ,p1∣x,R)=L(Θ∣x,R)P(R∣p1),(10)

where

P(R∣p1)=∏j=1m-1(n-m-∑k=1j-1RkRj)p1Rj(1-p1)n-m-∑k=1jRk.(11)

However, estimating p₁ is beyond the scope of the present study because p₁ is considered a known design quantity throughout the analysis.

More explicitly, the likelihood under the four termination cases is given by

L(Θ∣x)∝{[∏i=1qg(xi;Θ)][∏i=1q-1S(xi;Θ)Ri]S(xq;Θ)R,T=Yq:q:n<T1,[∏i=1d1g(xi;Θ)][∏i=1d1S(xi;Θ)Ri]S(T1;Θ)R,Yp:q:n<T1≤Yq:q:n,[∏i=1pg(xi;Θ)][∏i=1p-1S(xi;Θ)Ri]S(xp;Θ)R,T1<Yp:q:n≤T2,[∏i=1d2g(xi;Θ)][∏i=1d2S(xi;Θ)Ri]S(T2;Θ)R,T2<Yp:q:n.(12)

we obtain

L(Θ∣x)∝∏i=1νψκλxiκ-1[2ψ+(1-ψ)e-λxiκ]e-2ψλxiκ[1-(1-ψ)(1-e-λxiκ)]ψ+1 ∏i=1a[e-2λxiκ1-(1-ψ)(1-e-λxiκ)]ψRi ×[e-2λTκ1-(1-ψ)(1-e-λTκ)]ψR*.(13)

Thus, the log-likelihood function ℓ(*Θ) = logL(Θ∣x*) becomes

ℓ(Θ)=νlogψ+νlogκ+νlogλ+(κ-1)∑i=1νlogxi -2ψλ∑i=1νxiκ+∑i=1νlog[2ψ+(1-ψ)e-λxiκ] -(ψ+1)∑i=1νlog[1-(1-ψ)(1-e-λxiκ)] +ψ∑i=1aRilog(e-2λxiκ1-(1-ψ)(1-e-λxiκ)) +ψR*log(e-2λTκ1-(1-ψ)(1-e-λTκ)).(14)

The maximum likelihood estimates Θ^=(ψ^,κ^,λ^)⊤ are obtained by solving the system of score equations

∂ℓ∂ψ=0, ∂ℓ∂κ=0, ∂ℓ∂λ=0,(15)

Subject to positivity constraints. As closed-form solutions are unavailable, numerical optimization techniques (e.g., Newton-Raphson or quasi-Newton methods) must be employed. The observed Fisher information matrix I(Θ^), computed as the negative Hessian of ℓ(*Θ*) evaluated at the MLE, provides asymptotic standard errors and facilitates the construction of approximate confidence intervals. For further details on Fisher information and related uncertainty measures, see Husseiny et al. [16] and Mansour et al. [17].

3.2 Asymptotic interval estimators for maximum likelihood estimates

After obtaining the maximum likelihood estimates

Θ^=(ψ^,κ^,λ^)⊤,

Asymptotic confidence intervals can be constructed. So the MLE can be nearly normally distributed

Θ^~aN3(Θ,I-1(Θ^)),

such that I(Θ^) is the wellknown Fisher matrix evaluated at Θ^.

I(Θ^)=-[∂2ℓ(Θ)∂Θi∂Θj]Θ=Θ^, Θ1=ψ, Θ2=κ, Θ3=λ.

Accordingly, the approximate variance–covariance matrix of Θ^ is given by

Σ^=I-1(Θ^).

Thus, a two-sided 100(1−α)% asymptotic confidence interval for Θ_*j*, j = 1, 2, 3, is

Θ^j±zα/2Σ^jj,

where z_α/2 is the upper α/2 quantile of the standard normal distribution and Σ^jj is the jth diagonal element of Σ^.

Therefore, the ACIs for the model parameters are given by

ψ^±zα/2Σ^11, κ^±zα/2Σ^22, λ^±zα/2Σ^33.

For reliability measures such as S(t) and h(t), asymptotic intervals are obtained using the delta method. Let

ϕ(Θ)∈{S(t),h(t)}.

Then,

ϕ^=ϕ(Θ^)

is approximately normally distributed with variance

Var^(ϕ^)=∇ϕ(Θ^)⊤Σ^∇ϕ(Θ^),

where ∇ϕ(Θ^) is the gradient vector of ϕ(*Θ*) evaluated at Θ^.

Hence, the 100(1−α)% asymptotic confidence interval for ϕ is

ϕ^±zα/2∇ϕ(Θ^)⊤Σ^∇ϕ(Θ^).

3.3 Maximum product spacing estimation

The MPS method, originally introduced by Cheng and Amin [18], provides an effective alternative to maximum likelihood estimation, particularly in situations where the likelihood function is highly nonlinear or exhibits instability. In contrast to MLE, which relies on the maximization of the likelihood, the MPS approach is based on maximizing the geometric mean of spacings derived from the cumulative distribution function (CDF). This method has been shown to possess desirable properties such as consistency, asymptotic normality, and robustness, especially for heavy-tailed models and censored data.

Let x_1:ν:*n* < x_2:ν:*n* < ⋯ < x_ν:ν:*n* denote the ordered observed failure times under the UPHCS scheme. Define the spacings as let x_1:ν:*n*<x_2:ν:*n* < ⋯ < x_ν:ν:*n* denote the ordered observed failure times. The spacings are defined by

Di(Θ)=F(xi:ν:n;Θ)-F(xi-1:ν:n;Θ), i=1,…,ν,(16)

where x_0:ν:*n* = 0 and F(x_0:ν:*n*; *Θ*) = 0. The final spacing is

Dν+1(Θ)=1-F(xν:ν:n;Θ)=S(xν:ν:n;Θ).(17)

Consequently,

∑i=1ν+1Di(Θ)=1.

The product spacing function under the UPHCS framework takes the standard progressive censoring form

S(Θ∣x)∝[∏i=1ν+1Di(Θ)][∏i=1aS(xi;Θ)Ri]S(T;Θ)R*,(18)

where

a={ν-1,T=Yq:q:nor T=Yp:q:n,ν,T=T1or T=T2.

The final spacing D_ν+1(*Θ*) and the terminal censoring contribution S(T; *Θ)^R have different roles. The former completes the partition of the unit interval generated by the transformed failure times, whereas the latter represents the contribution of the units remaining at the termination time.

In particular, when termination occurs at T = T₁ or T = T₂, we have D_ν+1(*Θ*) = S(x_ν; *Θ*), while the contribution of the R^* units censored at the fixed time is S(T; *Θ)^R. Therefore, the final spacing must not be replaced by S(T; *Θ*), since x_ν and T generally differ under fixed-time termination.

For numerical optimization, the logarithm of the MPS objective is used:

Q(Θ)=∑i=1ν+1logDi(Θ)+∑i=1aRilogS(xi;Θ)+R*logS(T;Θ).(19)

The MPS estimator is then defined by

Θ^MPS=argmaxΘ>0Q(Θ).(20)

The MPS estimates of the parameter vector *Θ* = (ψ, κ, λ)^⊤ are obtained by maximizing (Equation 19), that is,

Θ^MPS=argmaxΘS(Θ∣x).

We cannot find a solution by solving these equations algebraically, so we will use a numerical solution.

3.4 Asymptotic intervals based on the MPS estimators

Under the standard regularity conditions for the MPS method, Θ^MPS is asymptotically normally distributed. The observed MPS information matrix is obtained from the negative Hessian of the log-MPS objective function in Equation 19:

IMPS(Θ^MPS)=-∂2Q(Θ)∂Θ∂Θ⊤|Θ=Θ^MPS.(21)

The estimated variance–covariance matrix is

Σ^MPS=IMPS-1(Θ^MPS).(22)

Therefore, an approximate 100(1−α)% confidence interval for Θ_*j*, where j = 1, 2, 3, is

Θ^j,MPS±zα/2[Σ^MPS]jj.(23)

For a reliability measure

ϕ(Θ)∈{S(t),h(t)},

the delta method gives

Var^[ϕ(Θ^MPS)]=∇ϕ(Θ^MPS)⊤Σ^MPS∇ϕ(Θ^MPS).(24)

Hence, the corresponding approximate confidence interval is

ϕ(Θ^MPS)±zα/2∇ϕ(Θ^MPS)⊤Σ^MPS∇ϕ(Θ^MPS).(25)

These intervals are used to calculate the average interval lengths and coverage probabilities of the MPS estimators reported in Tables 1–6.

Table 1. Simulation results for certain statistical measures for MLE, MPS, and Bayesian when large sample of Case I MLE, MPS, and Bayesian when a small sample (n = 40) of Case I.

Table 2. Simulation results for certain statistical measures for MLE, MPS, and Bayesian when large sample (n = 100) of Case I.

Table 3. Simulation results for certain statistical measures for MLE, MPS, and Bayesian when small sample (n = 40) of Case II.

Table 4. Simulation results for certain statistical measures for MLE, MPS, and Bayesian when large sample (n = 100) of Case II.

Table 5. Simulation results for certain statistical measures for MLE, MPS, and Bayesian when small sample (n = 40) of Case III.

Table 6. Simulation results for certain statistical measures for MLE, MPS, and Bayesian when large sample (n = 100) of Case III.

4 Bayesian inference

Bayesian estimation is a probabilistic framework for parameter estimation by integrating prior information with observed data. Unlike classical approaches that rely solely on the likelihood function, the Bayesian methodology allows for a full characterization of parameter uncertainty through the posterior distribution. This feature is particularly advantageous in reliability analysis under complex censoring schemes, where the available data may be limited or partially observed.

In the context of the TII-HTW distribution, the likelihood function exhibits a highly nonlinear structure due to the heavy-tailed behavior and the unified progressive hybrid censoring mechanism. As a result, classical estimators may suffer from instability or reduced efficiency, especially in small samples or under high censoring levels. Bayesian methods offer a robust alternative by stabilizing inference through prior incorporation and by enabling direct probabilistic interpretation of the model parameters. Recent work has explored Bayesian approaches for Weibull-type distributions with flexible tail behavior [19]; in this section, we adopt gamma priors for the TII-HTW parameters.

Let *Θ* = (ψ, κ, λ)^⊤ denote the unknown parameter vector of the model, where all parameters are strictly positive. To reflect this constraint and to maintain analytical flexibility, independent gamma prior distributions are assigned to each parameter. The gamma prior is particularly suitable for lifetime models as it is defined on the positive real line and can represent a wide range of prior beliefs through its shape and rate hyperparameters. Moreover, it facilitates computational implementation within MCMC frameworks.

Accordingly, the prior distributions are specified as

ψ~Gamma(a1,c1), κ~Gamma(a2,c2),λ~Gamma(a3,c3),(26)

leading to the joint prior density

π(ψ,κ,λ)∝ψa1-1κa2-1λa3-1exp[-(c1ψ+c2κ+c3λ)].(27)

The gamma prior distribution is widely adopted for positive-valued parameters due to its remarkable flexibility and computational tractability. By suitably selecting the hyperparameters, it can represent both informative and weakly informative prior beliefs, making it particularly attractive in reliability and lifetime modeling. The gamma priors employed here are informative simulation-based priors. Their hyperparameters are fixed before analyzing each main simulated sample and remain unchanged throughout all replications belonging to the same design configuration.

To elicit the hyperparameters of the independent gamma priors, an independent pilot simulation is conducted before the main Monte Carlo experiment. The pilot samples are generated using the same population parameters and censoring design as the corresponding main simulation scenario. However, they are generated independently and are not included in the replications used to evaluate the estimators.

Θ^j,pil(ℓ), ℓ=1,…,L0, j=1,2,3,

denote the MLE of the jth parameter obtained from the ℓth independent pilot sample, where Θ₁ = ψ, Θ₂ = κ, and Θ₃ = λ. The pilot mean and variance are calculated as

Θ^¯j,pil=1L0∑ℓ=1L0Θ^j,pil(ℓ),(28)

and

sj,pil2=1L0-1∑ℓ=1L0(Θ^j,pil(ℓ)-Θ^¯j,pil)2.(29)

By matching these pilot moments with the mean and variance of the gamma distribution, the hyperparameters are determined by

aj=Θ^¯j,pil2sj,pil2, cj=Θ^¯j,pilsj,pil2, j=1,2,3.(30)

Thus, the pilot samples used to construct the informative priors are independent of the main simulated samples used to evaluate the performance of the Bayesian estimators.

Here, L denotes the number of replications used to compute the empirical moments. This approach ensures that the prior distributions are centered around the classical estimates while preserving the variability observed in the data, thus yielding informative yet data-driven priors. In addition, its analytical form facilitates posterior construction, even in complex models; see Dey and Dey [20] and Almetwally [21] for further discussion.

By combining the likelihood function in Equation 13 with the joint prior in Equation 27, the resulting posterior density of the parameter vector *Θ* = (ψ, κ, λ)^⊤ is obtained as

π(ψ,κ,λ∣x) ∝ψa1-1κa2-1λa3-1{∏i=1νψκλxiκ-1[2ψ+(1-ψ)e-λxiκ]e-2ψλxiκ[1-(1-ψ)(1-e-λxiκ)]ψ+1 ×∏i=1a(e-2λxiκ1-(1-ψ)(1-e-λxiκ))ψRi} e-(c1ψ+c2κ+c3λ)(e-2λTκ1-(1-ψ)(1-e-λTκ))ψR.(31)

From Equation 31, the Bayes estimator of any arbitrary function g(*Θ*) under the squared-error loss (SEL) is given by its posterior expectation, i.e.,

g~=E(g(Θ)∣x)=∫∫∫g(Θ)π(Θ)L(Θ∣x)dΘ∫∫∫π(Θ)L(Θ∣x)dΘ,(32)

where *Θ* = (ψ, κ, λ)^⊤.

The integrals in Equation 32 are analytically intractable due to the highly nonlinear form of the posterior density. Consequently, Bayesian estimates are obtained numerically using Markov chain Monte Carlo (MCMC) methods. In particular, a Metropolis-Hastings (MH) algorithm is employed to generate dependent samples from the joint posterior distribution [22, 23].

To ensure positivity of the parameters and improve sampling efficiency, a random-walk MH algorithm is implemented on the logarithmic scale. The sampling scheme proceeds as follows:

• 1 Initialize *Θ*^[0] = (ψ^[0], κ^[0], λ^[0]) using the corresponding MLEs, and select proposal variances sψ2,sκ2,sλ2>0. Set ϱ = 1.
• 2 Propose new values on the log-scale: and set
• 3 Find this value where the product ψ^*κ^*λ^* arises from the Jacobian of the log-transformation.
• 4 Generate u~Unif(0, 1). If u ≤ q, set *Θ*^[ϱ] = *Θ; otherwise, set *Θ*^[ϱ] = *Θ*^[ϱ−1].
• 5 Compute for every x > 0, the reliability measures: where h(·) S(·) and represent both the hazard and survival functions of the TII-HTW model.
• 6 Assign ϱarrowϱ+1 and by making the same steps from 2-5 a number of 𝔇 iterations.

To reduce the effect of the starting real values, the first 𝔇₀ are removed from the calculations as burn-in. The remaining samples are used to approximate posterior expectations. So for the SEL, we can compute the Bayesian estimates φ(*Θ*) is given by

φ^BS=1𝔇∑ϱ=𝔇0+1𝔇φ(Θ[ϱ]), 𝔇=𝔇-𝔇0.

4.1 MCMC convergence diagnostics

To assess the convergence and mixing of the MCMC algorithm, multiple chains were generated from dispersed initial values. Convergence was evaluated using trace plots, autocorrelation function (ACF) plots, the Gelman–Rubin potential scale reduction factor R^, and the effective sample size (ESS). The Gelman–Rubin statistic is defined by

R^=V^W,(33)

where W denotes the within-chain variance and V^ is the estimated posterior variance based on the within- and between-chain variations.

Values of R^ close to one indicate convergence, while a sufficiently large ESS indicates adequate posterior sampling efficiency. In the present analysis, convergence was accepted when R^<1.05 for all parameters and when the retained samples provided an adequate ESS. The trace plots showed stable fluctuations around the posterior means without systematic trends, while the ACF plots demonstrated a satisfactory decline in serial dependence. Detailed diagnostic plots and numerical results are provided in the Supplementary material.

4.2 Bayesian interval estimators

To quantify the uncertainty associated with the Bayesian estimates, highest posterior density (HPD) credible intervals are constructed based on the generated MCMC samples. Unlike equal-tail intervals, HPD intervals possess the desirable property of having the shortest length among all intervals with the same posterior probability, making them particularly suitable for skewed posterior distributions.

Let {φ[ϱ]: ϱ=𝔇0+1,…,𝔇} denote the retained MCMC samples of a scalar function φ(*Θ*), where 𝔇₀ is the burn-in size and 𝔇⋆=𝔇-𝔇0. Arrange these samples in ascending order as

φ(1)≤φ(2)≤⋯≤φ(𝔇⋆).

Now as for the 100(1−α)% HPD credible interval, define

h=⌊(1-α)𝔇⋆⌋,

and determine the index

j*=argmin1≤j≤𝔇⋆-h{φ(j+h)-φ(j)}.

The corresponding HPD interval is then given by

(φ(j), φ(j+h)),

where ⌊·⌋ denotes the floor operator; see Chen and Shao [24] for further details.

5 Simulation results and interpretations

We examine the finite-sample efficacy of the suggested estimators for the model parameters *Θ* = (ψ, κ, λ)^⊤, along with the reliability metrics S(t) and h(t) under UPHCS with binomial removal.

Now we generate data using the following variables (n, m, r, T₁, T₂), all programing was done using R environment. Point estimation is carried out using maximum likelihood estimation (MLE), MPS, and Bayesian methods under informative priors. We estimated the asymptotic confidence intervals (ACIs) based on the delta method for MLEs, for Bayesian estimation, we also calculated HPD intervals.

The simulation settings are summarized as follows:

• We made 10, 000 Monte Carlo replications are conducted for each design configuration. Parameter settings are considered:
• Before conducting the main Monte Carlo experiment, an independent pilot simulation was performed for each population and censoring configuration. The MLEs obtained from the pilot samples were used to determine the gamma prior hyperparameters through the moment-matching procedure in Equation 30. The pilot samples were generated independently and were excluded from the 10, 000 main replications used to calculate the reported performance measures.
• For Bayesian estimation, three independently initialized MCMC chains were generated. Each chain consisted of 12, 000 iterations, with the first 2, 000 iterations discarded as burn-in. Convergence was assessed using trace plots, ACF plots, the Gelman–Rubin statistic, and the effective sample size.
• All numerical computations are performed in R. Numerical optimization for MLE is carried out using the Newton-Raphson technique, and we MCMC method for Bayesian estimation results.
• S(t) and h(t) are evaluated at the stopping time t=T (see Equation 1), which varies across replications. This approach verifies the estimator stability under different termination times.
• The sample sizes n∈{40, 100} are considered to represent small and large sample scenarios.
• The censoring time points T₁ and T₂ are determined through the quantile function of the TII-HTW distribution. Specifically, the quantile levels q₁∈{0.55, 0.70} and q₂∈{0.80, 0.90} are used to obtain T₁ = Q(q₁) and T₂ = Q(q₂), respectively. This setting helps preserve a balanced censoring design by avoiding excessively early or excessively late termination times.
• We used the following setup: for n = 40, m∈{28, 36} and r∈{23, 30}; for n = 100, m∈{70, 90} and r∈{60, 80}.
• Using a binomial scheme with parameter p₁∈{0.3, 0.8} we made the binomial removal. Specifically, the number of units removed at each failure time follows a binomial distribution. For the first removal, Conditionally, for j = 2, 3, …, m−1, where 0≤rj≤n-m-∑k=1j-1rk.
• For each target quantity φ∈{ψ, κ, λ, S(t), h(t)}, the performance of estimators is evaluated using the mean squared error (MSE). Interval estimators are compared in terms of average interval length (AIL) and coverage probability (CP).
• Unless otherwise stated, all interval estimates are constructed at the 95% confidence (or credibility) level.

The performance of the estimators is assessed using the mean squared error (MSE) and relative absolute bias (RAB). In addition, the lengths of the corresponding intervals estimators are reported, including the length of asymptotic confidence intervals (LACI), and the length of highest posterior density intervals (LCCI). The MSE is defined as

MSE(φ^)=1N∑j=1N(φ^(j)-φ)2,

where φ^(j) denotes the estimate (MLE or Bayesian) obtained from the jth Monte Carlo replication, φ is the true parameter value, and N = 10, 000 is the total number of replications.

The RAB for a generic estimator φ^ is defined as

RAB(φ^)=1N∑j=1N|φ^(j)-φ|φ,

For theses parameters or quantities φ∈{ψ, κ, λ, S(t), h(t)}.

To evaluate the interval estimators, the average confidence (or credible) length (ACL) and the CP at the 95% level are computed as

ACLφ95%=1N∑j=1N(Uφ^(j)-Lφ^(j)),CPφ95%=1N∑j=1N1{Lφ^(j)≤φ≤Uφ^(j)},

where 1{·} denotes the indicator function, and (Lφ^(j),Uφ^(j)) represent the lower and upper bounds of the 95% interval (asymptotic or Bayesian) for φ in the jth replication.

Maximum likelihood estimation is implemented using the maxLik package [25], while MCMC-based Bayesian inference is conducted and diagnosed using the coda package [26].

Results in the simulation tables from Tables 1–6 recoded results deduced from the replications, for the model parameters (ψ, κ, λ) as well as the reliability measures S(t) and h(t) under representative UPHCS designs.

5.1 Concluding results deduced from the simulation tables

• (1) Across all configurations, the estimators exhibit stable performance, With RAB and MSE decreasing as the available information increases (i.e., larger n and higher numbers of observed failures m and r), which is consistent with standard asymptotic behavior.
• (2) Under the independently elicited informative gamma priors considered in this study, the Bayesian estimators generally produce lower MSEs and shorter interval lengths than MLE and MPS. This advantage is conditional on the availability and relevance of the incorporated prior information. This highlights the benefit of incorporating prior information under progressive censoring.
• (3) The MPS method provides noticeable improvement over MLE in several settings, especially under moderate censoring, yielding smaller bias and more stable estimates. This confirms its effectiveness as a likelihood-based alternative in complex models.
• (4) Asymptotic confidence intervals (ACIs) based on MLE tend to be shorter in some cases; however, Bayesian HPD intervals are generally more accurate, achieving coverage probabilities closer to the nominal level with competitive interval lengths.
• (5) Increasing the censoring thresholds (T₁, T₂) improves parameter identifiability, leading to lower RAB and MSE and narrower intervals for (ψ, κ, λ) and S(t). The estimation of h(t) becomes more stable as the number of observed failures increases.
• (6) While the sample size gets larger, and by increasing (m, r) this enhances estimation accuracy across all methods, as reflected by decreasing RAB and MSE, shorter intervals, and CP values approaching the nominal level.
• (7) Higher removal probabilities p₁ reduce the effective sample size, leading to increased estimation error and wider intervals. Nevertheless, Bayesian estimation remains more robust compared to MLE and MPS under such conditions.
• (8) Overall, Bayesian estimation with informative priors is recommended for the UPHCS framework, followed by MPS as a competitive alternative to MLE. The results suggest that designs with larger n, m, and r and moderate removal probabilities provide the most reliable inference.

5.2 Simulation tables

Convergence was considered satisfactory when R^<1.05, the effective sample size exceeded 1,000, the acceptance rate was between 20% and 45%, and the absolute Geweke statistic was below 1.96. Across all investigated settings, these criteria were satisfied, and the trace and ACF plots showed stable mixing without noticeable trends or persistent autocorrelation.

6 Real data analysis

We evaluate the performance of the TII-HTW distribution using four real-world datasets spanning engineering, business, medical (diabetic nephropathy), and public health (COVID-19) contexts. Three of the four datasets air conditioning failures, pizza franchise costs, and COVID-19 mortality rates exhibit pronounced skewness and extreme values, making them suitable benchmarks for heavy-tailed modeling. The diabetic nephropathy dataset, by contrast, exhibits near-zero skewness and light-tailed behavior, providing a useful point of comparison. Summary statistics for all four datasets are presented in Table 7.

Table 7. Summary statistics for Data-I, Data-II, Data-III, and Data-IV.

It is important to note that all four datasets are originally complete samples. The UPHCS scheme with binomial removals is artificially imposed on these datasets solely for the purpose of illustrating the proposed methodology under a censored data framework.

For each dataset, model parameters are estimated via maximum likelihood estimation (MLE). Let ℓ^ denote the MLE, k indicates the parameters, and n represents the size of the sample under test. The following information criteria are computed: Akaike information criterion (AK), AIC with a correction (AKK), Hannan-Quinn information criterion (HQ), Bayesian information criterion (BK), and Consistent Akaike Information Criterion (CK).

AK=2k-2ℓ^, AKK=AK+2k(k+1)n-k-1, BK=klogn-2ℓ^HQ=2klog(logn)-2ℓ^, CK=-2ℓ^+k(logn+1)

To further evaluate the fit, we compute the Anderson-Darling (A²), Cramér-von Mises (W²), and Kolmogorov-Smirnov (D) statistics, along with their corresponding p-values. Smaller values of these statistics indicate a better fit to the data. To demonstrate the superior performance of the proposed model, we compare it against several well-known competitor distributions: the Inverse Power Lomax (IPL) [27, 28], Three-Parameter Burr XII (TBXII) [29, 30], Power Lomax (PL) [28], Singh-Maddala (SM) [31], Inverted Exponentiated Lomax (IEL) [32], Weibull (W) [33], and Kies (K) [34], which serve as baseline models.

Although the Three-Parameter Burr XII (TBXII), Power Lomax (PL), and Singh–Maddala (SM) distributions are mathematically equivalent through one-to-one parameter transformations, they are retained separately in this study because they are widely used under different parameterizations and names across the reliability, survival analysis, and applied statistics literature. Reporting results for all three models facilitates comparison with previous studies employing different formulations. As expected, these models generally yield identical or nearly identical likelihood-based information criteria, while any minor numerical differences arise from optimization and floating-point computations rather than genuine differences in the underlying distributional family.

The analysis of real datasets is a common practice in reliability and survival analysis for validating models [35–37].

6.1 Data-I: air conditioning system failures (Boeing 720 fleet)

The first dataset (hereafter Data-I) consists of n = 213 successive failure times (in hours) observed in the air-conditioning systems of a fleet of Boeing 720 jet aircraft. These data were originally reported by Proschan [38] and have since been widely used as a benchmark in reliability analysis, accelerated life testing, and heavy-tailed modeling studies; see Irshad et al. [39], Sirajo et al. [40], and Rahman et al. [41].

6.2 Data-II: pizza franchise start-up costs

The second dataset (hereafter Data-II) comprises n = 36 observations of annual start-up costs (in thousands of dollars) for a random sample of pizza franchises [42]. The sample exhibits positive skew and leptokurtosis, consistent with heavy-tailed behavior and with empirical findings in business start-up cost data.

6.3 Data-III: type 2 diabetic nephropathy patients data

The third dataset (hereafter Data-III) comprises n = 19 observations recorded from Type 2 diabetic nephropathy patients. These data have been previously analyzed in the context of bivariate modeling and copula-based dependence structures; see, for example El-Sherpieny et al. [43]. Similar copula-based approaches have also been studied by Alawady et al. [44].

The sample exhibits positive skewness, indicating the presence of heavy-tailed behavior commonly observed in clinical biomarkers and survival times for chronic disease patients.

6.4 Data-IV: COVID-19 mortality rates (Canada)

The fourth dataset (hereafter Data-IV) comprises n = 224 observations of mortality rates recorded during the COVID-19 pandemic in Canada. These observations were originally studied by Ahmad et al. [45]. A detailed analysis of the same dataset has also been conducted in recent work that proposes heavy-tailed and flexible Weibull-type models (see, for example, Nadir et al. [46]). The sample exhibits positive skew and leptokurtosis, consistent with heavy-tailed behavior and with empirical findings in pandemic mortality data.

Figures 1–4 present the visual assessment of the fitted distributions for all four datasets. Each figure displays three panels arranged vertically: (a) histogram with fitted PDFs and empirical vs. fitted CDFs, (b) PP plots, and (c) QQ plots for all competing models.

Table 8 summarizes the UPHCS results under three removal probabilities: low (0.3), moderate (0.5), and high (0.8), for each of the four real datasets.

Table 8. Complete UPHCS results: censoring schemes and parameter estimates.

For the Boeing data (n = 213), the estimates of ψ range from 0.030 to 0.033, κ from 1.114 to 1.136, and λ from 0.302 to 0.308. The reliability Ŝ(t) is approximately 0.68, while the hazard ĥ(t) is about 0.031 across all removal levels, indicating stable estimates regardless of the censoring intensity.

For the Pizza data (n = 36), ψ^ is roughly 0.006, κ^ varies between 1.957 and 2.030, and λ^ ranges from 0.010 to 0.013. The reliability Ŝ(t) falls between 0.73 and 0.75, and the hazard ĥ(t) is approximately 0.013 to 0.014. These values suggest that about 74% of pizza franchise costs exceed $75, which is consistent with the right-skewed nature of the data.

For the Diabetic data (n = 19), ψ^ is around 0.005 to 0.006, κ^ ranges from 2.910 to 3.167, and λ^ lies between 0.012 and 0.021. The reliability Ŝ(t) is approximately 0.67 to 0.71, while the hazard ĥ(t) is about 0.08 to 0.10. The wider confidence intervals observed here are expected due to the small sample size, yet the estimates remain reasonably stable.

For the COVID-19 data (n = 224), ψ^ is approximately 0.006 for the low and moderate removal probabilities, and 2.489 for the high removal probability. κ^ ranges from 2.179 to 2.503, and λ^ from 0.025 to 21.082. The reliability Ŝ(t) is approximately 0.69 to 0.70, and the hazard ĥ(t) ranges from 0.71 to 0.72. These estimates show stable mortality reliability across the different removal schemes.

Overall, the 95% confidence intervals are reasonably narrow, indicating precise estimation. As expected, higher removal probabilities (0.8) tend to produce slightly wider intervals due to the reduced effective sample size, though the estimates remain stable overall.

Table 9 provides the full observed failure times X_(*i*) and the corresponding removal vectors R_*i* for all datasets, ensuring complete reproducibility of the analysis. As expected, higher removal probabilities lead to more removals at earlier failure times, which is consistent with the nature of progressive censoring. For the Pizza data (ν = 16), removals are concentrated at the smallest cost values; for the Diabetic data (ν = 11), censoring is relatively light; and for the COVID-19 data (ν = 90), removals are spread across the range of mortality values.

Table 9. Complete removal vectors and observed failure times.

Together, Tables 8, 9 offer a comprehensive view of the UPHCS analysis. The first table summarizes the parameter estimates and reliability measures, while the second provides the detailed censoring information needed for replication. The consistency of the estimates across different removal probabilities supports the robustness of the TII-HTW model. The wider confidence intervals for the Diabetic dataset (n = 19) reflect the small sample size, yet the estimates remain stable. For the Pizza data, Ŝ(t)≈0.73–0.75 indicates that about 74% of start-up costs exceed $75, which aligns with the right-skewed distribution of the data. For the COVID-19 data, Ŝ(t)≈0.69–0.70 shows stable mortality reliability across the different removal schemes.

7 Conclusions and important results

Under UPHCS, maximum likelihood, MPS, and Bayesian estimation procedures for the TII-HTW distribution were carefully developed, and interval estimation methods were investigated to improve inferential accuracy under restricted or gradually censored data. Comprehensive Monte Carlo simulations were done. Also, the reliability measures were computed using both classical and Bayesian estimators. We obtained some important results from the simulation tables and concluded that the Under the independently elicited informative gamma priors considered in this study, the Bayesian estimators generally produce lower MSEs and shorter interval lengths than MLE and MPS. This advantage is conditional on the availability and relevance of the incorporated prior information. The MPS method provides noticeable improvement over MLE in several settings. Regarding the results from real data analysis, we can conclude that our model outperforms all its competitors, as we can see in Tables 10–13, as the TII-HTW distribution has the highest P-value and lowest information records, proving its superiority among all its rivals.

Table 10. Distribution comparison results for Data I.

Table 11. Distribution comparison results for Data II.

Table 12. Distribution comparison results for Data III.

Table 13. Distribution comparison results for Data IV.

Generative AI statement

The author(s) declared that Generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fams.2026.1932119/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Reza Pakyari, Qatar University, Qatar

Reviewed by: Ziaurrahman Ramaki, Persian Gulf University, Iran

Aman Prakash, Indian Institute of Management Lucknow, India