Change point detection for a skew normal distribution based on the Q-function

1.   Introduction

Studies on lifetime testing focused on the specific causes of failure. However, recent lifetime measurement models have been built upon the concept of multiple competing factors leading to the failure of lifetime testing units. These contemporary models for estimating lifetimes consider various competing causes of failure that challenge the durability of lifetime testing units.

Competing risks (CR) analysis is commonly used in medical research, reliability engineering, and econometric sciences to investigate outcomes in which numerous alternative events may affect the path to the primary event. For example, in a lifetime experiment conducted by Boag ^[1], the cause of death was recorded as either breast cancer or another cancer type, demonstrating that multiple factors can contribute to a subject's death or failure. In the current context, individuals infected with COVID-19 may also face other health complications, potentially resulting in death due to one of these underlying conditions. These models are referred to as CR models. For a more comprehensive understanding of CR, refer to Pintilie^[2] and Chen et al. ^[3].

In reliability engineering, the failure times of components with numerous possible failure modes are analyzed (for example, a machine that can fail owing to wear, electrical problems, or other factors). In another application, using the brake system, Liu et al.^[4] where found the brake system of a racing car consists of a pedal, a balance lever, a main cylinder, and a brake which is considered as CR.

In econometrics or social sciences, let us say that you are researching how long people work in a certain industry. Some conflicting hazards in this situation could end the employment period, including voluntary resignation, such as when someone leaves to look for work elsewhere, retirement, illness, or disability, meaning that situation can be treated using CR.

CR models are statistical tools designed to analyze scenarios where multiple events can occur concurrently or sequentially, potentially influencing each other's probabilities. These models find widespread application in fields, like medicine, engineering, and economics, where various factors can contribute to an outcome. CR models and uncertainty are closely intertwined. While CR models focus on the simultaneous or sequential occurrence of multiple events, uncertainty refers to the lack of complete knowledge about future events. Uncertainty can significantly impact CR models, arising from factors such as insufficient or biased data, model assumptions that may not accurately reflect reality, and unforeseen events that disrupt expected patterns. For further details, please refer to ^[5,6,7].

In various industries and businesses, lifetime tests are commonly employed. Because these tests can be costly and time-consuming, statisticians devised alternative procedures known as censored samples, which allow experimenters to stop work before all units fail. As a result, things are often checked using various censoring systems, with types I and II being the most popular. In such cases, removing live units from the testing process is not an option, resulting in the design of a more adaptable and progressive type-II censoring scheme (PTIICS), which marks a significant advancement in the practice of censorship and practical application. To understand the PTIICS approach and implementation comprehensively, it has been examined by Balakrishnan ^[8], and Ng and Chan ^[9].

1.1.   Generalized progressive hybrid censoring

This subsection shows one type of censoring scheme called the generalized progressive hybrid censoring scheme (GPHCS), introduced by Cho et al. ^[10] to overcome the problems with type-I progressive hybrid censoring schemes. It has also been proposed by Kundu and Joarder^[11], a mixture of type-I HCS and PTIICS since HCS combines type-I and type-II censoring schemes.

This technique would prefer to allow the experiment to continue after a predetermined time if very few failures have been noticed up to that time. So we assumed a predetermined minimum of $k$ failures, and the experimenter would ideally like to observe $m$ failures but is willing to accept a bare minimum of $k (k < m)$ failures.

The GPHCS can be depicted as follows: If $n\in \mathbb{N}$ identical items $X_{1}, X_{2}, ..., X_{n}$ are placed on a test and $\tau \in (0, \infty)$ is a predetermined time, with integers $k$ and $m$ being predetermined such that $k$ is the minimum number of failing items allowed since $k, m \in\{1, 2, ..., n\}$ and $(R_{1}, R_{2}, ..., R_{m})$ are predetermined integers such that

then at the time of each failure $X_{i:m:n}$, $R_i$, $i = 1, 2, \dots, m$, the remaining units are randomly removed. This process continues until the termination time

when all remaining units are removed from the experiment. Under this setup, the experimenter would ideally like to observe $m$ failures but would be willing to accept a minimum of $k$ failures. The number of observed failures until time $\tau$ is denoted by $d$. We noted the observed data and survival items removed in the GPHCS categories, which have one of the three possible types:

From Eq (1.1), we make one of the following observations:

For Case I, $X_{d:m:n} < \tau < X_{d+1:m:n} < ... < X_{k:m:n}, \ R_{d+1} = ... = R_{k-1} = 0$ and

For Case II, $X_{d:m:n} < \tau < X_{d+1:m:n} < ... < X_{m:m:n}$ and $X_{d+1:m:n}, ..., X_{m; m:n}$ are not observed and are at time $\tau$, which observed some survival items $R_{\tau}^*$ since

For Case III,

is the number of remaining units left at time $X_{m:m:n}$.

Figure 1 shows a graphical description of the GPHCS.

1.2.   Related literature

A survey of the literature on lifetime distributions with CR under various censorship schemes is suggested in this subsection. As shown in Table 1, some researchers considered the estimation problem based on CR under censoring data.

1.3.   Work motivation

This study aims to analyze the generalized Lomax distribution (GLD) under the GPHCS scheme in the context of the CR model. The GPHCS method is increasingly popular due to its potential to reduce time and cost in life-testing experiments significantly. Moreover, the GLD's versatility has led to its wide application across various fields. Figure 2 presents a visual breakdown of the methodological framework employed in this study. The core objectives of this study are:

$(1)$ Analyze the GLD's parameters through CR based on the GPHCS by using maximum likelihood (ML) and Bayesian approaches.

$(2)$ Acquire the asymptotic two-sided confidence intervals (CIs) and two bootstrap CIs for GLD parameters.

$(3)$ Explain the Bayesian estimation of the GLD's parameters using an independent gamma prior distribution under three loss functions.

$(4)$ Use the Markov chain Monte Carlo (MCMC) algorithm, to generate the Bayesian credible intervals (BCIs) and the highest posterior density (HPD) intervals.

$(5)$ Evaluate the performance of the suggested estimators with different GPHC schemes through a simulation study and real data.

The remaining sections of this article are arranged as follows: Section 2 formulates the model. In Section 3, the maximum likelihood estimates (MLEs) for the parameters, the reliability characteristics, and the relative risks are obtained. Additionally, the Bayesian estimates (BEs) based on the linear exponential (LINEX), generalized entropy (GE), and squared error loss functions are presented. The MCMC technique is used to create the BCI and HPD intervals. Asymptotic CIs (ACIs), Bootstrap CIs, and other significant functions of the model parameters are given in Section 4. Section 5 shows the performance of the suggested estimates by comparison through mean squared error (MSE). Section 6 describes a practical application. The study's conclusion and future research are offered in Section 7. Finally, two algorithms used and tables of numerical work results in the present article are included in Appendixes A and B, respectively.

2.   Model description

Suppose that $n\in \mathbb{N}$ identical components $X_{1}, X_{2}, ..., X_{n}$ are put on a GPHC life test in the presence of CR. Each component is exposed to two risks. We have

where $X_{1i}$ and $X_{2i}$ are the latent failure times and the pairs $(X_{1i}, X_{2i})$ are assumed to be distributed independently and identically (i.i.d). Now, the consequential data in the presence of CR under the GPHC experiment is shown:

where $\delta_{i} = 1$ when $X_{1i} < X_{2i}$ (failure from Cause 1) or $\delta_{i} = 2$ when $X_{1i} > X_{2i}$ (failure from Cause 2); $\forall i = 1, 2, ..., n$.

The GLD, a widely used power-type heavy-tailed model, has been employed in numerous applied problems. Specifically, it was utilized to analyze business failure data. The GLD is a Lomax generalization proposed by Maurya et al. ^[25]. Researchers have employed various generalizations of the Lomax distribution in studies for modeling business records to reliability and lifetime testing, as shown in Hassan et al.^[26] (for more details about the study of GLD, see Alghamdi^[27]).

The GLD is derived by using the power transformation method of Box and Cox ^[28]. This distribution has three parameters: only one shape parameter $\alpha_j; j = 1, 2$, and the common other parameters $\beta, \gamma > 0$ are scale and shape, respectively. The following are the probability density function (PDF) and the cumulative distribution function (CDF) for the GLD:

respectively. The GLD's PDF is denoted by $GLD(\alpha_{j}, \beta, \gamma),$ where $j = 1, 2$ and if $\gamma$ is less than or more than 1, the PDFs may exhibit a decrease in value or a unimodal form. For $\beta = 1$ and $\gamma = 1$, it leads to Burr XII and Lomax distributions, respectively. In numerous cases, the GLD often proves to be more appropriate than the Burr XII, transmuted Burr III, and Burr III distributions, as shown by Maurya et al. ^[25].

Remark 1. If the latent failure times $X_{1i}$ and $X_{2i}$ are i.i.d random variables following GLD$(\alpha_1, \beta, \gamma)$ and GLD$(\alpha_2, \beta, \gamma)$, respectively, then the random variable

follows GLD$(\alpha_1+\alpha_2, \beta, \gamma)$, where $\beta, \gamma$ are the scale and shape parameters and $(\alpha_1+\alpha_2)$ is the shape parameter.

Applying Remark 1, the reliability and hazard function of the random variable $X_i$, referred to as $\Re(t)$ and $\hslash(t),$ respectively, are provided by

Remark 2. The relative risks due to Causes 1 and 2 are determined using the independence of the latent failure times, and they are as follows:

where $j = 1, 2.$

Now, based on the observed sample in Eq (2.1), the likelihood function formula for three different cases is defined in a unified expression:

where

for simplicity of notation,

is the survival function and $j_1$ and $j_2$ are the numbers of the observed failure times due to reasons 1 and 2. The total number of failures observed at experiment terminated time is denoted by $W^*$; i.e.,

Also, $W^* = k$ for Case I, or $W^* = d$ for Case II, $W^* = m$ for Case III. The indicator function is denoted $D^*$; i.e., $D^* = 1$ when $W^* = d;$ otherwise, $D^* = 0.$

3.   Estimation methods

In this section, we propose the two methods of point estimation for the CR model using the GPHCS data from the GLD.

3.1.   MLE

In the presence of CR data, the GPHC schemes are used to acquire the ML estimators of the unknown parameters $\alpha_ 1$, $\alpha_2$, and $\beta$, the reliability characteristics $\Re(t)$ and $\hslash(t)$, and the relative risks $\pi_j, j = 1, 2$. Assuming that $\gamma$ is known, without any generality, $\gamma = \gamma_0$.

From Eqs (2.2) and (2.6), the likelihood function can be written as follows:

Therefore, taking the logarithms of likelihood function (3.1), say, $l^*$, then

Now, by differentiating $l^*$ with respect to $\alpha_{1}$, $\alpha_{2}$, and $\beta$, then

and

Equating (3.3) and (3.4) with zero, then

where

To solve Eq (3.6), we need to get the ML estimator of $\beta$ say $\widehat{\beta}_{ML}$. There is no closed-form solution for $\beta$, so by setting Eq (3.5) with zero and solving numerically, the $\widehat{\beta}_{ML}$ is produced. Now, by substituting (3.6) into (3.5), we get the following:

In this situation, one can utilize numerical methods such as the Newton-Raphson method to solve for $\beta$ and then use this value to determine $\widehat{\alpha}_{1ML}$ and $\widehat{\alpha}_{2ML}$, respectively. In the real data example, we used profile plots of the log-likelihood function of parameter $\beta$, as shown in Figure 3. These findings suggest that the required ML estimators may exist uniquely; see Section 7.

Using the invariance property, the ML estimator of the reliability characteristics $\widehat{\Re}(t)$ and $\widehat{\hslash}(t)$, at the various time $t$, can be acquired from (2.3) and (2.4) by replacing the original values $\alpha_{1}$, $\alpha_{2}$, and $\beta$ by their ML estimators $\widehat{\alpha}_{1ML}, \widehat{\alpha}_{2ML}$, and $\widehat{\beta}_{ML}$, respectively.

In the same way, the ML estimator of the relative risk resulting from Causes 1 and 2 may be determined from Eq (2.5) as follows:

3.2.   Bayesian estimation

In this subsection, we discuss Bayesian estimators for the unknown parameters $\alpha_1, \alpha_2$, and $\beta$ of the GLD based on the CR model under GPHCS. We assumed independent gamma with hyper-parameters $a_l > 0$ and $b_l > 0, \ l = 1, 2, 3$. The following forms can be acquired from the combined prior distribution of these parameters:

The joint posterior distribution of unknown parameters $\alpha_{1}$, $\alpha_{2}$, and $\beta$ is determined by utilizing the likelihood function (3.1) and the joint prior distributions (3.8) as follows:

where

In Bayesian analysis, a loss function must be specified to estimate an unknown parameter. There are no hard-and-fast guidelines for choosing the best loss function; it depends on the particular issue. The two main categories of loss functions are symmetric and asymmetric. The Bayesian estimator for the function $u(\alpha_1, \alpha_2, \beta)$ under symmetric (SE) and asymmetric (LINEX and GE) loss functions are driven. A common symmetric loss that penalizes overestimation and underestimating equally is the SE. Because of its convexity, which guarantees a unique minimum and makes optimization easier, it is frequently utilized in regression situations (Varde^[29]). The posterior mean is the Bayesian estimator for a function $u(\alpha_1, \alpha_2, \beta)$ under the SE loss function and is computed by

In machine learning, the LINEX is a commonly used statistic, especially for classification issues. It calculates how much the target variable's actual distribution differs from the predicated probability distribution (Varian^[30]). Although it provides flexibility in estimating location parameters, the LINEX loss function might not be the best option for estimating scale parameters (Basu and Ibrahimi^[31]. Based on the LINEX loss function, the Bayesian estimator for a function $u(\alpha_1, \alpha_2, \beta)$ is given by

A good alternative for the modified LINEX loss function is the GE, which Calabria and Pulcini ^[32] suggested. The GE loss function is a useful tool in information theory and machine learning. The Bayesian estimator for function $u(\alpha_1, \alpha_2, \beta)$ under the GE loss function is described below:

Moreover, from Eq (3.9), we can obtain the Bayesian estimators of the reliability characteristics at operation time $t$ and relative risks, respectively, as follows:

Since it is extremely difficult to solve the integrals from (3.10) to (3.15) analytically, we shall use the MCMC approach. The lower and upper limits of the BCI and HPD intervals for the unknown parameters, reliability characteristics, and relative risks will also be calculated using the MCMC approach's sample data. Therefore, the conditional posterior probability for $\alpha_1$ given $\alpha_2, \beta$, and $\underline{x}$, and $\alpha_2$ given $\alpha_1, \beta$, and $\underline{x}$ are, respectively, expressed as follows:

Also, the conditional posterior probability for $\beta$ given $\alpha_1, \alpha_2$, and $\underline{x}$ is given by:

Because the conditional probability distribution in Eq (3.18) is not well-known, we will apply the Metropolis-Hastings sampler to obtain the values of $\beta$ that follow the distribution in (3.18). Using the Metropolis-Hastings technique, we can produce random samples from the normal proposal distribution; see Metropolis et al.^[33]. The conditional posterior distributions $\alpha_1$, $\alpha_2$, and $\beta$ are generated using Algorithm A.2, as shown in Appendix A.

The Bayesian estimators based on the three loss functions and the generated values $\varphi_g^{(z)}$; $g = 1, 2, ..., 7$; and

respectively, are obtained as follows:

If the values of the hyper-parameters tend to be zero, then the Bayesian estimator for vector $\varphi_g$ can be computed using the same approach as before. In this case, the prior is a uniform prior.

4.   CIs

This section discusses two methods of constructing CIs for $\alpha_1$, $\alpha_2$, $\beta$, reliability characteristics and relative risks. The first method uses the asymptotic distributions of the estimators to obtain CIs of the different parameters of interest. The second method is the parametric bootstrap CIs, constructed using the percentile bootstrap (Boot-p) and the bootstrap-t (Boot-t) methods.

4.1.   ACIs

This subsection presents the ACIs of the unknown parameters, reliability characteristics and relative risks. The ACIs can be acquired by inverting the Fisher information matrix (FIM) with the negative elements of expected values of the second-order derivatives of logarithms of the likelihood functions. Cohen ^[34] replaced expected values with their ML estimators to obtain the approximation variance-covariance matrix. Let

be the vector of unknown parameters, where $\Theta_1 = \alpha_1, \Theta_2 = \alpha_2$, and $\Theta_3 = \beta$. Let $I^{-1}(\underline{\Theta})$ denote the inverse of FIM of the parameters,

From Eqs (3.3)–(3.5), the second derivatives of $l^*$ with respect to $\alpha_{1}, \alpha_{2}$ and $\beta$ are

where

Therefore, we obtain the ACIs of the unknown parameters $\alpha_{1}, \alpha_2$ and $\beta$ for $j_{1}, j_{2} > 0$. Depending upon the asymptotic distribution of the ML estimators of the parameters, it is considered as the following:

Now, the two-sided $100(1-\delta)\%$, $0 < \delta < 1$, ACIs for the vector of unknown parameters can be obtained as follows:

where $Z_{\delta/2}$ is the percentile of the standard normal distribution with right-tail probability $\delta/2$, and $Var(\widehat{\Theta}_{l})$ is the element of the main diagonal of $I^{-1}(\widehat{\Theta}_{l})$ for $l = 1, 2, 3$.

The obtained ACIs may occasionally have negative lower bounds. To overcome this issue, the logarithmic transformation and delta techniques can be used to construct the asymptotic normality distribution of $\ln \widehat{\Theta}_l$ as

As a result, the formula for a $100(1-\delta)\%$, $0 < \delta < 1$, ACI of $\Theta_l$ obtained in this way is

where

Furthermore, based on the ML estimators' asymptotic normality, it is recognized that $\widehat{\Re}(t) \thicksim N(\Re(t)$, $\widehat{\sigma}_{\Re}^2), \widehat{\hslash}(t) \thicksim N(\hslash(t), \widehat{\sigma}_{\hslash}^2)$, and $\widehat{\pi}_j \thicksim N(\pi_j, \widehat{\sigma}_{\pi_j}^2); j = 1, 2$. It is possible to create the ACIs for $\Re(t), \hslash(t)$, and $\pi_j$ by employing the corresponding normality. Therefore, Greene's ^[35] delta approach must be used to approximate and represent the variances of the estimators of $\Re(t), \hslash(t)$, and $\pi_j$. The three partial derivative vectors for $\Re(t), \hslash(t)$, and $\pi_j$ with respect to unknown parameters are known as $\triangle_{\Re}, \triangle_{\hslash}$, and $\triangle_{\pi_j}; j = 1, 2$, respectively, which are described as follows:

where

Consequently, it is possible to acquire the approximate variances of $\Re(t), \hslash(t)$, $\pi_j; j = 1, 2$ as follows:

Now, the two-sided $100(1-\delta)\%$, $0 < \delta < 1$, ACIs for $\Re(t), \hslash (t)$, and $\pi_j$ can be obtained as follows:

4.2.   Boot-p CIs

This subsection constructs two parametric Boot-p CIs for unknown parameters, reliability characteristics, and relative risks of the GLD in the presence of two CR based on GPHC schemes, which are known as the Boot-p (Efron^[36]) and the Boot-t (Hall^[37]) methods.

Let

be a vector of the unknown parameters, reliability characteristics, and relative risks, where

By following Kundu et al.^[38] and Kundu and Joarder^[11], we generate Boot-p CIs by using Algorithm 1, as shown in Appendix A.

5.   Simulation study

A Monte Carlo analysis is used in this section to compare the performance of the MLEs and the BEs for gamma distribution under different GPHC schemes. In addition, the BEs using the gamma prior distribution are compared. We generated GPHC data for each failure sample. In addition, we identified the cause of failure as 1 or 2 with a probability of $\dfrac{ \alpha_1}{ \alpha_1+ \alpha_2}$ and $\dfrac{ \alpha_2}{ \alpha_1+ \alpha_2}$, according to an algorithm proposed by Balakrishnan and Sandhu^[39]. With the following assumptions, one can generate 10,000 GPHC samples, 11,000 MCMC samples with 1000 nburn, and 1000 Boot-p samples in the presence of CR from the GLD.

$(1)$ Under the prefixed time $\tau = 1.2$, assume the sample sizes for each cause of failure, and the number of failed units with the various minimum number of failure units are

$(2)$ Assume that the following censoring schemes will be used to remove the remaining units:

– Scheme 1: $R_{1} = n-m$; otherwise, $R_{i} = 0$.

– Scheme 2: $R_{i} = 1$ if $i = 1, 2, ..., n-m$; otherwise, $R_{i} = 0$.

– Scheme 3: $R_{m} = n-m$ and $R_{i} = 0$ if $i$ otherwise.

$(3)$ In this study, we assumed fixed parameters $\gamma_0 = 3$ and $\gamma_0 = 1$ to use for choosing the GLD parameters

It is important to mention that in both states 1 and 2, $X_i = \min(X_{1i}, X_{2i})$ follows the Burr type XII($\beta = 1$) and Lomax($\gamma_0 = 1$) distributions, as special states of GLD, respectively.

$(4)$ From BE, the three loss functions, SE, LINEX $(h = 1.5$ and $-1.5)$, and GE $(c = 0.5$ and $-0.5)$, are used. The hyper-parameters of gamma prior to BE are obtained in these states:

$(a)$ State 1: $(a_l, b_l) = {(5.00, 10.00), (11.25, 15.00), (20.00, 20.00)}; l = 1, 2, 3$, when $\gamma_0 = 3$.

$(b)$ State 2: $(a_l, b_l) = {(12.8, 16), (45.0, 30.0), (80.0, 40.0)}; l = 1, 2, 3$, when $\gamma_0 = 1$.

Let the mean of the marginal prior distribution be $\Theta_{true}$ and its variance be 0.05 for each prior see El-Din et al.^[40].

We used MSE, absolute bias (AB), and average estimate (AE) to compare each estimate. The values of AB and MSE of MLEs and BEs for the parameters $\alpha_1$, $\alpha_2$ and $\beta$ depending on gamma priors are presented in Tables B.1–B.6. In addition, for the values of the AE of the MLE and BEs for the reliability characteristics $\Re(t)$ and $\hslash(t)$ at different times see Tables B.7–B.10). From these tables, we noted that

● The BEs under the LINEX loss function (BL) with $h = -1.5$ have the highest MSE values among all BEs for the various estimates.

● The BEs under BL with $h = 1.5$ for $\widehat{\alpha}_{1BL}$ and $\widehat{\alpha}_{2BL}$ exhibit the lowest MSE values, as shown in Figure 4a, b.

● The BEs under the SE loss function (BS) for $\widehat{\beta}_{BS}$ had the lowest MSE values when $\gamma_0 = 1$; see Figure 4c.

● When $\gamma_0 = 3$, the MSE values of BEs exhibit similar behavior for different schemes, but the MSE values of MLEs decrease.

● For a fixed $m$, the MSE values of various estimates of $\alpha_1$ and $\beta$ decrease for different schemes as $k$ increases (see Figure 5a, c), whereas the MSE values of $\alpha_2$ increase and vice versa.

● The MSE values of $\alpha_2$ and $\beta$ estimates increase and vice versa for $\alpha_1$, for fixed $k$ and change $m$, as shown in Figure 5.

● Tables B.3–B.5 show that the MSE values of all estimates decrease for different schemes when $m$ increases with constant $k$ in most cases.

● The MSE values for different estimates of $\Re(t)$ increase as time $t$ increases, and vice versa for $\hslash(t)$, as shown in Figure 6 under $\gamma_0 = 1$.

● When $\gamma_0 = 3$, the MSE values of $\Re(t)$ and $\hslash(t)$ increased with time $t$.

● The BEs of $\Re(t)$ under BL with $h = -1.5$ have the minimum MSE values at $t = 0.1$, whereas the MSE values of BE under BL (with $h = 1.5$) are the smallest at $t = 0.5$, as shown in Figure 6a.

● The BEs exhibit the lowest MSE values for $\hslash(t)$ at time $t = 0.1$ when the loss function is the BS; see Figure 6b.

● The MSE values are the highest under BL (with $h = -1.5$) at various times $t$.

● Scheme 3 has the lowest MSE values under the ML method, as shown in Figures 4a, b and 7 concerning both $\Re(t)$ and $\hslash(t)$.

● In most instances, the MSE values of $\Re(t)$ and $\hslash(t)$ are approximately equal to those of the BEs; see Figure 7.

● At the time $t = 0.5$ and $\gamma_0 = 3$, the BEs for $\Re(t)$ under BL with $h = -1.5$ had the smallest MSE, as shown in Figure 7a, b.

● For $\gamma_0 = 3$ and at time $t = 0.1$, the BEs for $\hslash(t)$ under the GE loss function (BG) with $c = 0.5$ had the smallest MSE; however, at time $t = 0.5$, the BEs for $\hslash(t)$ under the BL with $h = 1.5$ had the smallest MSE; see Figure 7c, d.

Furthermore, Tables B.11–B.14 illustrate the average width (AW) and coverage probability (CP) for the $95\%$ ACI, BCIs, and HPDs, respectively, using the MCMC-generated values of the parameters and reliability characteristics. According to the tables, we noted that

$(1)$ The ACIs of all parameters have the largest CP; therefore, the AWs of these parameters are also the highest.

$(2)$ Boot-t has the smallest CP for $\alpha_1$ and $\alpha_2$, whereas the opposite is true for $\beta$. Furthermore, Boot-p has a smaller AW than Boot-t for all the parameters.

$(3)$ The HPDs had smaller AWs than BCI, which had a higher CP than the HPDs.

$(4)$ When $k = 30$ and $m$ increase, the CP for $\Re(t)$ increases over time in most cases and vice versa for $k = 25$.

$(5)$ For fixed $k$, the CP for $\hslash(t)$ increases over time in most cases as $m$ increases.

$(6)$ At $\gamma_0 = 3$ and for fixed $k$ and $m$, the AW for $\Re(t)$ increases over time, similar to that for $\hslash(t)$.

$(7)$ For fixed $k$, $m$ and at $\gamma_0 = 1$, the AW for $\Re(t)$ increases over time in most cases and vice versa for $\hslash(t)$, except in case ACI.

6.   Real data applications

This section examines the adaptability and validity of GLD using real-time data analysis. It also demonstrated the practical use of the proposed CR model based on the GPHCS. We display the original dataset from the jute fiber breaking strength experiment conducted (see Xia et al. ^[41]). In this experiment, the breaking strength failure data of the jute fiber was affected by two different gauge lengths: 5 mm and 10 mm. The following data is displayed using the transformation $Y = X/200$:

Time of failure due to Cause 1: 0.6454, 0.83935, 0.841, 0.89125, 0.9271, 0.9384, 1.0943, 1.13265, 1.27145, 1.30485, 1.341, 1.35395, 1.5242, 1.53495, 1.57665, 1.804, 1.8385, 1.8501, 2.20935, 2.47755, 2.4814, 2.5824, 2.68725, 2.73055, 2.77305, 2.83155, 2.91985, 3.09285, 4.04615, 4.11515.

Time of failure due to Cause 2: 0.21965, 0.2508, 0.50575, 0.5447, 0.6153, 0.7069, 0.7574, 0.817, 0.88625, 0.9158, 1.06065, 1.2872, 1.3145, 1.45635, 1.5195, 1.61915, 1.7662, 1.8821, 1.91715, 2.11055, 2.533, 2.65275, 2.9524, 3.1883, 3.35745, 3.46865, 3.5037, 3.5233, 3.63615, 3.89085.

We conducted the Kolmogorov-Smirnov (KS) test to evaluate whether the data conformed to the GLD, and the $P$-values obtained from the KS test for Causes 1 and 2 were 0.560938 and 0.651419, respectively. The GLD provided a reasonable fit for the data, see Figure 8. Subsequently, we employ a censored dataset and various values that differ from the original data (both big and small).

6.1.   Big data

In this subsection, we create a PTIIC sample of $m = 40$ out of $n = 60$. The PTIIC sample is presented in Table 2. This table displays the failure lifetime $y_{i:m:n}$, cause of failure $\delta_i$, and censoring scheme $R$ used to remove the remaining units. Three alternative GPHCS were produced when the known parameter $\gamma_0 = 3$ was selected, assuming a value of $k = 25$ and various values of $\tau$.

$1)$ Case 1: Suppose $\tau = 1.88259$, since $\tau < Y_{25:40:60}$, then the experiment would have terminated at $y_{25:40:60}$, with $J_1 = 10$, $J_2 = 15$, $R^* = (2, 0, 0, 2, 0, 0, 2, 0, 0, 1, 0, 0, 0, 1, 0, 2, 0, 1, 0, 1, 0, 0, 0, 0, 23$), $R_{k}^* = 23$; and we would have the following data: 0.50575, 0.5447, 0.7069, 0.7574, 0.817, 0.83935, 0.841, 0.88625, 0.89125, 0.9271, 0.9384, 1.06065, 1.0943, 1.13265, 1.27145, 1.341, 1.35395, 1.45635, 1.5242, 1.53495, 1.57665, 1.61915, 1.7662, 1.8821, 1.91715.

$2)$ Case 2: Suppose $\tau = 2.7$, since $Y_{25:40:60} < \tau < Y_{40:40:60}$, then the experiment would have terminated at $\tau = 2.7$, with $J_1 = 12, J_2 = 18, R^{**} = (, 0, 0, 2, 0, 0, 2, 0, 0, 1, 0, 0, 0, 1, 0, 2, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1$), $R_{\tau}^{**} = 15$; and we would have the following data: 0.50575, 0.5447, 0.7069, 0.7574, 0.817, 0.83935, 0.841, 0.88625, 0.89125, 0.9271, 0.9384, 1.06065, 1.0943, 1.13265, 1.27145, 1.341, 1.35395, 1.45635, 1.5242, 1.53495, 1.57665, 1.61915, 1.7662, 1.8821, 1.91715, 2.20935, 2.47755, 2.533, 2.5824, 2.65275.

$3)$ Case 3: Suppose $\tau = 4.5$, since $Y_{40:40:60} < \tau$, then the experiment would have terminated at $Y_{40:40:60} = 4.11515$, with $J_1 = 15$, $J_2 = 25$, $R^{***} = R$, $R_{m}^{***} = 0$; and we would have the following data: 0.50575, 0.5447, 0.7069, 0.7574, 0.817, 0.83935, 0.841, 0.88625, 0.89125, 0.9271, 0.9384, 1.06065, 1.0943, 1.13265, 1.27145, 1.341, 1.35395, 1.45635, 1.5242, 1.53495, 1.57665, 1.61915, 1.7662, 1.8821, 1.91715, 2.20935, 2.47755, 2.533, 2.5824, 2.65275, 2.77305, 2.83155, 2.9524, 3.09285, 3.1883, 3.5037, 3.5233, 3.63615, 4.04615, 4.11515.

Next, the estimate for parameter $\beta$ can be obtained numerically by setting $g(\beta) = 0$ and graphically by plotting the graph of $g(\beta)$ in Eq (3.7) and finding the point of intersection with the $\beta$-axis. Figure 9 illustrates the behavior of $g(\beta)$ for cases 1–3. Additionally, Figure 3 confirms that the MLE of $\beta$ is unique since it is the only one that maximizes the log-likelihood for each scheme. As a result, the MLEs for the parameters $\alpha_1, \alpha_2,$ and $\beta$ are unique and exist. Consequently, we can calculate the reliability characteristics $\Re(t)$ and $\hslash(t)$ at different times and relative risks due to Causes 1 and 2. Therefore, we utilized ML and Bayesian approaches to estimate parameter values. Since we have no prior knowledge of the unknown parameters, we compute BEs with a non-informative (uniform) prior, as seen in Kundu and Pradhan ^[42].

6.2.   Small data

This subsection details the creation of a PTIIC sample consisting of $m = 25$ items from a total of $n = 60$. Table 3 presents the PTIIC sample, showing the failure lifetime $y_{i:m:n}$, failure cause $\delta_i$, and the censoring scheme $R$ employed to eliminate the remaining units. Using a known parameter value of $\gamma_0 = 3$ and setting $k = 15$, three different GPHCS were generated with varying values of $\tau$.

$1)$ Case 1*: Suppose $\tau = 1.88432$, since $\tau < Y_{15:25:60}$, then the experiment would have terminated at $y_{15:25:60}$, with $J_1 = 6$, $J_2 = 9$, $R^* = (, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 24$), $R_{k}^* = 24$; and we would have the following data: 0.21965, 0.50575, 0.6153, 0.6454, 0.817, 0.83935, 0.841, 0.9271, 1.13265, 1.30485, 1.5242, 1.57665, 1.8385, 1.8501, 1.91715.

$2)$ Case 2*: Suppose $\tau = 2.5$, since $Y_{15:25:60} < \tau < Y_{25:25:60}$, then the experiment would have terminated at $\tau = 2.7$, with $J_1 = 8, J_2 = 9, R^{**} = (, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3, 0, 3$), $R_{\tau}^{**} = 16$; and we would have the following data: 0.21965, 0.50575, 0.6153, 0.6454, 0.817, 0.83935, 0.841, 0.9271, 1.13265, 1.30485, 1.5242, 1.57665, 1.8385, 1.8501, 1.91715, 2.20935, 2.47755.

$3)$ Case 3*: Suppose $\tau = 4.5$, since $Y_{25:25:60} < \tau$, then the experiment would have terminated at $Y_{25:25:60} = 4.11515$, with $J_1 = 11$, $J_2 = 14, R^{***} = R, R_{m}^{***} = 0$; and we would have the following data: 0.21965, 0.50575, 0.6153, 0.6454, 0.817, 0.83935, 0.841, 0.9271, 1.13265, 1.30485, 1.5242, 1.57665, 1.8385, 1.8501, 1.91715, 2.20935, 2.47755, 2.68725, 2.77305, 2.83155, 2.9524, 3.09285, 3.5233, 3.89085, 4.11515.

The parameter $\beta$ can be estimated numerically by solving $g(\beta) = 0$ and visually by identifying where the graph of $g(\beta)$ in Eq (3.7) crosses the $\beta$-axis. The behavior of $g(\beta)$ for cases 1*–3* is depicted in Figure 10. Furthermore, Figure 11 demonstrates that the MLE of $\beta$ is unique, as it is the sole value that optimizes the log-likelihood for each scenario.

Tables B.15 and B.18 display the standard error (SER), which is the root mean square error for the MLEs and BEs of the uniform prior distribution. Tables B.15 and B.18 show the result of real data when $(n, m, k) = (60, 40, 25)$ and $(n, m, k) = (60, 25, 15)$, respectively. From these tables, the SER for MLE was the smallest concerning most estimates. Moreover, the reliability characteristic $\Re(t)$ exhibits a decline from $t = 0.5$ to $t = 1.5$, whereas the opposite trend is observed for $\hslash(t)$. Concurrently, the SER for both reliability characteristics exhibited an upward trend over time. Across all scenarios, the relative risk associated with Cause 1 was consistently lower than that of Cause 2.

Tables B.16 and B.17 show the lower and upper limits for the AWs of each parameter scheme for ACI, Boot-t, Boot-p, and BCI. It is evident from this table that the AWs of the HPD for all schemes were smaller than those of the BCI. The AWs for estimates of $\alpha_1$, $\alpha_2$, and $\beta$ of Boot-t were found to be smaller than those of Boot-p, whereas the reliability characteristic $\Re(t)$ and $\hslash(t)$ at different times $t$ and relative risks due to Causes 1 and 2 showed the opposite trend. As time $t$ increased, the AWs for $\Re(t)$ and $\hslash(t)$ also increased. It is worth noting that BCI generally had the largest AWs in Table B.16, whereas in Table B.17, the AWs of ACI were the largest in most cases.

7.   Concluding remarks and future research

This study aimed to investigate the application of CR models for GLD in the context of GPHC. Within the CR model that incorporated the two failure mechanisms, it was postulated that the GLD independently determined the lifetime.

The GLD introduces a new parameter to the traditional Lomax distribution, making it more flexible for applications in diverse fields such as agriculture, finance, and actuarial sciences. This increased flexibility is particularly beneficial in lifetime studies and reliability analysis. For instance, the GLD can be effectively applied to analyze datasets like the failure times of Kevlar 373/epoxy under constant pressure. By incorporating CR and the GPHCS, we can extend this analysis to more complex scenarios.

In this paper, the point and interval estimators for the unknown parameters, reliability characteristics, and relative risks due to two causes were derived using both ML and Bayesian methods. Moreover, two Boot-p CIs based on MLEs were obtained. The BEs and their associated BCI and HPD intervals were calculated using the MCMC technique, which utilized three loss functions and gamma priors. A simulation analysis was conducted to assess the accuracy of the estimators. The results show the BEs of $\alpha_1$ and $\alpha_2$ under $(h = 1.5)$, whereas the BEs of $\beta$ based on the BS loss function have the lowest MSE values when $\gamma = 1$. BEs for the gamma prior outperformed the MLEs based on MSE values. The simulation study also revealed that as time $t$ increased and $\gamma_0 = 1$, the MSE values of the $\hslash(t)$ estimates decreased, while the MSE values of the $\Re(t)$ estimates increased. When $\gamma_0 = 3$, the MSE values of $\Re(t)$ and $\hslash(t)$ increased with time $t$ as shown in real data which the performance of the proposed estimators was evaluated using a real dataset, suggesting that the required MLEs may exist uniquely.

Future research directions include investigating lifetime distribution in the presence of CR and can be considered in the following scenarios as follows:

$1)$ Exploring scenarios with dependent competing risks.

$2)$ Employing alternative estimation techniques such as E-Bayesian method.

$3)$ Examining various censoring schemes, including those with random removals.

Author contributions

Amal Hassan: conceptualization, validation and review; Sudhansu Maiti: visualization and validation; Rana Mousa: writing, editing and methodology.; Najwan Alsadat: conceptualization, methodology, funding acquisition and review; Mahmoued Abu-Moussa: methodology, software, and editing. All authors have read and agreed to the published version of the manuscript.

Acknowledgments

This research is supported by researchers supporting project number (RSPD2024R548), King Saud University, Riyadh, Saudi Arabia.

Conflicts of interest

The authors declare no conflicts of interest.

Appendix

Appendix A

First, Appendix A consists of two algorithms used.

Next, all tables of numerical work results and real data outcomes are included in Appendix B.

Appendix B