Tampered Random Variable Analysis in Step-Stress Testing: Modeling, Inference, and Applications

Abstract

This study explores a new dimension of accelerated life testing by analyzing competing risk data through Tampered Random Variable (TRV) modeling, a method that has not been extensively studied. This method is applied to simple step-stress life testing (SSLT), and it considers multiple causes of failure. The lifetime of test units under changeable stress levels is modeled using Power Rayleigh distribution with distinct scale parameters and a constant shape parameter. The research introduces unique tampering coefficients for different failure causes in step-stress data modeling through TRV. Using SSLT data, we calculate maximum likelihood estimates for the parameters of our model along with the tampering coefficients and establish three types of confidence intervals under the Type-II censoring scheme. Additionally, we delve into Bayesian inference for these parameters, supported by suitable prior distributions. Our method’s validity is demonstrated through extensive simulations and real data application in the medical and electrical engineering fields. We also propose an optimal stress change time criterion and conduct a thorough sensitivity analysis.

1. Introduction

With ongoing improvements in the manufacturing sector, numerous industrial products, known for their high reliability and complex designs, are becoming increasingly usable in everyday life. Accelerated life testing (ALT) addresses the challenge of evaluating such products by exposing them to stress levels higher than their usual operating conditions, producing rapid failures in turn. These growing stress factors—such as temperature, voltage, and humidity—significantly influence the lifespan of electronic equipment, including electric bulbs, fans, computers, toasters, and more. By employing these high-stress factors in ALT experiments, valuable insights concerning product reliability can be rapidly acquired within a condensed experimental time frame. Analyzing reliability and making inferences from it have gained significant interest in the literature, as illustrated by references [1,2,3].

ALT experiments can be conducted in two ways: with a starting constant high-stress level or with a changeable stress factor that can be varied during different time intervals. In the realm of ALT, there exists a specific class known as step-stress life testing (SSLT). This method permits experimenters to incrementally increase the stress levels at predetermined time points during the experiment. A basic form of SSLT is exemplified in experiments involving only two stress levels, denoted as $s_1$ and $s_2$, along with a single, pre-determined point in time, $\tau$, at which the stress level shifts.

To understand how lifetime distributions vary under different stress levels, some basic modeling assumptions are typically discussed:

• Cumulative Exposure models (CE). In this method, specific restrictions are applied to ensure that the lifetime distributions at each progressive stress level align at their designated transition points, maintaining continuity. This approach is detailed in works by Sedyakin [4] and Nelson [5].
• Tampered Failure Rate (TFR) modeling. This technique involves adjusting the failure rates, increasing them at each subsequent stress level. Key references for TFR modeling include Bhattacharyya and Soejoeti [6], and Madi [7].
• The Tampered Random Variable (TRV) model. Here, the focus is on reducing the remaining lifetime for each new stress level. For more information on this approach, you can refer to Goel [8] and DeGroot and Goel [9].
• Step-stress partially accelerated life testing with a large amount of censored data. This approach addresses the gap in estimating non-homogeneous distribution and acceleration factor parameters under multiple censored data conditions. For more details, one can refer to Khan and Aslam [10].

Additionally, Sultana and Dewanji [11] explored the relationships between the TRV model and the two other models, TFR and CE, within a multi-step stress environment. They noted that TRV modeling aligns with CE and TFR when the fundamental lifetime distribution is exponential and the distributions at each stress level adhere to a scale-based parametric family. Thus, it is observed that the above three models converge when the fundamental distribution is exponential. TRV modeling stands out for its ability to be generalized to multiple-step-stress situations more effectively than the other two models. It also offers advantages in terms of modeling discrete and multivariate lifetimes, which are more complex tasks for the CE and TFR models.

Comparing factors that lead to risk model failures is essential for comprehending the contributing factors, detecting common changes, assessing model performance, and influencing decision making and risk management. It assists in identifying important issues that must be resolved to increase the precision and dependability of the model. The development of more reliable models can be made possible by recognizing similar patterns among various occurrences or outcomes that can be identified by understanding these components. Additionally, it offers useful information for model creators and validators, enabling them to improve model development processes, assumptions, and validation processes for more accurate and dependable models. The competing risks concept refers to the possibility of individual failure in a specific field owing to distinct factors. The cause-of-failure indication and the individual failure time are examples of observable data in this approach. When examining data on competing risks, the failure variables are typically unrelated to one another which means that the two risk factors are statistically independent. In the industrial and mechanical domains, fatigue and aging deterioration can lead to an assembly device failing due to electrical/optical signal (voltage, current, or light intensity) falling to an intolerable level. Numerous studies in the existing literature utilize CEM and TFR modeling within competing risk scenarios. However, to the best of our knowledge, research incorporating TRV modeling into the context of competing risk data is notably scarce. See for example Sultana et al. [12], Ramadan et al. [13] and Tolba et al. [14].

In this work, TRV is used with the SSLT model under two independent competing risk factors where the failure times follow the Power Rayleigh distribution. The sample is observed under the Type-II censored scheme. The censoring schemes have been introduced to solve the lack of information in lifetime experiments, saving time and costs. Type-I censoring has a predetermined time, while Type-II censoring has predetermined failure units.

The main goals of this study are summarized below:

• Performing an inferential analysis to obtain point and interval estimation of the unknown parameters of the distribution and the acceleration factor using both the maximum likelihood estimator and the Bayesian method.
• Applying numerical methods like Monte Carlo simulation to assess the performance of estimators obtained from Maximum Likelihood Estimation (MLE) and Bayesian methods, focusing on their bias, mean squared error, and the coverage probability (CP) for the confidence intervals.
• Evaluating real-world data sets from the medical field concerning AIDS infection, alongside another study from electrical engineering involving the causes of the failure of electronic components, serves to empirically assess the effectiveness of the newly proposed model.

The structure of the remainder of this document is as follows: Section 2 outlines the SSLT model under the TRV framework with the Power Rayleigh distribution. Section 3 details the methodologies used for point estimation, specifically using maximum likelihood and Bayesian methods. Section 4 is dedicated to interval estimation, exploring three distinct methods. Section 5 focuses on simulation analysis and presents the results in tabular form. The determination of the optimal time for stress change and an analysis of sensitivity are discussed in Section 6. An application using real-world data is examined in Section 7. The paper concludes with a summary of the findings in Section 8.

2. Model Description

In this study, we consider the SSLT model with random failure time variables denoted by $U_1$ and $U_2$ along with the stress levels $s_1$ and $s_2$ that are assumed to follow a Power Rayleigh distribution with a common shape parameter $\gamma$ and distinct scale parameters ${\lambda }_1$ and ${\lambda }_2$. The two risk factors are called cause I and cause II and both are performed using Type-II censored samples. At a prefixed time $\tau$, the stress level moves from $s_1$ to $s_2$. During the first stress level, the $s_1$ units will operate until a specific time $\tau$, following which any remaining survivals that have not failed by time $\tau$ are moved to be tested under accelerated conditions with an acceleration factor $\beta$. Consequently, the system will operate under the second stress level $s_2$ until we obtain the required failure times. The effect of stress transition from the first stress to the accelerated condition may be explained by multiplying the remaining lifetime by the acceleration factor $\beta$. Hence the TRV for $U_1$ and $U_2$ is expressed as

$$\widetilde{U_1}=\left\{\begin{array}{ll}{U}_1& ,0<t\le \tau \\ \tau +\beta (U_1-\tau )& ,t>\tau \end{array}\right.$$

(1)

and

$$\widetilde{U_2}=\left\{\begin{array}{ll}{U}_2& ,0<t\le \tau \\ \tau +\beta (U_2-\tau )& ,t>\tau ,\end{array}\right.$$

(2)

where $\tau$ is the time at which the stress changes and the acceleration parameter is $0<\beta <1$.

We consider Power Rayleigh distribution as a lifetime model. The Rayleigh distribution, a continuous distribution of significant practical relevance, has been the subject of extensive study by various authors who have explored its statistical properties, inference methods, and reliability analysis. Additionally, a variety of extended versions of the Rayleigh distribution have been introduced. For example, Rosaiah and Kantam [15] applied the inverse Rayleigh to failure times data. Merovci [16] introduced the transmuted Rayleigh and modeled the amount of nicotine in blood. Cordeiro et al. [17] studied the beta-generalized Rayleigh distribution and its application. More generalizations of Rayleigh distribution can be found in the literature and one may refer to [18,19,20,21,22,23,24].

The Power Rayleigh (PR) distribution was first introduced by Neveen et al. [25]. It is a versatile and flexible statistical model known for its ability to handle a wide range of data types. This distribution is particularly useful due to its capability to model data that exhibit a skewed pattern, which is common in many practical situations. The Power Rayleigh distribution is characterized by two parameters that allow it to adapt to various data shapes and sizes, making it more flexible than the standard Rayleigh distribution. Its applications are diverse, ranging from reliability engineering and survival analysis to modeling wind speed and signal processing. The flexibility in shape and scale provided by the Power Rayleigh distribution makes it a valuable tool for analyzing and interpreting real-world data in various scientific and engineering fields. We assume that the Power Rayleigh distribution has a shape parameter $\gamma$ and scale parameter $\lambda$, where both have positive support, and then the cumulative distribution function (CDF) becomes

$$F\left(t\right)=1-e^{-\frac{t^{2\gamma }}{2{\lambda }^2}},$$

and the probability density function (PDF) is

$$f\left(t\right)=\frac{\gamma }{{\lambda }^2}{t}^{2\gamma -1}{e}^{-\frac{t^{2\gamma }}{2{\lambda }^2}}.$$

Consider a set of n units subjected to a life test starting at stress level $s_1$. Failures and their corresponding risks are documented over time. At a designated moment $\tau$ the stress level shifts from $s_1$ to $s_2$, and the test runs until r (with $r<n$) failures are noted. If r equals n, a complete dataset is collected as in a simple SSLT without data truncation. We assume that each unit’s failure is attributable to one of two competing risks, each described by a Power Rayleigh distribution with a consistent shape parameter $\gamma$ but distinct scale parameters ${\lambda }_j$ for $j=1,2$, aligned with the TRV model.

The CDF for the lifetime $U_j$ associated with risk j for $j=1,2$ is then expressed as follows:

$$F_j\left(t\right)=F_j(t;\gamma ,{\lambda }_j)=\left\{\begin{array}{cc}1-e^{-\frac{t^{2\gamma }}{2{\lambda }_j^2}}\hfill & \mathrm{if}0<t\le \tau \hfill \\ 1-e^{-\frac{{(\tau +{\beta }^{-1}(t-\tau ))}^{2\gamma }}{2{\lambda }_j^2}}\hfill & \mathrm{if}t>\tau \hfill \end{array}\right.$$

(3)

and the corresponding PDF of $U_j$ is given by

$$f_j\left(t\right)=f_j(t;\gamma ,{\lambda }_j)=\left\{\begin{array}{cc}\frac{\gamma }{{\lambda }_j^2}{t}^{2\gamma -1}{e}^{-\frac{t^{2\gamma }}{2{\lambda }_j^2}}\hfill & \mathrm{if}0<t\le \tau \hfill \\ \frac{\gamma }{{\lambda }_j^2\beta }{(\tau +{\beta }^{-1}(t-\tau ))}^{2\gamma -1}{e}^{-\frac{{(\tau +{\beta }^{-1}(t-\tau ))}^{2\gamma }}{2{\lambda }_j^2}}\hfill & \mathrm{if}t>\tau \hfill \end{array}\right.$$

(4)

Let us denote the overall failure time of a unit under test as U, which is obtained by $U=min\{U_1,U_2\}$. Then, the CDF and PDF are easily obtained as

$$F\left(t\right)=F(t;\gamma ,\lambda )=1-\left(1-F_1\left(t\right)\right)\left(1-F_2\left(t\right)\right)=\left\{\begin{array}{cc}1-e^{-\Lambda t^{2\gamma }}\hfill & \mathrm{if}0<t\le \tau \hfill \\ 1-e^{-\Lambda {(\tau +{\beta }^{-1}(t-\tau ))}^{2\gamma }}\hfill & \mathrm{if}t>\tau \hfill \end{array}\right.$$

(5)

and

$$f\left(t\right)=f(t;\gamma ,\lambda )=\left\{\begin{array}{cc}\gamma \Lambda t^{2\gamma -1}{e}^{-\Lambda t^{2\gamma }}\hfill & \mathrm{if}0<t\le \tau \hfill \\ \frac{\gamma }{\beta }\Lambda {(\tau +{\beta }^{-1}(t-\tau ))}^{2\gamma -1}{e}^{-\Lambda {(\tau +{\beta }^{-1}(t-\tau ))}^{2\gamma }}\hfill & \mathrm{if}t>\tau ,\hfill \end{array}\right.$$

(6)

respectively, where $\lambda =({\lambda }_1,{\lambda }_2)$ and $\Lambda =\frac{1}{2{\lambda }_2^2}+\frac{1}{2{\lambda }_1^2}$. Furthermore, let C denote the indicator for the cause of failure. Then, under our assumptions, the joint PDF of $(U,C)$ is given by

$$f_{U,C}(t;\gamma ,\lambda )=f_j\left(t\right)\left[1-F_k\left(t\right)\right]=\left\{\begin{array}{cc}\frac{\gamma }{{\lambda }_j^2}{t}^{2\gamma -1}{e}^{-\Lambda t^{2\gamma }}\hfill & \mathrm{if}0<t\le \tau \hfill \\ \frac{\gamma }{{\lambda }_j^2\beta }{(\tau +{\beta }^{-1}(t-\tau ))}^{2\gamma -1}{e}^{-\Lambda {(\tau +{\beta }^{-1}(t-\tau ))}^{2\gamma }}\hfill & \mathrm{if},t>\tau ,\hfill \end{array}\right.$$

(7)

for $j,k=1,2,j\ne k$.

In competing risk models, the assumption of independence is often considered to be impractical. Identifiable issues may emerge if dependencies exist within the model or due to a lack of covariates in the data. To mitigate these issues, we postulate a latent failure time model and treat the risks $U_1$ and $U_2$ as independent. Let $N_{j1}$ represent the number of units failing from risk j before time $\tau$ and $N_{j2}$ after $\tau$, with $N_j=N_{j1}+N_{j2}$, ensuring $N_1+N_2\le r$. The sequence of observed failure times is $0<t_{1:n_1}<t_{2:n_2}<\cdots <t_{r:n}$. Let ${\widehat{n}}_1$ denote the observed value for $N_1$, ${\widehat{n}}_2$ denote the observed value for $N_2$, and let $N=(N_1,N_2)$ be the vector of these counts.

In the next section, classical and Bayesian estimation methods are constructed to estimate the unknown parameters for the Power Rayleigh and the accelerated constant $\beta$ under the two competing risk factors with the Type-II censoring scheme.

3. Point Estimation

In this study, two approaches to estimation are examined: the frequentist maximum likelihood estimation (MLE) and the Bayesian estimation method. Section 4 is dedicated to conducting a simulation analysis and applying numerical techniques to evaluate the efficacy of these estimation strategies.

3.1. Maximum Likelihood Estimation

In this section, the maximum likelihood estimation MLE method is employed to determine the unknown parameters of the Power Rayleigh distribution within the TRV framework. Numerical methods, including the renowned Newton–Raphson technique, are utilized to compute the necessary estimators. Subsequently, assuming the TRV model, we construct the likelihood function of $\psi =(\gamma ,{\lambda }_1,{\lambda }_2,\beta )$ based on Type-II censored data as

$$L\left(\psi \right|(t,c))=\frac{n!}{(n-r)!}\prod _{i=1}^{\widehat{n_1}}{f}_{U,C}(t_{i:n},c_i)\prod _{i=\widehat{n_1}+1}^r{f}_{U,C}(t_{i:n},c_i){\left[1-F\left(t_{r:n}\right)\right]}^{n-r}.$$

Here $r={\widehat{n}}_1+{\widehat{n}}_2=n_{11}+n_{12}+n_{21}+n_{22}$. By substituting Equations (5) and (7) into the above likelihood equation we obtain

$$L\left(\psi \right|(t,c))=\frac{n!}{(n-r)!}\left(\frac{{\gamma }^r{\beta }^{-\widehat{n_2}}}{{\lambda }_1^{2\widehat{n_1}}{\lambda }_2^{2\widehat{n_2}}}\right)\prod _{i=1}^{\widehat{n_1}}{t}_i^{2\gamma -1}{e}^{-\Lambda t_i^{2\gamma }}\prod _{i=\widehat{n_1}+1}^r{(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma -1}{e}^{-\Lambda {(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma }}{e}^{-\Lambda (n-r)t_{r:n}^{2\gamma }}.$$

(8)

The log-likelihood function can be written as

$$\begin{array}{cc}\hfill \mathcal{l}\left(\psi \right)=& r\log \gamma -\widehat{n_2}\log \beta -2\widehat{n_1}\log {\lambda }_1-2\widehat{n_2}\log {\lambda }_2+(2\gamma -1)\left[\sum _{i=1}^{\widehat{n_1}}\log \left(t_i\right)+\sum _{i=n_1+1}^r\log (\tau +{\beta }^{-1}(t_i-\tau ))\right]-\hfill \\ \hfill & \Lambda \left[\sum _{i=1}^{\widehat{n_1}}{t}_i^{2\gamma }+\sum _{i=\widehat{n_1}+1}^r{(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma }\right]+(n-r)t_{r:n}^{2\gamma }.\hfill \end{array}$$

(9)

The maximum likelihood estimations of the parameters $(\gamma ,{\lambda }_1,{\lambda }_2,\beta )$ are obtained by differentiating the log-likelihood function $\mathcal{l}\left(\psi \right)$ with respect to the parameters $(\gamma ,{\lambda }_1,{\lambda }_2,\beta )$ and setting the result to zero, so we have the following normal equations.

$$\begin{array}{cc}\hfill \frac{\partial \mathcal{l}\left(\psi \right)}{\partial \gamma }=& \frac{r}{\gamma }+2\left[\sum _{i=1}^{\widehat{n_1}}\log \left(t_i\right)+\sum _{i=\widehat{n_1}+1}^r\log (\tau +{\beta }^{-1}(t_i-\tau ))\right]-\hfill \\ & \Lambda \left[\sum _{i=1}^{\widehat{n_1}}{t}_i^{2\gamma }\log \left(t_i\right)+\sum _{i=\widehat{n_1}+1}^r{(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma }\log (\tau +{\beta }^{-1}(t_i-\tau ))\right]+(n-r)t_{r:n}^{2\gamma }\log \left(t_{r:n}\right),\hfill \end{array}$$

(10)

$$\begin{array}{c}\hfill \frac{\partial \mathcal{l}\left(\psi \right)}{\partial \beta }=-\frac{\widehat{n_2}}{\beta }-(2\gamma -1)\sum _{i=\widehat{n_1}+1}^r\frac{{\beta }^{-2}(t_i-\tau ))}{(\tau +{\beta }^{-1}(t_i-\tau ))}+\Lambda \sum _{i=\widehat{n_1}+1}^r\gamma {\beta }^{-2}(t_i-\tau ){(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma -1},\end{array}$$

(11)

$$\begin{array}{c}\hfill \frac{\partial \mathcal{l}\left(\psi \right)}{\partial {\lambda }_1}=\frac{-2\widehat{n_1}}{{\lambda }_1}+{\lambda }_1^{-3}\left[\sum _{i=1}^{\widehat{n_1}}{t}_i^{2\gamma }+\sum _{i=\widehat{n_1}+1}^r{(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma }+(n-r)t_{r:n}^{2\gamma }\right]\end{array}$$

(12)

and

$$\begin{array}{c}\hfill \frac{\partial \mathcal{l}\left(\psi \right)}{\partial {\lambda }_2}=\frac{-2\widehat{n_2}}{{\lambda }_2}+{\lambda }_2^{-3}\left[\sum _{i=1}^{\widehat{n_1}}{t}_i^{2\gamma }+\sum _{i=\widehat{n_1}+1}^r{(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma }+(n-r)t_{r:n}^{2\gamma }\right].\end{array}$$

(13)

For known $\gamma$ and $\beta$, the MLEs of ${\lambda }_1$ and ${\lambda }_2$ are given by

$${\widehat{\lambda }}_1=\frac{1}{\sqrt{2}{\left(\widehat{n_1}\right)}^{\frac{1}{2}}}{\left[\sum _{i=1}^{\widehat{n_1}}{t}_i^{2\gamma }+\sum _{i=\widehat{n_1}+1}^r{(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma }+(n-r)t_{r:n}^{2\gamma }\right]}^{\frac{1}{2}}$$

and

$${\widehat{\lambda }}_2=\frac{1}{\sqrt{2}{\left(\widehat{n_2}\right)}^{\frac{1}{2}}}{\left[\sum _{i=1}^{\widehat{n_1}}{t}_i^{2\gamma }+\sum _{i=\widehat{n_1}+1}^r{(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma }+(n-r)t_{r:n}^{2\gamma }\right]}^{\frac{1}{2}}.$$

To address the system of nonlinear equations presented in Equations (10)–(13), numerical approaches are essential. Various numerical methods have been applied in existing research; in this instance, we employ the Newton–Raphson method. The outcomes of this application are detailed in Section 5.

3.2. Bayesian Inference

In this section, we apply the Bayesian estimation technique to determine the unknown parameters of the Power Rayleigh distribution. The fundamental principle of the Bayesian approach posits that the model’s parameters are random variables with a predefined distribution, referred to as the prior distribution. Given the availability of prior knowledge, selecting an appropriate prior is crucial. We opt for the gamma conjugate prior distribution for the parameters for many reasons, such as the flexibility in its nature with a non-informative domain and the calculations’ simplicity making analytical or computational updates to the posterior easier. Also, the positive of the domain makes it suitable for modeling parameters. We perform the Bayesian inference method for estimating the unknown parameters $\psi$ $=(\gamma ,{\lambda }_1,{\lambda }_2,\beta )$. We assume independent gamma priors for $\gamma , {\lambda }_1$, and ${\lambda }_2$ and a uniform prior for $\beta$. That is, $\gamma , {\lambda }_1$ and ${\lambda }_2$ have $Gamma(c_i,d_i),$ where $c_i,d_i$$>0$,$i=1,2,3$, are non-negative hyperparameters, and $\beta$ follows uniform prior as follows:

$$\pi \left(\beta \right)=1,0<\beta <1.$$

The estimates have been developed under the square error loss function (SELF) and the linear exponential loss function (LLF). Hence, the joint prior density of the independent parameters is given by

$$\pi \left(\psi \right)=\pi \left(\gamma \right)\pi \left({\lambda }_1\right)\pi \left({\lambda }_2\right)\pi \left(\beta \right),$$

$$\pi \left(\psi \right)={\gamma }^{c_1-1}{e}^{-d_1\gamma }{\lambda }_1^{c_2-1}{e}^{-d_2{\lambda }_1}{\lambda }_2^{c_3-1}{e}^{-d_3{\lambda }_2},\gamma >0,{\lambda }_1>0,{\lambda }_2>0,0<\beta <1.$$

(14)

The joint posterior density function for the parameters can be derived by incorporating the observed censored samples, and the prior distributions of these parameters are as follows:

$$\begin{array}{cc}\hfill {\pi }^{\ast }(\psi ,t,c)=& \pi \left(\psi \right)L\left(\psi \right|t,c)\hfill \\ \hfill =& {\gamma }^{c_1-1}{e}^{-d_1\gamma }{\lambda }_1^{c_2-1}{e}^{-d_2{\lambda }_1}{\lambda }_2^{c_3-1}{e}^{-d_3{\lambda }_2}\frac{{\gamma }^r{\beta }^{-\widehat{n_2}}}{{\lambda }_1^{2\widehat{n_1}}}{\lambda }_2^{2\widehat{n_2}}\prod _{i=1}^{{\widehat{n}}_1}{t}_i^{2\gamma -1}{e}^{-\Lambda t_i^{2\gamma }}\times \hfill \\ & \prod _{i={\widehat{n}}_1+1}^r{(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma -1}{e}^{-\Lambda {(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma }}{e}^{-\Lambda (n-r)t_{r:n}^{2\gamma }}.\hfill \end{array}$$

(15)

Thus, the conditional posterior densities of the parameters $\gamma ,{\lambda }_1,{\lambda }_2$, and $\beta$ can be obtained by simplifying Equation (15) as follows

$$\begin{array}{cc}\hfill {\pi }_1^{\ast }\left({\lambda }_1\right|{\lambda }_2,\gamma ,\beta ,t,c)=& {\lambda }_1^{-2\widehat{n_1}+c_2-1}{e}^{-d_2{\lambda }_1}\prod _{i=1}^{{\widehat{n}}_1}{e}^{-\Lambda t_i^{2\gamma }}\prod _{i={\widehat{n}}_1+1}^r{e}^{-\Lambda {(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma }}{e}^{-\Lambda (n-r)t_{r:n}^{2\gamma }},\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill {\pi }_2^{\ast }\left({\lambda }_2\right|{\lambda }_1,\gamma ,\beta ,t,c)=& {\lambda }_2^{2\widehat{n_2}+c_3-1}{e}^{-d_3{\lambda }_2}\prod _{i=1}^{{\widehat{n}}_1}{e}^{-\Lambda t_i^{2\gamma }}\prod _{i={\widehat{n}}_1+1}^r{e}^{-{(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma }\Lambda }{e}^{-\Lambda (n-r)t_{r:n}^{2\gamma }},\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill {\pi }_3^{\ast }\left(\gamma \right|{\lambda }_1,{\lambda }_2,\beta ,t,c)=& {\gamma }^{r+c_1-1}{e}^{-d_1\gamma }\prod _{i=1}^{{\widehat{n}}_1}{t}_i^{2\gamma -1}{e}^{-\Lambda t_i^{2\gamma }}\prod _{i={\widehat{n}}_1+1}^r{(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma -1}{e}^{-\Lambda {(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma }}{e}^{-\Lambda (n-r)t_{r:n}^{2\gamma }}\hfill \end{array}$$

(18)

and

$$\begin{array}{cc}\hfill {\pi }_4^{\ast }\left(\beta \right|{\lambda }_1,{\lambda }_2,\gamma ,t,c)=& {\beta }^{-\widehat{n_2}}\prod _{i={\widehat{n}}_1+1}^r{(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma -1}{e}^{-\Lambda {(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma }}.\hfill \end{array}$$

(19)

Since the Equations (16)–(19) cannot be computed explicitly, numerical techniques are employed. One of the most powerful numerical techniques in Bayesian estimation is the Monte Carlo Markov Chain method (MCMC). In this scenario, we suggest employing the Metropolis–Hastings (M-H) sampling method within the Gibbs algorithm, utilizing a normal proposal distribution as recommended by Tierney [26]. The procedure for Gibbs sampling incorporating the (M-H) approach is outlined as follows:

(1)

Set initial values $\left({\lambda }_1^{\left(0\right)},{\lambda }_2^{\left(0\right)},{\gamma }^{\left(0\right)},{\beta }^{\left(0\right)}\right).$

(2)

Set $j=1.$

(3)

Using the following M-H algorithm, from ${\pi }_1^{\ast }\left({\lambda }_1^{(j-1)}\right|{\lambda }_2^{(j-1)},{\gamma }^{(j-1)},{\beta }^{(j-1)},t,c){\pi }_2^{\ast }\left({\lambda }_2^{(j-1)}\right|{\lambda }_1^{\left(j\right)},{\gamma }^{(j-1)},{\beta }^{(j-1)},t,c)$, ${\pi }_3^{\ast }\left({\gamma }^{(j-1)}\right|{\lambda }_1^{\left(j\right)},{\lambda }_2^{\left(j\right)},{\beta }^{(j-1)},t,c)$, and ${\pi }_4^{\ast }\left({\beta }^{(j-1)}\right|{\lambda }_1^{\left(j\right)},{\lambda }_2^{\left(j\right)},{\gamma }^{\left(j\right)},t,c)$ generate ${\lambda }_1^{\left(j\right)},{\lambda }_2^{\left(j\right)},{\gamma }^{\left(j\right)}$, and ${\beta }^{\left(j\right)}$ with the normal proposal distributions

$$N\left({\lambda }_1^{(j-1)},var\left({\lambda }_1\right)\right),N\left({\lambda }_2^{(j-1)},var\left({\lambda }_2\right)\right),N\left({\gamma }^{(j-1)},var\left(\gamma \right)\right),\mathrm{and}N\left({\beta }^{(j-1)},var\left(\beta \right)\right),$$

and from the main diagonal in the inverse Fisher information matrix we obtained $var\left({\lambda }_1\right),var\left({\lambda }_2\right),var\left(\gamma \right)$, and $var\left(\beta \right)$.

(4)

Generate a proposal for ${\lambda }_1^{\ast }$ from $N\left({\lambda }_1^{(j-1)},var\left({\lambda }_1\right)\right)$, ${\lambda }_2^{\ast }$ from $N\left({\lambda }_2^{(j-1)},var\left({\lambda }_2\right)\right)$, ${\gamma }^{\ast }$ from $N\left({\gamma }^{(j-1)},var\left(\gamma \right)\right)$, and ${\beta }^{\ast }$ from $N\left({\beta }^{(j-1)},var\left(\beta \right)\right).$

(i)

The acceptance probabilities are

$$\left.\begin{array}{c}{\mu }_{{\lambda }_1}=\min \left[1,\frac{{\pi }_1^{\ast }\left({\lambda }_1^{\ast }\right|{\lambda }_2^{(j-1)},{\gamma }^{(j-1)},{\beta }^{(j-1)},t,c)}{{\pi }_1^{\ast }\left({\lambda }_1^{(j-1)}\right|{\lambda }_2^{(j-1)},{\gamma }^{(j-1)},{\beta }^{(j-1)},t,c)}\right],\hfill \\ \\ {\mu }_{{\lambda }_2}=\min \left[1,\frac{{\pi }_2^{\ast }\left({\lambda }_2^{\ast }\right|{\lambda }_1^{\left(j\right)},{\gamma }^{(j-1)},{\beta }^{(j-1)},t,c)}{{\pi }_2^{\ast }\left({\lambda }_2^{(j-1)}\right|{\lambda }_1^{\left(j\right)},{\gamma }^{(j-1)},{\beta }^{(j-1)},t,c)}\right],\hfill \\ \\ {\mu }_{\gamma }=\min \left[1,\frac{{\pi }_3^{\ast }\left({\gamma }^{\ast }\right|{\lambda }_1^{\left(j\right)},{\lambda }_2^{\left(j\right)},{\beta }^{(j-1)},t,c)}{{\pi }_3^{\ast }\left({\gamma }^{(j-1)}\right|{\lambda }_1^{\left(j\right)},{\lambda }_2^{\left(j\right)},{\beta }^{(j-1)},t,c)}\right]\hfill \\ \\ and\hfill \\ \\ {\mu }_{\beta }=\min \left[1,\frac{{\pi }_4^{\ast }\left({\beta }^{\ast }\right|{\lambda }_1^{\left(j\right)},{\lambda }_2^{\left(j\right)},{\gamma }^{\left(j\right)},t,c)}{{\pi }_4^{\ast }\left({\beta }^{(j-1)}\right|{\lambda }_1^{\left(j\right)},{\lambda }_2^{\left(j\right)},{\gamma }^{\left(j\right)},t,c)}\right].\hfill \end{array}\right\}$$

(ii)

From a Uniform $(0,1)$ distribution $u_1,u_2,$$u_3$, and $u_4$ are generated.

(iii)

If $u_1<{\mu }_{{\lambda }_1}$, accept the proposal and set ${\lambda }_1^{\left(j\right)}={\lambda }_1^{\ast }$, otherwise set ${\lambda }_1^{\left(j\right)}={\lambda }_1^{(j-1)}$.

(iv)

If $u_2<{\mu }_{{\lambda }_2}$, accept the proposal and set ${\lambda }_2^{\left(j\right)}={\lambda }_2^{\ast }$, otherwise set ${\lambda }_2^{\left(j\right)}={\lambda }_2^{(j-1)}$.

(v)

If $u_3<{\mu }_{\gamma }$, accept the proposal and set ${\gamma }^{\left(j\right)}={\gamma }^{\ast }$, otherwise set ${\gamma }^{\left(j\right)}={\gamma }^{(j-1)}$.

(vi)

If $u_4<{\mu }_{\beta }$, accept the proposal and set ${\beta }^{\left(j\right)}={\beta }^{\ast }$, otherwise set ${\beta }^{\left(j\right)}={\beta }^{(j-1)}.$

(5)

Set $j=j+1.$

(6)

Steps (3)–(5), are repeated N times to obtain ${\lambda }_1^{\left(j\right)},{\lambda }_2^{\left(j\right)},{\gamma }^{\left(j\right)}$, and ${\beta }^{\left(j\right)}$, j = 1, 2, …N.

To guarantee convergence and eliminate the impact of initial value selection, the first M simulated variants are eliminated. For a sufficiently high N, the chosen samples are then ${\psi }_k^{\left(j\right)},j=M+1,\dots N$. The SEL function-based approximate BEs of ${\psi }_k$ are generated using

$${\widehat{\psi }}_k^{\left(j\right)}=\frac{1}{N-M}\sum _{j=M+1}^N{\psi }^{\left(j\right)},k=1,2,3,4.$$

(20)

The approximate Bayes estimates for ${\psi }_k$, under the Entropy loss function are given as

$${\widehat{\psi }}_k^{\left(j\right)}={\left[\frac{1}{N-M}\sum _{j=M+1}^N{\left({\psi }^{\left(j\right)}\right)}^{-q}\right]}^{\frac{-1}{q}},k=1,2,3,4.$$

(21)

4. Interval Estimation

Confidence interval estimation is a fundamental statistical method used to indicate the reliability of an estimate. It provides a range of values, derived from sample data, that is likely to contain the true value of an unknown population parameter. The concept is central to inferential statistics and has numerous applications across various fields such as engineering, economics, medicine, and the social sciences. Among its key properties, the asymptotic interval is notable for its reliance on large sample sizes, where the distribution of the estimate approaches a normal distribution, making it increasingly accurate as the sample size grows. This property is particularly useful for electrical engineering projects where large data sets are analyzed for reliability and performance assessments.

Credible intervals, on the other hand, are used in Bayesian statistics and represent the range within which a parameter lies with a certain probability, given the observed data and a prior belief about the parameter’s distribution. This approach is valuable in research and development projects within electrical engineering, where prior knowledge or expert opinions can be quantitatively incorporated into the analysis, offering a more nuanced understanding of uncertainty.

Bootstrap intervals utilize resampling techniques to generate an empirical distribution of the estimator by drawing samples with replacements from the original dataset. This method does not assume a specific distribution, making it versatile and robust, especially in complex engineering studies where theoretical distributions are hard to justify. The bootstrap approach is particularly important for evaluating the uncertainty of estimates derived from small or non-standard datasets, providing a powerful tool for uncertainty quantification in both academic research and practical applications.

The application and importance of these intervals lie in their ability to quantify the uncertainty in estimates, guiding decision making and hypothesis testing. In electrical engineering, for example, they can be used to assess the reliability of system parameters, evaluate the performance of new designs, or validate models against empirical data. By understanding and applying these concepts, researchers, and practitioners can enhance the rigor and credibility of their findings, contributing to more reliable and effective solutions in their respective fields. The following subsections work out the previously mentioned interval estimations.

4.1. Asymptotic Confidence Interval

This subsection presents the observed Fisher information matrix, commonly employed for the construction of asymptotic confidence intervals (ACIs).

The MLEs $(\widehat{{\lambda }_1},\widehat{{\lambda }_2},\widehat{\gamma },\widehat{\beta })$ are approximately normal with a mean of $(\widehat{{\lambda }_1},\widehat{{\lambda }_2},\widehat{\gamma },\widehat{\beta })$ and a variance matrix $I^{-1}(\widehat{{\lambda }_1},\widehat{{\lambda }_2},\widehat{\gamma },\widehat{\beta })$. Here, $\widehat{I}\left({\lambda }_1,{\lambda }_2,\gamma ,\beta \right)$ is the observed Fisher information matrix, and it is defined as

$$\widehat{I}\left({\lambda }_1,{\lambda }_2,\gamma ,\beta \right)={\left(\begin{array}{cccc}-\frac{{\partial }^2\mathcal{l}}{\partial {\lambda }_1^2}& -\frac{{\partial }^2\mathcal{l}}{\partial {\lambda }_1\partial {\lambda }_2}& -\frac{{\partial }^2\mathcal{l}}{\partial {\lambda }_1\partial \gamma }& -\frac{{\partial }^2\mathcal{l}}{\partial {\lambda }_1\partial \beta }\\ -\frac{{\partial }^2\mathcal{l}}{\partial {\lambda }_2\partial {\lambda }_1}& -\frac{{\partial }^2\mathcal{l}}{\partial {\lambda }_2^2}& -\frac{{\partial }^2\mathcal{l}}{\partial {\lambda }_2\partial \gamma }& -\frac{{\partial }^2\mathcal{l}}{\partial {\lambda }_2\partial \beta }\\ -\frac{{\partial }^2\mathcal{l}}{\partial \gamma \partial {\lambda }_1}& -\frac{{\partial }^2\mathcal{l}}{\partial \gamma \partial {\lambda }_2}& -\frac{{\partial }^2\mathcal{l}}{\partial {\gamma }^2}& -\frac{{\partial }^2\mathcal{l}}{\partial \gamma \partial \beta }\\ -\frac{{\partial }^2\mathcal{l}}{\partial \beta \partial {\lambda }_1}& -\frac{{\partial }^2\mathcal{l}}{\partial \beta \partial {\lambda }_2}& -\frac{{\partial }^2\mathcal{l}}{\partial \beta \partial \gamma }& -\frac{{\partial }^2\mathcal{l}}{\partial {\beta }^2}\end{array}\right)}_{\left({\lambda }_1,{\lambda }_2,\gamma ,\beta \right)=\left(\widehat{{\lambda }_1},\widehat{{\lambda }_2},\widehat{\gamma },\widehat{\beta }\right)},$$

(22)

where the second partial derivatives are as follows:

$$\frac{{\partial }^2\mathcal{l}}{\partial {\gamma }^2}=\frac{-r}{{\gamma }^2}-2\Lambda \left[\sum _{i=1}^{\widehat{n_1}}{t}_i^{2\gamma }{\left(\log \left(t_i\right)\right)}^2+\sum _{i=\widehat{n_1}+1}^r{(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma }{\left(\log (\tau +{\beta }^{-1}(t_i-\tau ))\right)}^2\right]+(n-r)t_{r:n}^{2\gamma }{\left(\log \left(t_{r:n}\right)\right)}^2,$$

$$\begin{array}{cc}\hfill \frac{{\partial }^2\mathcal{l}}{\partial {\beta }^2}=& (2\gamma -1)\sum _{i=\widehat{n_1}+1}^r\frac{{\beta }^{-4}{(t_i-\tau )}^2-2{\beta }^{-3}(t_i-\tau )(\tau +{\beta }^{-1}(t_i-\tau ))}{{(\tau +{\beta }^{-1}(t_i-\tau ))}^2}-\hfill \\ & \Lambda \sum _{i=\widehat{n_1}+1}^r\gamma {\beta }^{-3}(t_i-\tau ){(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma -1}\left[2+(2\gamma -1)(t_i-\tau ){\beta }^{-1}{(\tau +{\beta }^{-1}(t_i-\tau ))}^{-1}\right],\hfill \end{array}$$

$$\frac{{\partial }^2\mathcal{l}}{\partial {\lambda }_1^2}=\frac{2\widehat{n_1}}{{\lambda }_1^2}-3{\lambda }_1^{-4}\left[\sum _{i=1}^{n_1}{t}_i^{2\gamma }+\sum _{i=n_1+1}^r{(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma }+(n-r)t_{r:n}^{2\gamma }\right],$$

$$\frac{{\partial }^2\mathcal{l}}{\partial {\lambda }_2^2}=\frac{2\widehat{n_2}}{{\lambda }_2^2}-3{\lambda }_2^{-4}\left[\sum _{i=1}^{n_1}{t}_i^{2\gamma }+\sum _{i=n_1+1}^r{(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma }+(n-r)t_{r:n}^{2\gamma }\right],$$

$$\begin{array}{cc}\hfill \frac{{\partial }^2\mathcal{l}}{\partial \gamma \partial \beta }=& -2\sum _{i=n_1+1}^r\frac{{\beta }^{-2}(t_i-\tau )}{{(\tau +{\beta }^{-1}(t_i-\tau ))}^2}+\hfill \\ & \Lambda \left[\sum _{i=n_1+1}^r{\beta }^{-2}(t_i-\tau )\left[2\gamma {(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma -1}\log (\tau +{\beta }^{-1}(t_i-\tau ))-{(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma -2}\right]\right],\hfill \end{array}$$

$$\frac{{\partial }^2\mathcal{l}}{\partial \gamma \partial {\lambda }_1}=\frac{-2}{{\lambda }_1^3}\left[\sum _{i=1}^{n_1}{t}_i^{2\gamma }\log \left(t_i\right)+\sum _{i=n_1+1}^r{(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma }\log (\tau +{\beta }^{-1}(t_i-\tau ))\right]+(n-r)t_{r:n}^{2\gamma }\log \left(t_{r:n}\right)$$

$$\frac{{\partial }^2\mathcal{l}}{\partial \gamma \partial {\lambda }_2}=\frac{-2}{{\lambda }_2^3}\left[\sum _{i=1}^{n_1}{t}_i^{2\gamma }\log \left(t_i\right)+\sum _{i=n_1+1}^r{(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma }\log (\tau +{\beta }^{-1}(t_i-\tau ))\right]+(n-r)t_{r:n}^{2\gamma }\log \left(t_{r:n}\right)$$

$$\frac{{\partial }^2\mathcal{l}}{\partial \beta \partial {\lambda }_1}=-2\gamma {\lambda }_1^{-3}\sum _{i=n_1+1}^r{\beta }^{-2}(t_i-\tau ){(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma -1},$$

$$\frac{{\partial }^2\mathcal{l}}{\partial \beta \partial {\lambda }_2}=-2\gamma {\lambda }_2^{-3}\sum _{i=n_1+1}^r{\beta }^{-2}(t_i-\tau ){(\tau +{\beta }^{-1}(t_i-\tau ))}^{2\gamma -1}$$

$$\frac{{\partial }^2\mathcal{l}}{\partial {\lambda }_1\partial {\lambda }_2}=0.$$

Consequently, the estimated asymptotic variance–covariance matrix $\widehat{V}$ for the MLEs can be obtained by taking the inverse of the observed information matrix $\widehat{I}({\lambda }_1,{\lambda }_2,\gamma ,\beta )$ which is given by

$$\left[\widehat{V}\right]={\widehat{I}}^{-1}=\left(\begin{array}{cccc}\widehat{Var}\left(\widehat{{\lambda }_1}\right)& cov(\widehat{{\lambda }_1},\widehat{{\lambda }_2})& cov(\widehat{{\lambda }_1},\widehat{\gamma })& cov(\widehat{{\lambda }_1},\widehat{\beta })\\ cov(\widehat{{\lambda }_1},\widehat{{\lambda }_2})& \widehat{Var}\left(\widehat{{\lambda }_2}\right)& cov(\widehat{{\lambda }_2},\widehat{\gamma })& cov(\widehat{{\lambda }_2},\widehat{\beta })\\ cov(\widehat{{\lambda }_1},\widehat{\gamma })& cov(\widehat{{\lambda }_2},\widehat{\gamma })& \widehat{Var}\left(\widehat{\gamma }\right)& cov(\widehat{\gamma },\widehat{\beta )}\\ cov(\widehat{{\lambda }_1},\widehat{\beta })& cov(\widehat{{\lambda }_2},\widehat{\beta })& cov(\widehat{\gamma },\widehat{\beta })& \widehat{Var}\left(\widehat{\beta }\right)\end{array}\right).$$

(23)

The $100(1-\zeta )\%$ two-sided confidence interval can be written as

$$\widehat{{\lambda }_1}\pm Z_{\frac{\zeta }{2}}\sqrt{\widehat{Var}\left(\widehat{{\lambda }_1}\right)},\widehat{{\lambda }_2}\pm Z_{\frac{\zeta }{2}}\sqrt{\widehat{Var}\left(\widehat{{\lambda }_2}\right)},\widehat{\gamma }\pm Z_{\frac{\zeta }{2}}\sqrt{\widehat{Var}\left(\widehat{\gamma }\right)},\mathrm{and}\widehat{\beta }\pm Z_{\frac{\zeta }{2}}\sqrt{\widehat{Var}\left(\widehat{\beta }\right)},$$

(24)

where $Z_{\frac{\zeta }{2}}$ is the percentile of the standard normal distribution with right-tail probability $\frac{\zeta }{2}$.

4.2. Credible Interval

Using the Metropolis–Hastings algorithm within the Gibbs sampling framework, we determined the credible confidence interval (CCI). For clarity, we refer to Section 3.2, and the algorithm steps mentioned there. Proceeding after step (6), the $100(1-\zeta )\%$ CCIs of ${\psi }_k$ where $\left({\psi }_1,{\psi }_2,{\psi }_3,{\psi }_4\right)=\left({\lambda }_1,{\lambda }_2,\gamma ,\beta \right)$ with ${\psi }_k^{\left(1\right)}<{\psi }_k^{\left(2\right)}\dots <{\psi }_k^{\left(N\right)},$ is given by

$$\left({\psi }_{k\left(N(\zeta /2\right))},{\psi }_{k\left(N\left(1-\zeta /2\right)\right)}\right).$$

4.3. Bootstrap Interval

Bootstrap confidence intervals offer a versatile approach to estimating the uncertainty of an estimator when the underlying distribution is unknown or complex. There are two main types: the bootstrap-t and the bootstrap percentile (bootstrap-p) methods.

4.3.1. Parametric Boot-p

The bootstrap percentile (p) method involves generating a large number of bootstrap samples from the original data. For each sample, the statistic of interest is calculated, creating a distribution of these statistics. The confidence interval is then directly obtained by taking percentiles from this empirical distribution. The following steps describe the algorithm of this method:

(1)

Based on $x=x_{1:n},x_{2:n},\dots ,x_{m:n},$ obtain $\widehat{{\lambda }_1},\widehat{{\lambda }_2},\widehat{\gamma }$, and $\widehat{\beta }$ by maximizing Equations (10)–(13).

(2)

Generate $x^{\ast }=x_{1:n}^{\ast },x_{2:n}^{\ast },\dots ,x_{m:n}^{\ast }$ from the PR distribution with parameters $\widehat{{\lambda }_1},\widehat{{\lambda }_2},\widehat{\gamma }$, and $\widehat{\beta }$ based on Type-II censoring under TRV, by considering the algorithm presented in [27].

(3)

Obtain the bootstrap parameter estimation ${\widehat{\psi }}_i^{\ast }=\left({\widehat{{\lambda }_1}}_i^{\ast },{\widehat{{\lambda }_2}}_i^{\ast },{\widehat{\gamma }}_i^{\ast },{\widehat{\beta }}_i^{\ast },\right)$, with $i=1,2,3,\dots ,N$ boots using the MLEs under the bootstrap sampling.

(4)

Repeat steps (2) and (3) $N$ boot times, and obtain ${\widehat{\psi }}_1^{\ast },{\widehat{\psi }}_2^{\ast },\dots ,{\widehat{\psi }}_{Nboot}^{\ast }$

(5)

Obtain ${\widehat{\psi }}_{\left(1\right)}^{\ast },{\widehat{\psi }}_{\left(2\right)}^{\ast },\dots ,{\widehat{\psi }}_{\left(Nboot\right)}^{\ast }$ by arrange ${\widehat{\psi }}_i^{\ast },i=1,2,3,\dots ,N$$boot$ in ascending orders.

Define ${\widehat{\psi }}_{boot-p}=G_1^{-1}\left(z\right)$ for a given z, where $G_1\left(z\right)=P({\widehat{\psi }}^{\ast }\le z)$ denotes the cumulative distribution function of ${\widehat{\psi }}^{\ast }.$ The $100(1-\zeta )\%$ approximate bootstrap-p CI of $\widehat{\psi }$ is given by

$$\left[{\widehat{\psi }}_{boot-p}\left(\frac{\zeta }{2}\right),{\widehat{\psi }}_{boot-p}\left(1-\frac{\zeta }{2}\right)\right].$$

(25)

4.3.2. Parametric Boot-t

The bootstrap-t method is an adaptation of the traditional t-interval, designed to handle situations where the sample size is small or the data do not meet the assumptions of normality. It involves resampling the original data with replacements to generate a large number of bootstrap samples. These are used to calculate a t-statistic, analogous to the one used in traditional t-tests but derived from the bootstrap distribution. This collection of t-statistics forms a distribution from which confidence intervals can be derived, the bootstrap-t algorithm is itemized as follows:

(1)

Repeat the initial three steps of the parametric Boot-p procedure.

(2)

Calculate the variance–covariance matrix $I^{-1\ast }$ utilizing the delta method.

(3)

Define the statistic $T^{\ast \psi }$ as

$$T^{\ast \psi }=\frac{\left({\widehat{\psi }}^{\ast }-\widehat{\psi }\right)}{\sqrt{\widehat{var\left({\widehat{\psi }}^{\ast }\right)}}}$$

(4)

Generate $T_1^{\ast \psi },T_2^{\ast \psi },\dots ,T_{Nboot}^{\ast \psi }$ from repeating steps $2-5$ N Boot times

(5)

Sort the sequence $T_{\left(1\right)}^{\ast \psi },T_{\left(2\right)}^{\ast \psi },\dots ,T_{\left(Nboot\right)}^{\ast \psi }$ by arranging ${\widehat{\psi }}_i^{\ast },i=1,2,3,\dots ,N$$boot$ in $T_1^{\ast \psi },T_2^{\ast \psi },\dots ,T_{Nboot}^{\ast \psi }$ in ascending order.

Define ${\widehat{\psi }}_{boot-t}=\widehat{\psi }+G_2^{-1}\left(z\right)\sqrt{\widehat{var\left({\widehat{\psi }}^{\ast }\right)}}$, where $G_2\left(z\right)=P(T^{\ast }\le z)$ is the cumulative distribution function of $T^{\ast }$ for a given z.

Then, the approximate bootstrap-t $100(1-\zeta )\%$ CI of $\widehat{\psi }$ is obtained by

$$\left[{\widehat{\psi }}_{boot-t}\left(\frac{\zeta }{2}\right),{\widehat{\psi }}_{boot-t}\left(1-\frac{\zeta }{2}\right)\right].$$

(26)

5. Simulation Analysis

In this section, we present various simulation methods to demonstrate the theoretical results. Initially, we create accelerated PR datasets using the inverse transformation technique. To achieve this, we employ a quantile function derived from the equation where V represents a random sample from the uniform distribution. Consequently, we generate random samples of sizes 40, and 100 using Equation (27).

$$F_j^{-1}(v,{\lambda }_1,{\lambda }_2,\gamma ,\beta )=\left\{\begin{array}{cc}\hfill & {\left[-2{\lambda }_j^{2}ln\left(1-v\right)\right]}^{2\gamma }\hfill \\ & \tau +\left(\beta {\left[-2{\lambda }_j^{2}ln\left(1-v\right)\right]}^{2\gamma }-\tau \right)\hfill \end{array}\right.$$

(27)

where $j=1,2$. Moreover, within the Type-II censoring framework, we employed two distinct predetermined numbers of failures for each sample size. Thus, we selected m = 25 and m = 35 for n = 40, and r = 75 and r = 90 for n = 100, respectively. We examined two different sets of actual parameter values in this context. In the initial approach, we set $({\lambda }_1,{\lambda }_2,\gamma ,\beta )$ = (1.5, 1.8, 1.2, 0.8), (1.5, 1.8, 1.2, 0.3), (0.6, 0.7, 2, 0.3), and (0.6, 0.7, 2, 0.8) with two distinct stress transition points: $\tau$ = 0.60 and $\tau$ = 0.90. In all scenarios, we determined the stress transition points based on the ranges of the generated samples, which varied depending on the chosen actual parameter values.

We employed the software developed by R Team et al. [28] for computational tasks. For MLE computations, we utilized the “L-BFGS-B” method within the “optim” function to optimize the profile log-likelihood function described in Equation (9) within the restricted area of $0<\beta <1$. We set the significance level to $0.05$ for approximate confidence intervals. Subsequently, we conducted simulations repeatedly for 5000 iterations. Observing that the means of the gamma priors yield the real parameter values with the given hyper-parameters. The determination of hyper-parameters relies on informative priors, which are derived from the Maximum Likelihood Estimates (MLEs) of $({\lambda }_1,{\lambda }_2,\gamma )$ by aligning the mean and variance of $({\widehat{{\lambda }_1}}^j,{\widehat{{\lambda }_2}}^j,{\widehat{\gamma }}^j)$ with those of specified priors (Gamma priors). Here, $j=1,2,3,\dots ,k$, where k denotes the number of available samples from the PR distribution. By equating the moments of $({\widehat{{\lambda }_1}}^j,{\widehat{{\lambda }_2}}^j,{\widehat{\gamma }}^j)$ with those of the gamma priors, we derive the following set of equations:

$$\begin{array}{cc}\hfill \frac{1}{k}\sum _{j=1}^k{\widehat{\gamma }}^j& =\frac{c_1}{d_1},\frac{1}{k-1}\sum _{j=1}^k{\left({\widehat{\gamma }}^j-\frac{1}{k}\sum _{j=1}^k{\widehat{\gamma }}^j\right)}^2=\frac{c_1}{d_1^2},\hfill \\ \hfill \frac{1}{k}\sum _{j=1}^k{\widehat{{\lambda }_1}}^j& =\frac{c_2}{d_2},\frac{1}{k-1}\sum _{j=1}^k{\left({\widehat{{\lambda }_1}}^j-\frac{1}{k}\sum _{j=1}^k{\widehat{{\lambda }_1}}^j\right)}^2=\frac{c_2}{d_2^2},\hfill \\ \hfill \frac{1}{k}\sum _{j=1}^k{\widehat{{\lambda }_2}}^j& =\frac{c_3}{d_3}\mathrm{and}\frac{1}{k-1}\sum _{j=1}^k{\left({\widehat{{\lambda }_2}}^j-\frac{1}{k}\sum _{j=1}^k{\widehat{{\lambda }_2}}^j\right)}^2=\frac{c_3}{d_3^2}.\hfill \end{array}$$

By solving the aforementioned system of equations, the estimated hyper-parameters can be expressed as follows:

$$\begin{array}{cc}\hfill c_2& =\frac{{\left(\frac{1}{k}{\sum }_{j=1}^k{\widehat{{\lambda }_1}}^j\right)}^2}{\frac{1}{k-1}{\sum }_{j=1}^k({\widehat{{\lambda }_1}}^j-\frac{1}{k}{\sum }_{j=1}^k{\widehat{{\lambda }_1}}^j{)}^2},d_2=\frac{\frac{1}{k}{\sum }_{j=1}^k{\widehat{{\lambda }_1}}^j}{\frac{1}{k-1}{\sum }_{j=1}^k{\left({\widehat{{\lambda }_1}}^j-\frac{1}{k}{\sum }_{j=1}^k{\widehat{{\lambda }_1}}^j\right)}^2}\hfill \\ \hfill c_3& =\frac{{\left(\frac{1}{k}{\sum }_{j=1}^k{\widehat{{\lambda }_2}}^j\right)}^2}{\frac{1}{k-1}{\sum }_{j=1}^k({\widehat{{\lambda }_2}}^j-\frac{1}{k}{\sum }_{j=1}^k{\widehat{{\lambda }_2}}^j{)}^2},d_3=\frac{\frac{1}{k}{\sum }_{j=1}^k{\widehat{{\lambda }_2}}^j}{\frac{1}{k-1}{\sum }_{j=1}^k{\left({\widehat{{\lambda }_2}}^j-\frac{1}{k}{\sum }_{j=1}^k{\widehat{{\lambda }_2}}^j\right)}^2}\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill c_1& =\frac{(\frac{1}{k}{\sum }_{j=1}^k{\widehat{\gamma }}^j{)}^2}{\frac{1}{k-1}{\sum }_{j=1}^k({\widehat{\gamma }}^j-\frac{1}{k}{\sum }_{j=1}^k{\widehat{\gamma }}^j{)}^2},d_1=\frac{\frac{1}{k}{\sum }_{j=1}^k{\widehat{\gamma }}^j}{\frac{1}{k-1}{\sum }_{j=1}^k{\left({\widehat{\gamma }}^j-\frac{1}{k}{\sum }_{j=1}^k{\widehat{\gamma }}^j\right)}^2}.\hfill \end{array}$$

We executed the MCMC algorithm 12,000 times for each of the 5000 replications. We then discarded the initial 2000 values during the burn-in period. Given that Markov chains inherently produce samples with autocorrelation, we opted for a thinning strategy, selecting every third variate to achieve uncorrelated samples from the post-burn-in sample pool. As a result, we generated 1000 uncorrelated samples from Markov chains by repeating this thinning process 5000 times.

In the simulation scenario, we present bias values and mean squared errors (MSEs) for the point estimates, along with average lengths (ALs) and corresponding coverage probabilities (CPs) of the approximate confidence intervals.

Table 1. Some simulation measures from MLE, bootstrap, Bayesian based on SELF, and ELF for parameters of PR distribution based on TRV: ${\lambda }_1=1.5, \gamma =1.2,$ ${\lambda }_2=1.8$, $\beta =0.8$.

				MLE				Bootstrap		SELF			ELF c = −1.25			ELF c = 1.25
n	$\mathit{\tau }$	m		Bias	MSE	LACI	CP	LBP	LBT	Bias	MSE	LCCI	Bias	MSE	LCCI	Bias	MSE	LCCI
40	0.6	25	${\lambda }_1$	−0.2909	0.2916	1.7842	96.1%	0.0561	0.0569	−0.0752	0.0274	0.6008	−0.0508	0.0272	0.5998	−0.0631	0.0298	0.6240
			$\gamma$	0.5299	0.6215	2.2892	95.2%	0.0704	0.0691	0.0470	0.0206	0.5408	0.0382	0.0207	0.5390	0.0264	0.0200	0.5441
			${\lambda }_2$	0.2413	0.7516	3.2657	94.8%	0.0957	0.0962	0.0232	0.0256	0.6184	0.0241	0.0256	0.6195	0.0146	0.0250	0.6224
			$\beta$	0.3789	0.5924	2.6275	95.0%	0.0799	0.0795	−0.0552	0.0169	0.4607	−0.0507	0.0167	0.4602	−0.0656	0.0190	0.4694
		35	${\lambda }_1$	0.2035	0.2412	1.3889	96.3%	0.0726	0.0730	0.0741	0.0210	0.4968	0.0468	0.0211	0.4973	0.0607	0.0198	0.4895
			$\gamma$	0.4141	0.4885	2.2082	95.4%	0.0698	0.0699	−0.0460	0.0151	0.3970	−0.0359	0.0150	0.3951	−0.0168	0.0166	0.4061
			${\lambda }_2$	−0.0342	0.4384	2.5934	95.6%	0.0805	0.0800	0.0157	0.0169	0.5044	0.0063	0.0169	0.5042	−0.0065	0.0170	0.5058
			$\beta$	0.2641	0.4174	2.4269	96.6%	0.1088	0.1088	0.0491	0.0147	0.4050	0.0410	0.0152	0.4059	0.0590	0.0179	0.3967
	0.9	25	${\lambda }_1$	−0.1263	0.1698	1.5383	95.3%	0.0493	0.0494	−0.0609	0.0244	0.5561	−0.0498	0.0242	0.5546	−0.0597	0.0261	0.5569
			$\gamma$	0.4162	0.6104	2.1593	95.4%	0.0684	0.0583	0.0461	0.0192	0.5148	0.0355	0.0200	0.5048	0.0235	0.0192	0.5524
			${\lambda }_2$	0.2353	0.6907	2.5138	95.1%	0.0777	0.0780	0.0176	0.0238	0.5957	0.0185	0.0239	0.5949	0.0094	0.0234	0.5919
			$\beta$	0.1743	0.3717	2.2913	94.7%	0.0734	0.0725	−0.0528	0.0152	0.4527	−0.0509	0.0152	0.4523	−0.0609	0.0182	0.5492
		35	${\lambda }_1$	0.1045	0.1523	1.1130	95.9%	0.0671	0.0673	0.0607	0.0150	0.4670	0.0373	0.0200	0.4675	0.0565	0.0185	0.4511
			$\gamma$	0.4030	0.4685	2.0774	95.8%	0.0587	0.0609	−0.0457	0.0146	0.3842	−0.0327	0.0136	0.3419	−0.0163	0.0160	0.4292
			${\lambda }_2$	0.1219	0.1666	1.4238	95.3%	0.0428	0.0429	0.0113	0.0147	0.4899	0.0051	0.0157	0.4901	0.0057	0.0169	0.4883
			$\beta$	0.1530	0.3199	2.1296	95.8%	0.0982	0.1079	0.0489	0.0132	0.3944	0.0393	0.0122	0.3544	0.0561	0.0159	0.4309
100	0.6	75	${\lambda }_1$	−0.0995	0.1601	1.5198	95.2%	0.0486	0.0484	0.0158	0.0182	0.5072	0.0167	0.0182	0.5077	0.0081	0.0179	0.5052
			$\gamma$	0.4028	0.3386	1.6471	94.3%	0.0541	0.0550	0.0348	0.0154	0.4495	0.0346	0.0155	0.4499	0.0236	0.0145	0.4477
			${\lambda }_2$	−0.0497	0.2607	1.9928	94.8%	0.0630	0.0643	0.0291	0.0212	0.5632	0.0172	0.0213	0.5625	0.0081	0.0206	0.5556
			$\beta$	0.4319	0.5261	2.2852	94.5%	0.0714	0.0715	0.0508	0.0152	0.4397	0.0452	0.0145	0.4405	0.0407	0.0139	0.4305
		90	${\lambda }_1$	0.0882	0.1522	1.1671	95.9%	0.0529	0.0533	0.0148	0.0124	0.4386	0.0087	0.0172	0.4391	0.0071	0.0162	0.4289
			$\gamma$	0.3253	0.2541	1.5100	94.5%	0.0508	0.0497	−0.0325	0.0118	0.3807	−0.0315	0.0117	0.3081	−0.0152	0.0124	0.3726
			${\lambda }_2$	−0.0392	0.2081	1.7523	95.0%	0.0532	0.0533	0.0039	0.0136	0.4509	0.0045	0.0136	0.4514	−0.0041	0.0137	0.4506
			$\beta$	0.4057	0.4709	1.9445	95.0%	0.0770	0.0787	0.0480	0.0132	0.3861	0.0349	0.0113	0.3287	0.0388	0.0121	0.3783
	0.9	75	${\lambda }_1$	0.0911	0.0899	1.0850	94.4%	0.0314	0.0318	0.0128	0.0174	0.4962	0.0153	0.0175	0.4961	0.0079	0.0167	0.4928
			$\gamma$	0.2835	0.1965	1.3365	95.2%	0.0422	0.0419	0.0158	0.0129	0.4219	0.0167	0.0130	0.4221	0.0075	0.0125	0.4221
			${\lambda }_2$	0.0419	0.0923	0.9197	95.2%	0.0295	0.0298	0.0274	0.0171	0.4803	0.0162	0.0172	0.4819	0.0079	0.0163	0.4717
			$\beta$	0.2108	0.1767	1.4266	95.2%	0.0444	0.0447	0.0328	0.0149	0.3453	0.0341	0.0141	0.4054	0.0207	0.0134	0.4050
		90	${\lambda }_1$	0.0810	0.0733	0.9337	94.7%	0.0416	0.0413	0.0126	0.0120	0.4317	0.0081	0.0163	0.4332	0.0071	0.0154	0.4184
			$\gamma$	0.2098	0.1822	0.9443	95.6%	0.0471	0.0468	−0.0237	0.0114	0.3597	−0.0274	0.0103	0.2994	−0.0108	0.0120	0.3628
			${\lambda }_2$	0.0317	0.0827	0.9078	95.4%	0.0279	0.0280	0.0029	0.0120	0.4102	0.0039	0.0112	0.4095	0.0039	0.0118	0.4081
			$\beta$	0.1410	0.1238	1.0821	95.9%	0.0606	0.0604	0.0413	0.0127	0.3042	0.0263	0.0102	0.3042	0.0116	0.0115	0.3401

,

Table 2. Some simulation measures from MLE, bootstrap, Bayesian based on SELF, and ELF for parameters of PR distribution based on TRV: ${\lambda }_1=1.5,\gamma =1.2$, ${\lambda }_2=1.8$, $\beta =0.3$.

				MLE				Bootstrap		SELF			ELF c = −1.25			ELF c = 1.25
n	$\tau$	m		Bias	MSE	LACI	CP	LBP	LBT	Bias	MSE	LCCI	Bias	MSE	LCCI	Bias	MSE	LCCI
40	0.6	25	${\lambda }_1$	−0.3308	0.2386	1.4870	94.7%	0.0451	0.0459	−0.0420	0.0254	0.6011	−0.0409	0.0252	0.5991	−0.0523	0.0274	0.6088
			$\gamma$	0.4901	0.4934	1.9734	95.2%	0.0623	0.0621	0.0215	0.0228	0.5649	0.0228	0.0229	0.5642	0.0101	0.0225	0.5605
			${\lambda }_2$	0.1578	0.3716	2.3093	93.5%	0.0742	0.0741	0.0267	0.0260	0.6403	0.0276	0.0261	0.6444	0.0184	0.0256	0.6265
			$\beta$	0.1282	0.0660	0.8731	94.3%	0.0284	0.0285	−0.0988	0.0034	0.2074	−0.0220	0.0034	0.2078	−0.0320	0.0039	0.2060
		35	${\lambda }_1$	0.0776	0.2268	1.0065	95.8%	0.0629	0.0632	0.0337	0.0190	0.5233	0.0344	0.0190	0.5232	0.0271	0.0184	0.5246
			$\gamma$	0.4732	0.4537	1.8195	95.5%	0.0711	0.0705	−0.0076	0.0151	0.4924	−0.0067	0.0151	0.4915	−0.0101	0.0155	0.4939
			${\lambda }_2$	−0.1569	0.3127	2.1050	95.6%	0.0707	0.0707	−0.0073	0.0177	0.5031	−0.0066	0.0177	0.5032	−0.0133	0.0181	0.5082
			$\beta$	0.1128	0.0520	0.6350	95.4%	0.0411	0.0411	0.0852	0.0021	0.2006	0.0209	0.0031	0.2006	0.0298	0.0036	0.1925
	0.9	25	${\lambda }_1$	0.3061	0.2147	1.3998	95.4%	0.0409	0.0410	−0.0416	0.0220	0.5343	−0.0384	0.0219	0.5342	−0.0496	0.0234	0.5436
			$\gamma$	0.4320	0.4638	1.8635	95.7%	0.0608	0.0606	0.0202	0.0203	0.5583	0.0213	0.0213	0.5183	0.0101	0.0226	0.5476
			${\lambda }_2$	0.1485	0.3400	1.2760	95.6%	0.0680	0.0678	0.0225	0.0261	0.6143	0.0234	0.0262	0.5061	0.0140	0.0256	0.6064
			$\beta$	0.1196	0.0611	0.8163	94.6%	0.0263	0.0262	−0.0208	0.0032	0.2053	−0.0206	0.0031	0.2001	−0.0313	0.0031	0.2048
		35	${\lambda }_1$	−0.0674	0.1845	1.0016	95.9%	0.0613	0.0628	0.0314	0.0171	0.4842	0.0304	0.0172	0.4832	0.0231	0.0164	0.4777
			$\gamma$	0.3999	0.4152	1.7359	95.8%	0.0704	0.0684	−0.0070	0.0146	0.4748	−0.0056	0.0146	0.4742	−0.0091	0.0151	0.4795
			${\lambda }_2$	0.1474	0.3048	1.1014	96.5%	0.0397	0.0397	−0.0061	0.0152	0.4716	−0.0055	0.0152	0.4713	−0.0125	0.0156	0.4710
			$\beta$	0.0734	0.0486	0.5815	96.1%	0.0370	0.0372	0.0201	0.0021	0.2005	0.0201	0.0014	0.1925	0.0283	0.0029	0.1920
100	0.6	75	${\lambda }_1$	−0.0969	0.1007	1.1855	93.9%	0.0378	0.0378	−0.0105	0.0181	0.5247	−0.0096	0.0181	0.5244	−0.0187	0.0186	0.5230
			$\gamma$	0.3355	0.2237	1.3073	95.2%	0.0399	0.0401	0.0195	0.0203	0.5065	0.0206	0.0205	0.5063	0.0094	0.0187	0.4921
			${\lambda }_2$	−0.0389	0.1561	1.2542	93.5%	0.0481	0.0485	0.0139	0.0224	0.5605	0.0148	0.0225	0.4956	0.0064	0.0222	0.5523
			$\beta$	0.1133	0.0494	0.6976	94.5%	0.0217	0.0216	0.0200	0.0017	0.1926	0.0200	0.0029	0.1933	0.0268	0.0031	0.1871
		90	${\lambda }_1$	0.0614	0.0912	0.9262	95.2%	0.0388	0.0388	0.0308	0.0162	0.4449	0.0277	0.0168	0.4452	0.0092	0.0159	0.4358
			$\gamma$	0.3313	0.2126	1.2581	95.9%	0.0407	0.0405	0.0069	0.0130	0.4336	0.0043	0.0131	0.4327	0.0089	0.0123	0.4290
			${\lambda }_2$	−0.0314	0.1277	1.0928	95.0%	0.0385	0.0384	−0.0060	0.0150	0.4520	−0.0053	0.0153	0.4511	−0.0114	0.0146	0.4598
			$\beta$	0.0692	0.0382	0.4354	94.7%	0.0240	0.0240	0.0201	0.0015	0.1872	0.0191	0.0011	0.1877	0.0195	0.0011	0.1126
	0.9	75	${\lambda }_1$	0.0911	0.0791	1.0200	94.9%	0.0334	0.0334	0.0092	0.0163	0.4957	0.0092	0.0163	0.4963	0.0084	0.0158	0.4896
			$\gamma$	0.2783	0.1912	1.1032	94.8%	0.0384	0.0384	0.0184	0.0187	0.4951	0.0192	0.0188	0.4951	0.0083	0.0178	0.4507
			${\lambda }_2$	0.0199	0.0909	0.9317	95.0%	0.0297	0.0288	0.0124	0.0185	0.5052	0.0137	0.0186	0.4504	0.0053	0.0178	0.5045
			$\beta$	0.0787	0.0246	0.5317	95.3%	0.0171	0.0170	0.0201	0.0010	0.1823	0.0193	0.0020	0.1823	0.0189	0.0031	0.1802
		90	${\lambda }_1$	0.0510	0.0691	0.6013	95.8%	0.0351	0.0344	0.0090	0.0152	0.4206	0.0082	0.0152	0.4223	0.0021	0.0149	0.4148
			$\gamma$	0.2090	0.1696	1.0311	95.2%	0.0401	0.0394	0.0024	0.0111	0.4026	0.0032	0.0111	0.4033	−0.0043	0.0109	0.4046
			${\lambda }_2$	0.0115	0.0632	0.7802	96.5%	0.0264	0.0259	−0.0058	0.0129	0.4485	−0.0051	0.0129	0.4483	−0.0051	0.0091	0.4473
			$\beta$	0.0615	0.0230	0.4061	95.7%	0.0187	0.0185	0.0200	0.0010	0.1722	0.0181	0.0010	0.1722	0.0110	0.0010	0.1522

,

Table 3. Some simulation measures from MLE, bootstrap, Bayesian based on SELF, and ELF for parameters of PR distribution based on TRV: ${\lambda }_1=0.6,\gamma =2$, ${\lambda }_2=0.7$, $\beta =0.3$.

				MLE				Bootstrap		SELF			ELF c = −1.25			ELF c = 1.25
n	$\tau$	m		Bias	MSE	LACI	CP	LBP	LBT	Bias	MSE	LCCI	Bias	MSE	LCCI	Bias	MSE	LCCI
40	0.6	25	${\lambda }_1$	−0.1998	0.1449	0.7621	96.0%	0.0327	0.0329	−0.0477	0.0182	0.3608	−0.0746	0.0152	0.3615	−0.0619	0.0134	0.3561
			$\gamma$	0.2789	0.4008	1.8532	93.8%	0.0814	0.0818	0.0447	0.0279	0.6292	0.0531	0.0279	0.6298	0.0297	0.0276	0.6339
			${\lambda }_2$	0.1927	0.1262	1.1712	95.4%	0.0514	0.0519	0.0688	0.0197	0.4464	0.0705	0.0201	0.4454	0.0532	0.0166	0.4370
			$\beta$	−0.1116	0.0305	0.3889	94.0%	0.0166	0.0164	−0.0692	0.0084	0.2207	−0.0682	0.0084	0.2224	−0.0584	0.0069	0.2084
		35	${\lambda }_1$	0.1260	0.1380	0.7400	96.4%	0.0442	0.0443	0.0465	0.0128	0.3536	0.0663	0.0130	0.3536	0.0550	0.0111	0.3459
			$\gamma$	0.2355	0.3055	1.7076	95.4%	0.0904	0.0910	0.0342	0.0202	0.5307	0.0348	0.0202	0.5309	0.0287	0.0199	0.5275
			${\lambda }_2$	0.0517	0.0590	0.8553	95.6%	0.0391	0.0388	0.0256	0.0092	0.3572	0.0267	0.0093	0.3585	0.0160	0.0086	0.3555
			$\beta$	0.0581	0.0195	0.3499	96.4%	0.0220	0.0225	0.0662	0.0075	0.2163	0.0672	0.0077	0.2162	0.0573	0.0061	0.2019
	0.9	25	${\lambda }_1$	−0.0919	0.0314	0.5944	96.6%	0.0262	0.0264	−0.0348	0.0048	0.2355	−0.0340	0.0048	0.2361	−0.0422	0.0052	0.2291
			$\gamma$	0.1929	0.2398	1.7663	95.2%	0.0754	0.0758	0.0094	0.0205	0.5689	0.0052	0.0205	0.5686	−0.0092	0.0205	0.5646
			${\lambda }_2$	0.1410	0.0346	0.4753	95.8%	0.0215	0.0213	0.0385	0.0095	0.3348	0.0398	0.0097	0.3357	0.0268	0.0081	0.3274
			$\beta$	−0.1529	0.0290	0.2949	94.2%	0.0131	0.0132	−0.0216	0.0022	0.2143	−0.0593	0.0032	0.2042	−0.0301	0.0032	0.1944
		35	${\lambda }_1$	0.0812	0.0291	0.5708	96.9%	0.0371	0.0366	0.0326	0.0041	0.2031	0.0326	0.0031	0.2131	0.0409	0.0041	0.2193
			$\gamma$	−0.1912	0.1355	1.2343	96.6%	0.0557	0.0559	0.0081	0.0147	0.4481	0.0049	0.0147	0.4485	0.0239	0.0147	0.4489
			${\lambda }_2$	0.0491	0.0325	0.2734	96.5%	0.0126	0.0126	0.0246	0.0085	0.2672	0.0206	0.0087	0.2680	0.0145	0.0072	0.2571
			$\beta$	−0.0499	0.0149	0.2827	94.6%	0.0127	0.0127	0.0204	0.0021	0.2032	0.0256	0.0031	0.2033	0.0075	0.0026	0.1831
100	0.6	75	${\lambda }_1$	−0.0888	0.0181	0.5272	96.6%	0.0238	0.0238	0.0219	0.0046	0.2130	0.0229	0.0047	0.2304	0.0130	0.0051	0.2197
			$\gamma$	0.1727	0.1880	1.3465	94.2%	0.0584	0.0584	0.0085	0.0192	0.5551	0.0051	0.0203	0.5559	0.0082	0.0202	0.5564
			${\lambda }_2$	0.0570	0.0276	0.4613	96.0%	0.0207	0.0207	0.0369	0.0091	0.2033	0.0370	0.0091	0.3304	0.0258	0.0071	0.3163
			$\beta$	0.0301	0.0059	0.2908	94.8%	0.0127	0.0127	−0.0028	0.0014	0.1484	−0.0518	0.0021	0.1482	−0.0293	0.0014	0.1462
		90	${\lambda }_1$	0.0712	0.0126	0.5264	95.4%	0.0294	0.0292	0.0203	0.0040	0.2013	0.0231	0.0030	0.2031	0.0124	0.0041	0.2030
			$\gamma$	0.1683	0.1234	1.2051	94.4%	0.0547	0.0507	0.0075	0.0138	0.4283	0.0045	0.0138	0.4383	0.0072	0.0137	0.4380
			${\lambda }_2$	0.0370	0.0242	0.2592	94.2%	0.0121	0.0122	0.0203	0.0056	0.2047	0.0202	0.0057	0.2570	0.0125	0.0050	0.2470
			$\beta$	0.0300	0.0051	0.2338	94.8%	0.0126	0.0116	0.0019	0.0012	0.1317	0.0236	0.0020	0.1349	0.0054	0.0012	0.1410
	0.9	75	${\lambda }_1$	0.0756	0.0152	0.3813	96.0%	0.0171	0.0170	0.0193	0.0043	0.2013	0.0203	0.0044	0.2015	0.0129	0.0038	0.1981
			$\gamma$	−0.0977	0.0595	0.8773	96.0%	0.0383	0.0387	0.0081	0.0164	0.4916	0.0050	0.0164	0.4919	0.0079	0.0164	0.4920
			${\lambda }_2$	0.0416	0.0274	0.2007	95.2%	0.0092	0.0092	0.0275	0.0080	0.1807	0.0356	0.0081	0.1815	0.0247	0.0067	0.1732
			$\beta$	−0.0080	0.0042	0.1561	95.2%	0.0070	0.0067	−0.0018	0.0014	0.1381	−0.0509	0.0018	0.1381	−0.0290	0.0010	0.1372
		90	${\lambda }_1$	0.0632	0.0120	0.3506	96.2%	0.0249	0.0250	0.0173	0.0040	0.1928	0.0203	0.0028	0.1985	0.0092	0.0032	0.1827
			$\gamma$	−0.0817	0.0493	0.7108	96.4%	0.0313	0.0314	0.0072	0.0129	0.4136	0.0042	0.0129	0.4371	0.0061	0.0127	0.4319
			${\lambda }_2$	0.0230	0.0206	0.1993	96.0%	0.0089	0.0088	0.0188	0.0049	0.1587	0.0189	0.0053	0.1587	0.0118	0.0046	0.1564
			$\beta$	−0.0074	0.0031	0.1471	96.2%	0.0076	0.0074	0.0010	0.0012	0.1248	0.0113	0.0013	0.1285	0.0052	0.0009	0.1284

and

Table 4. Some simulation measures from MLE, bootstrap, Bayesian based on SELF, and ELF for parameters of PR distribution based on TRV: ${\lambda }_1=0.6,\gamma =2$, ${\lambda }_2=0.7$, $\beta =0.8$.

				MLE				Bootstrap		SELF			ELF c = −1.25			ELF c = 1.25
n	$\tau$	m		Bias	MSE	LACI	CP	LBP	LBT	Bias	MSE	LCCI	Bias	MSE	LCCI	Bias	MSE	LCCI
40	0.6	25	${\lambda }_1$	−0.0942	0.0975	1.8407	94.4%	0.0379	0.0382	−0.0613	0.0136	0.3633	−0.0599	0.0134	0.3638	−0.0732	0.0149	0.3561
			$\gamma$	0.3890	0.5133	2.3610	96.2%	0.1002	0.1019	0.0160	0.0285	0.6367	0.0168	0.0285	0.6383	0.0085	0.0284	0.6491
			${\lambda }_2$	0.2001	0.1825	1.4812	95.0%	0.0679	0.0678	0.0559	0.0153	0.4141	0.0576	0.0156	0.4147	0.0404	0.0127	0.3986
			$\beta$	−0.1482	0.2855	1.1319	96.4%	0.0486	0.0491	−0.0486	0.0159	0.4415	−0.0471	0.0157	0.4391	−0.0625	0.0177	0.4456
		35	${\lambda }_1$	0.0820	0.0967	1.1768	95.0%	0.0523	0.0526	0.0593	0.0128	0.3374	0.0426	0.0128	0.3276	0.0682	0.0135	0.3458
			$\gamma$	0.3482	0.5068	2.2636	97.2%	0.1198	0.1191	0.0148	0.0188	0.5247	0.0162	0.0188	0.5251	0.0081	0.0185	0.5277
			${\lambda }_2$	0.0160	0.0805	1.1115	95.6%	0.0485	0.0488	0.0539	0.0106	0.3311	0.0550	0.0108	0.3326	0.0394	0.0092	0.3227
			$\beta$	0.1352	0.2146	1.1057	96.8%	0.0723	0.0729	0.0410	0.0150	0.3672	0.0375	0.0152	0.3695	0.0615	0.0132	0.3601
	0.9	25	${\lambda }_1$	−0.0822	0.0417	0.7350	88.2%	0.0358	0.0359	−0.0486	0.0056	0.2287	−0.0478	0.0056	0.2278	−0.0557	0.0063	0.2278
			$\gamma$	0.1871	0.2469	1.8064	95.4%	0.0837	0.0843	0.0053	0.0219	0.5513	0.0061	0.0219	0.5507	−0.0072	0.0220	0.5530
			${\lambda }_2$	0.1577	0.0430	0.5280	94.8%	0.0249	0.0249	0.0311	0.0088	0.3443	0.0325	0.0089	0.3446	0.0182	0.0077	0.3279
			$\beta$	−0.0414	0.0722	0.8377	95.7%	0.0383	0.0383	−0.0459	0.0133	0.4165	−0.0456	0.0132	0.4165	−0.0609	0.0164	0.3715
		35	${\lambda }_1$	0.0793	0.0371	0.6949	95.8%	0.0428	0.0423	0.0381	0.0052	0.2133	0.0391	0.0052	0.2033	0.0481	0.0051	0.2070
			$\gamma$	−0.1791	0.1439	1.2860	95.6%	0.0567	0.0565	0.0041	0.0169	0.5039	0.0046	0.0169	0.5032	−0.0018	0.0170	0.5084
			${\lambda }_2$	0.0104	0.0405	0.2861	95.4%	0.0128	0.0126	0.0305	0.0087	0.2515	0.0305	0.0088	0.2516	0.0141	0.0073	0.2444
			$\beta$	−0.0392	0.0681	0.7567	95.8%	0.0336	0.0339	0.0092	0.0106	0.3479	0.0106	0.0130	0.3048	−0.0533	0.0122	0.3482
100	0.6	75	${\lambda }_1$	−0.0155	0.0225	0.5852	95.4%	0.0263	0.0265	0.0223	0.0051	0.2183	0.0382	0.0054	0.2185	0.0144	0.0056	0.2177
			$\gamma$	0.1731	0.2335	1.7403	95.4%	0.0771	0.0773	0.0047	0.0214	0.5357	0.0059	0.0204	0.5371	0.0062	0.0203	0.5419
			${\lambda }_2$	0.0443	0.0370	0.4735	94.6%	0.0232	0.0239	0.0297	0.0081	0.2934	0.0317	0.0081	0.2942	0.0176	0.0071	0.2850
			$\beta$	0.0413	0.0556	0.6791	95.8%	0.0340	0.0377	−0.0133	0.0086	0.3675	−0.0412	0.0086	0.3666	−0.0523	0.0090	0.3623
		90	${\lambda }_1$	0.0147	0.0213	0.5674	95.9%	0.0323	0.0320	0.0213	0.0042	0.2086	0.0315	0.0052	0.2019	0.0141	0.0050	0.1928
			$\gamma$	0.1328	0.1320	1.1807	96.4%	0.0379	0.0478	0.0039	0.0147	0.4574	0.0043	0.0148	0.4569	0.0017	0.0144	0.4606
			${\lambda }_2$	0.0102	0.0345	0.2722	96.0%	0.0099	0.0113	0.0283	0.0079	0.2453	0.0300	0.0081	0.2503	0.0136	0.0071	0.2408
			$\beta$	0.0318	0.0510	0.6704	96.4%	0.0305	0.0305	0.0090	0.0079	0.3334	0.0102	0.0082	0.3004	0.0483	0.0081	0.3347
	0.9	75	${\lambda }_1$	0.0149	0.0184	0.4049	94.0%	0.0188	0.0188	0.0217	0.0042	0.2022	0.0280	0.0043	0.2022	0.0123	0.0038	0.2130
			$\gamma$	−0.1066	0.0595	0.8613	95.8%	0.0386	0.0386	0.0046	0.0177	0.4874	0.0054	0.0178	0.4883	0.0053	0.0171	0.4817
			${\lambda }_2$	0.0416	0.0274	0.1928	95.6%	0.0093	0.0092	0.0260	0.0060	0.1798	0.0306	0.0061	0.1802	0.0153	0.0050	0.1755
			$\beta$	−0.0396	0.0414	0.4319	95.2%	0.0189	0.0190	−0.0116	0.0083	0.3528	−0.0401	0.0083	0.3525	−0.0514	0.0074	0.3550
		90	${\lambda }_1$	0.0132	0.0172	0.3954	95.8%	0.0237	0.0232	0.0202	0.0043	0.1926	0.0250	0.0042	0.1826	0.0114	0.0032	0.1821
			$\gamma$	−0.0912	0.0495	0.7302	96.5%	0.0332	0.0323	0.0031	0.0124	0.4160	0.0040	0.0124	0.4165	0.0017	0.0123	0.4177
			${\lambda }_2$	0.0092	0.0256	0.1896	96.5%	0.0090	0.0089	0.0228	0.0048	0.1562	0.0249	0.0058	0.1565	0.0128	0.0048	0.1494
			$\beta$	−0.0265	0.0408	0.4249	95.6%	0.0208	0.0207	0.0066	0.0070	0.3079	0.0041	0.0071	0.3179	−0.0305	0.0069	0.3287

display all results from these simulation schemes. The performance of the point and interval estimations can be itemized as follows:

• Our observations consistently show reduced biases, MSEs, and ALs as sample sizes increase.
• The CPs mostly align closely with their anticipated 95% level.
• In general, the informative Bayes estimation method outperforms MLE, with the disparity between the two estimators decreasing as the sample size grows. This highlights the Bayesian methods’ advantage for smaller samples.
• In particular, confidence intervals based on the Highest Posterior Density (HPD) method tend to be smaller than those based on the Approximate Confidence Interval (ACI) method, while still providing similar CPs.
• Altering the pre-determined number of failures or stress change time yields comparable performances across all cases, demonstrating the consistent efficiency and productivity of the theoretical findings.
• Increasing the sample size generally leads to improvements in bias, MSE, and the precision of confidence intervals across all methods. This is expected because larger samples provide more information about the population. The number of bootstrap samples m also influences the Bootstrap method’s accuracy and precision, with a higher m usually leading to better estimates.
• changing the stress transition time point $\tau$ affects the estimation, especially for Bayesian estimation under ELF that adjusts based on the distribution’s tail properties. Different $\tau$ values can lead to variations in bias and MSE, suggesting the importance of choosing an appropriate $\tau$ value for accurate estimation.

6. The Optimal Stress Change Time and Sensitivity Analysis

In this section, we describe an optimal method based on asymptotic variances in maximum likelihood estimators. The inverse Fisher information matrix’s diagonals can be used to compute the parameters’ asymptotic variances. In this section, we used the sum of coefficients of variations (SVCs) as the optimal function instead of the sum of parameter variances, as recommended and implemented by Samanta et al. ([29,30]). Samanta et al. [29] proposed a method to calculate an optimal solution by minimizing the predicted value of the SVC. Since the sum of variances can be calculated using the variance of any specific parameter if the parameter values are on a different scale. That is why we employed the expected value of the SVC by maximizing $E\left(\varphi \right(\tau \left)\right)$, where

$$\varphi \left(\tau \right)=\frac{\sqrt{F_{11}^{-1}}}{\widehat{{\lambda }_1}}+\frac{\sqrt{F_{22}^{-1}}}{\widehat{{\lambda }_2}}+\frac{\sqrt{F_{33}^{-1}}}{\widehat{\gamma }}+\frac{\sqrt{F_{44}^{-1}}}{\widehat{\beta }},$$

(29)

where $F_{ii}^{-1}$ is the element in the main diagonal of the inverse Fisher information matrix that was described by Equation (22). However, the closed forms of the parameters’ posterior variances may be imprecisely estimated. Samanta et al. [30] recommend adopting the Gibbs sampling technique for computation.

• Step 1: Obtain the samples $U_1$, $U_2$ and $U=\min \{U_1,U_2\}$ using given $\tau$, n, r and parameter values.
• Step 2: The objective function $\varphi \left(\tau \right)$ is calculated.
• Step 3: For N times, repeat Step 1 to Step 2, and obtain ${\varphi }^1\left(\tau \right),{\varphi }^2\left(\tau \right),\dots ,{\varphi }^N\left(\tau \right)$.
• Step 4: The median of the objective functions is obtained and applied to ${\varphi }^m\left(\tau \right)$.
• Step 5: For all possible values of $\tau$ repeat Step 1 to Step 4.
• Step 6: The optimal $\tau$ for which ${\varphi }^m\left(\tau \right)$ is the minimum is obtained.

Optimal stress change time $\tau$ values, indicated by ${\tau }^{\ast }$ are determined for given $n,r$, and ${\psi }_i$ for $i=1,\dots ,4$ and are reported in

Table 5. Optimal stress change time $\tau$ for different sample sizes and parameter values by SVC $\varphi \left(\tau \right)$.

n	$\mathit{\tau }$	m	Table 1	Table 2	Table 3	Table 4
40	0.6	25	0.4225	0.3565	0.5263	0.5119
		35	0.3850	0.3477	0.3439	0.3934
	0.9	25	0.4215	0.3301	0.2896	0.2695
		35	0.3384	0.3605	0.2680	0.2661
100	0.6	75	0.3600	0.2797	0.1999	0.2476
		90	0.3199	0.2737	0.1683	0.2312
	0.9	75	0.2850	0.2638	0.1686	0.2050
		90	0.2349	0.2352	0.1567	0.2008

.

From Table 5, it is evident that the optimal stress change times, denoted as $\tau$, fall within the range of 0.6 to 0.9 for the first parameter set. As the range of the generated dataset is not extensive, there is not a significant deviation in the range of $\tau$ in this initial case. It is noticeable that the stress change times utilized in the simulations closely align with the optimal stress change times. Hence, the consistency and effectiveness of the simulation outcomes are contingent upon accurately determining the stress change time.

7. Real Data Examples

In this section, two real data sets are examined for the suitability of the PR model with tampered random variables under the Type-II censoring framework.

7.1. HIV Infection to AIDS

This example discusses the application of a real-life dataset focusing on male individuals and their progression from HIV infections to AIDS over nearly 15 years. According to the United States Center for Disease Control and Prevention, 54% of total diagnosed adult AIDS cases in the U.S. up to 1996 were due to intimate contact with a person who was HIV positive, also, an additional 40% of incident cases occurred in that same year. A subset of the 54% who also engaged in injection drug use accounted for an additional 7% of cumulative and 5% of incident cases in 1996. These data were collected during the era of antiretroviral combined therapy in 1996. For further background information, readers are directed to studies by Dukers et al. [31] and Xiridou et al. [32], while Putter et al. [33] and Geskus et al. [34] cite this dataset as an example for competing risk analysis. The dataset encompasses instances where some patients either remained uninfected or their outcomes were censored in the study.

We focused on a pre-determined number of failures, setting r as 150 from a complete dataset of $n=222$. We also examined stress change times: $\tau =4.6$. For clarity, we present the competing risk data as follows in

Table 6. Data from HIV Infection to AIDS dataset.

$t_i$	c	$t_i$	c	$t_i$	c	$t_i$	c	$t_i$	c	$t_i$	c	$t_i$	c	$t_i$	c	$t_i$	c	$t_i$	c
0.112	1	2.048	1	2.798	1	3.373	0	3.8	0	4.389	1	5.018	0	5.566	0	5.982	1	6.461	0
0.137	1	2.053	1	2.814	1	3.439	1	3.817	1	4.394	1	5.021	1	5.574	0	6.018	1	6.511	1
0.474	1	2.155	0	2.866	1	3.477	0	3.819	0	4.4	1	5.082	1	5.582	0	6.042	0	6.516	1
0.824	1	2.177	0	2.875	0	3.477	1	3.88	1	4.52	1	5.106	1	5.618	0	6.042	1	6.579	0
0.884	1	2.234	0	2.891	0	3.486	0	3.94	1	4.523	1	5.12	1	5.667	0	6.045	0	6.733	0
0.969	1	2.283	0	2.982	1	3.513	0	3.953	0	4.583	0	5.224	1	5.678	0	6.054	0	6.801	0
1.013	1	2.322	1	3.039	1	3.535	0	3.975	0	4.608	0	5.251	0	5.7	1	6.177	0	6.82	1
1.101	1	2.513	1	3.064	0	3.584	1	4.033	1	4.69	0	5.314	1	5.703	1	6.195	0	6.85	1
1.205	1	2.533	0	3.064	1	3.592	0	4.079	1	4.734	1	5.336	1	5.723	0	6.199	0	6.866	0
1.44	0	2.565	1	3.195	0	3.639	0	4.099	0	4.811	0	5.374	1	5.73	1	6.218	1	6.943	1
1.462	1	2.571	1	3.214	0	3.647	0	4.219	0	4.854	1	5.454	1	5.736	1	6.224	0	6.955	0
1.503	1	2.631	1	3.22	1	3.663	0	4.219	0	4.909	1	5.478	1	5.886	1	6.267	0	6.979	1
1.593	1	2.672	0	3.242	0	3.707	0	4.23	1	4.966	1	5.525	1	5.889	1	6.311	1	7.006	0
1.837	0	2.683	0	3.258	1	3.724	0	4.334	1	4.981	0	5.555	0	5.908	1	6.412	0	7.17	1
1.889	1	2.705	0	3.315	0	3.797	1	4.375	1	5.013	0	5.563	1	5.938	0	6.439	1	7.302	0

, where the black color is $t_i<\tau$ and the gold color is $t_i>\tau$.

Table 7. MLE and different measures of HIV Infection to AIDS data.

${\mathit{t}}_{\mathit{i}}<\mathit{\tau }$	Size				KSD	PVKS	AIC	BIC	CAIC	HQIC
$t_i<max\left(t\right)$	220	$\lambda$	5.6859	0.9111	0.0890	0.1853	602.8644	608.8856	602.9460	605.3106
		$\gamma$	1.3300	0.0913
$t_i<\tau$	121	$\lambda$	3.9018	0.7312	0.1217	0.1813	260.7185	265.5074	260.8724	262.6399
		$\gamma$	1.4365	0.1381
$\tau <t_i<max\left(t\right)$	101	$\lambda$	44.2060	20.9781	0.0612	0.9585	145.3549	152.0572	145.7241	148.0139
		$\gamma$	1.1368	0.1974
		$\beta$	0.0432	0.0353

showcases the MLE alongside various fit metrics for the HIV Infection to AIDS dataset, utilizing both the baseline model and SSLT as complete datasets. The analysis derived from Table 7 indicates an adequate fit of the model to the data, evidenced by a Kolmogorov–Smirnov P-value (KSPV) exceeding $0.05$. Furthermore, the table provides a range of fit indices, including the Consistent Akaike Information Criterion (CAIC), Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), and Hannan–Quinn Information Criterion (HQIC), all of which serve as measures of goodness-of-fit.

Table 8. MLE and Bayesian estimation for the parameters of PR based on SSLT under TRV.

	MLE				Bayesian
	Estimates	StEr	Lower	Upper	Estimates	StEr	Lower	Upper
${\lambda }_1$	6.0925	1.0225	4.0884	8.0966	6.1924	0.7669	4.6570	7.6641
$\gamma$	1.1113	0.1052	0.9050	1.3175	1.1176	0.0774	0.9679	1.2601
${\lambda }_2$	6.6004	1.1181	4.4089	8.7919	6.7538	0.8359	5.3139	8.4690
$\beta$	0.5539	0.0991	0.3596	0.7482	0.5801	0.1257	0.3312	0.8086
$1-F_1(\overline{t};\gamma ,{\lambda }_1)$	0.70965				0.71305
$1-F_2(\overline{t};\gamma ,{\lambda }_2)$	0.74660				0.75253
$1-F_1(\overline{t};\gamma ,{\lambda }_1,\beta )$	0.74056				0.74077
$1-F_2(\overline{t};\gamma ,{\lambda }_2,\beta )$	0.77422				0.77705
$1-\left(1-F_1\left(t_r\right)\right)\left(1-F_2\left(t_r\right)\right)$	0.97519				0.97895

presents the maximum likelihood and Bayesian point estimation in addition to the interval estimates for the PR parameters derived from step-stress life testing using the Tampered Random Variable model. Table 8 presents a reliability analysis that evaluates the reliability function of various models through maximum likelihood and Bayesian methods for estimating parameters. The models analyzed include those with a risk factor from cause I, from cause II, and both under standard conditions, followed by an examination under an accelerated framework. Additionally, the reliability of the TRV model is analyzed in the context of two competing risk factors. The findings suggest that the TRV model exhibits the greatest reliability among the models assessed, underscoring the robustness of our proposed model. Figure 1 depicts the likelihood profile for the PR parameters based on SSLT under the TRV model which indicates the existence of the maximum value for the log-likelihood function. Figure 2 illustrates the trace plots and marginal posterior probability density functions of the parameters for the PR distribution, employing SSLT under the TRV model, as obtained via Bayesian estimation.

7.2. Electrical Appliances Data

The real-world dataset analyzed in reference [35] (p. 441) examines 36 small electronic components subjected to an automated life test, where failures are categorized into 18 types. However, out of the 33 identified failures, only seven modes were observed, with modes 6–9 recurring more than twice. Mode 9 failure is particularly significant. Consequently, the dataset is categorized into two failure causes, c = 0 (mode 9 failure) and c = 1 (all other modes). The provided data presents the failure times in sequence along with the respective cause of each failure, the stress change time is selected to be $\tau =2500$ as detailed in

Table 9. Electrical appliances data.

$t_i$	c	$t_i$	c	$t_i$	c	$t_i$	c	$t_i$	c
11	1	381	1	1594	1	2400	0	2694	0
35	1	708	1	1925	0	2451	1	2702	1
49	1	958	1	1990	0	2471	0	2761	1
170	1	1062	1	2223	0	2551	0	2831	1
329	1	1167	0	2327	1	2568	0	3034	0

.

Table 10. MLE and different measures for electrical appliances data.

${\mathit{x}}_{\mathit{i}}<\mathit{\tau }$	Size		Estimate	StEr	KSD	PVKS	AIC	BIC	CAIC	HQIC
$x_i<max\left(x_i\right)$	25	$\lambda$	53.3420	43.2098	0.2403	0.0938	423.9574	426.3951	424.5028	424.6335
		$\gamma$	0.5803	0.1044
$x_i<\tau$	18	$\lambda$	22.3780	16.9145	0.1772	0.5644	296.2920	298.0727	297.0920	296.5375
		$\gamma$	0.4859	0.1002
$\tau <x_i<max\left(x_i\right)$	7	$\lambda$	2814.8644	5262.1741	0.2066	0.8722	94.8192	94.6569	102.8192	92.8135
		$\gamma$	0.7401	0.1678
		$\beta$	0.0036	0.0012

discusses MLE and different measures used for the electrical appliances data in baseline model and SSLT model as complete data. From the results in Table 10, we note that the data are fitting of this model where the KSPV is greater than 0.05. Also, some different measures have been obtained as CAIC, AIC, BIC goodness-of-ft measures, and HQIC.

Table 11. MLE and Bayesian estimation for the parameters of PR based on SSLT under TRV: electrical appliances data.

	MLE				Bayesian
	Estimates	StEr	Lower	Upper	Estimates	StEr	Lower	Upper
${\lambda }_1$	30.436	11.750	7.407	73.065	30.475	2.933	24.635	36.016
$\gamma$	0.456	0.091	0.277	0.635	0.453	0.022	0.412	0.495
${\lambda }_2$	37.277	16.856	4.240	89.915	37.313	3.475	31.022	44.277
$\beta$	0.086	0.043	0.0015	0.179	0.096	0.042	0.020	0.171
$1-F_1(t_r;\gamma ,{\lambda }_1)$	0.44465				0.46349
$1-F_2(t_r;\gamma ,{\lambda }_2)$	0.58258				0.59872
$1-F_1(t_r;\gamma ,{\lambda }_1,\beta )$	0.12041				0.15590
$1-F_2(t_r;\gamma ,{\lambda }_2,\beta )$	0.24386				0.28945
$1-\left(1-F_1\left(t_r\right)\right)\left(1-F_2\left(t_r\right)\right)$	0.97064				0.97549

presents the maximum likelihood and Bayesian point estimation in addition to the interval estimates for the PR parameters derived from step-stress life testing using the Tampered Random Variable model for the electrical appliances data. Similar to the discussion of reliability analysis in the first data set in Table 8, the reliability analysis presented in Table 11 indicates that the TRV model outperformed the other models. Figure 3 depicts the likelihood profile of PR parameters based on SSLT under the TRV model for electrical appliances data. From Figure 3, we can conclude that the parameters of PR distribution based on SSLT under TRV have maximum value for the log-likelihood function for electrical appliances data. Figure 4 shows the trace plots and marginal posterior probability density functions of the parameters for the PR distribution based on SSLT under the TRV model, derived through Bayesian estimation for the electrical appliances data.

8. Conclusions

In conclusion, this work has significantly contributed to the field of reliability engineering through the application of the Tampered Random Variable (TRV) model within the step-stress life testing (SSLT) framework, particularly focusing on the Power Rayleigh distribution in the context of competing risks. By integrating TRV with SSLT under such complex scenarios, the study has addressed critical gaps in current research, particularly the various applications of TRV modeling in competing risk analyses.

The methodological advancements presented in this paper, including the use of maximum likelihood estimation and the Bayesian methods for inferential analysis, as well as Monte Carlo simulations for estimator performance evaluation, represent a robust approach to understanding and improving product reliability under varied stress conditions. These techniques have been validated through empirical analysis of real-world datasets from the medical sector, regarding AIDS infection, and the electrical engineering domain, focusing on electronic component failures. The reliability evaluations underscore the model’s empirical suitability and the potential for broader application.

Furthermore, the study’s exploration of Type-II censoring schemes as a solution to information shortage in lifetime experiments highlights the practical value of the research, offering other options for more cost-effective and efficient testing methodologies. The comparison of TRV modeling with other established models (CEM and TFR) within a competing risks framework not only clarifies the conditions under which these models converge but also showcases the unique advantages of TRV in handling complex, multi-step-stress situations and discrete or multivariate lifetime data.

The comprehensive analysis and the resulting insights into model precision, reliability, and risk management presented in this study provide a solid foundation for future research in this area. It opens up new ways for the development of more accurate and dependable models, enhancing the decision making process and risk management strategies in the medical, industrial, and mechanical domains.