Reliability analysis of constant partially accelerated life tests under progressive first failure type-II censored data from Lomax model: EM and MCMC algorithms

1.   Introduction

The Pareto system of distributions includes four types (I–IV) of cumulative distribution functions that produce different density shapes, see Johnson et al. ^[1]. The type II Pareto distribution is known in the literature by the Lomax distribution, see Lomax ^[2]. Lomax distribution is considered a serious distribution of lifetime models because it belongs to the family of decreasing failure rate models. When the population has a heavy-tailed distribution and the experimenter seeks a good probability distribution describing this case, the Lomax model is an excellent alternative to most common models such as the exponential, Weibull, or gamma distributions, see Bryson ^[3]. This distribution could be a flexible model because it's used for analyzing different lifetime data. It has various applications in medicine, biological sciences, engineering, operations research, natural sciences, queuing theory, and internet traffic modeling. One can see, among many others, Johnson et al. ^[1]. Recently, David et al. ^[4] found many applications for this distribution in business, economics, and actuarial modelling. Additionally, it's noted that this distribution has some characteristics that make it a preferred distribution. One in every of these, the statistical and reliability properties of this distribution may be expressed in closed forms, which help statisticians model different types of data.

Due to these features, several authors used this model in modeling data, especially, censored observations. Hence our motivation to verify this distribution throughout this paper is because of its practical importance in the many different alternative areas mentioned above, as mentioned in many references. Among these references, Harris ^[5] and Atkinson and Harrison ^[6] found this distribution very convenient when examining data on income and wealth. Dubey ^[7] reported that this distribution is one of the special cases of the specific compound gamma distribution. Bryson ^[3] suggested that when the available data are heavy-tailed, the Lomax distribution is another suitable distribution to the exponential distribution. Tadikamalla ^[8] listed that the Lomax distribution can be derived from a family of Burr distributions. In addition to any or all of the above, this distribution has itself been used as a basis for several generalizations, see Al-Awadhi and Ghitany ^[9], Ghitany et al. ^[10], and Punathumparambath ^[11].

If $X$ is a continuous random variable that follows a Lomax distribution with shape and scale parameters $\theta$ and $\beta$, respectively, then the probability density function (PDF) and its corresponding cumulative distribution function (CDF) can be formulated as

and

respectively. Further, the Lomax failure rate function can be expressed as

where$\ \beta > 0$ and$\ \theta > 0.\ $In industrial operations, under normal operating conditions, we need a very long time to obtain the failure times of these products, which also causes the average failure times of these products to be long, this is not in line with the industrial revolution and modern technology. To provide quick solutions in such cases, experimenters perform the ALT, where the units we want to test, are placed under several stress levels whose values are higher than the normal stress values in order to reduce the time it takes to fail. Under typical conditions, the experimenters use the data obtained from the accelerated test to estimate the failure distribution of the units.

The types of stress loading in ALT are generally classified as constant stress, step stress, or progressive stress. The constant-stress loading is a time-independent test setting where the stress remains fixed until an item is desert. It has different advantages when compared to time-dependent stress loads. This is simply, because most real products operate under constant stress. There are several major references in the area of ALT, among the most important of these references, the follower can look at Nelson ^[12], Meeker and Escobar ^[13], Bagdonavicius and Nikulin ^[14], among others.

ALT first assumption is that there is a known or assumed relationship between life stresses, therefore, data generated from accelerated conditions can be relied upon for normal use conditions. However, it has been found that in some cases this relationship cannot be known or even adequately assumed, especially when the units under test are new. Therefore, partially accelerated life tests (PALT) are oftentimes used in such cases. In general, PALT can be divided into two types, the first is constant-stress PALT and the other is step-stress PALT. In constant-stress PALT, all groups of test units are placed separately in use conditions and accelerated conditions. For step-stress PALT, the test conditions for the remaining items in the experiment shift from conditions of use to conditions of higher stress at a given time or when a specified number of failures occur.

Constant-stress PALT is our area of interest in this paper. In this type, all test units are divided into two groups, where the first group is assigned to work under normal conditions while the second group is assigned to work under accelerated conditions. In recent years, a lot of studies have been proposed in the field of PALT. For instance, we can refer to Ismail ^[15], Bhattacharyya and Soejoeti ^[16], Gouno ^[17], El-Morshedy et al. ^[18], Nassar et al. ^[19], Rao ^[20], Bai and Chung ^[21], Hassan and Al-Ghamdi ^[22], Wu and Kus ^[23], Wu et al. ^[24], Fan et al. ^[25], and Lio and Tsai ^[26]. Although the main objective of PALT is to shorten the test time of the experiment, the experimenter loses a lot of time waiting for all test units to fail. For this, it is necessary to deal with censored data. In literature, type-I and type-II censoring schemes are the main censoring schemes. These do not permit the withdrawal of intermediate components. Because of this constraint, these schemes are not flexible enough. To avoid the difficulties of these schemes, various generalizations have been proposed. One of these is the progressive censoring scheme. It allows for live units to be removed from the test at time points (controlled conditions) other than the final termination point of the experiment. The article presented by Balakrishnan and Aggarwala ^[27] has many applications in progressive type-II censoring. Sometimes, some problems arise concerning the age of the products, and therefore the experimental time of control is simply too long. To resolve this issue, the statistician provided other censoring schemes. Among them was Johnson ^[28], who introduced the scheme of the first failure censoring.

During this censoring, the test units whose number is $N = n\times k$ are distributed over the $n$ number of groups where each group has $k$ number of test units, after that, a life test is performed for all groups at the same time and conditions and the first failure times are recorded in each of the test groups. With this life-telling design, a great deal of time and cost is saved. Thus, it is clear that both first failure and first failure progressive censoring schemes significantly improve the efficiency of the life test. From this point, Wu and Kus ^[23] presented the progressive first-failure censoring scheme which profits the practical advantages of the two schemes. There are several advantages that distinguish this censoring, including the test time is short as well as the resources are saved, and a few units at risk will be removed during the experiment at each of several ordered failure times. Another advantage of this censoring is that it allows the experimenters to remove one or more of the groups under test even without noticing the first failure in the groups that will be removed from the experiment. Hence, this censoring is more efficient in reliability studies. Recently, more literature is available on progressive first-failure censoring, including Lio and Tsai ^[26], Soliman et al. ^[29], Ahmed ^[30], Krishna et al. ^[31], Kumar et al. ^[32], El-Din et al. ^[33], Ahmadi and Doostparast ^[34], and Saini et al. ^[35].

Now, we briefly discuss the first-failure progressive type-II censoring scheme. Let $n$ independent groups with $k$ items in each group are put on a life test. Start dropping $R_{1}$ number of groups as soon as the first failure time $X_{1:m:n:k}$ has occurred. When the second failure time $X_{2:m:n:k}$ occurs, randomly delete $R_{2}$ groups from the experiment and also delete the group in which the second failure was observed and repeat. The experimenter continues in the same manner until all remaining live $R_{m}$ groups and the group in which the $m$th failure has occurred are removed at the time of the $m$th failure. The observed failures $X_{1:m:n:k} < X_{2:m:n:k} < ... < X_{m:m:n:k}$ called progressive first failure censored order statistics and while $\Re = (R_{1}, ..., R_{m})$ is called a progressive censoring scheme. This scheme includes various kinds of censoring schemes. For example, if $k = 1$, it reduces to a progressive censoring scheme. Also, it contains the usual order statistics ($k = 1,$ $\Re = (0, 0, ..., 0)$ and $m$ $= n$) and type II right censoring scheme ($k = 1, R_{1}, ..., R_{m-1} = 0$ and $R_{m} = n-m$) as special cases.

Before completing this section, a brief presentation on the parameter estimation of the Lomax distribution should be made. There are several studies presented by many authors focusing on estimating the parameters of the Lomax distribution when observations are completed or censored. The Bayesian and classical estimators for the same sample size from Lomax distribution based on the progressive type-I censoring, are discussed by Elfattah et al. ^[36]. Raqab et al. ^[37] presented different predictors of failure times based on multistage progressive censoring from Pareto distribution. The optimal censoring scheme for estimating the parameters of the Lomax distribution based on progressive type-II censoring is presented by Cramer and Schmiedt ^[38]. Empirical Bayes estimators of reliability performances using progressive type-II censoring from the Lomax model are discussed by El-Din et al. ^[33]. Al-Zahrani and Al-Sobhi ^[39] have concerned with the estimation problem of the $P\left(Y < X\right)$ of the Lomax distribution based on general progressive censoring. When the available data is progressive type-II censoring data, the MLEs of the Lomax distribution are derived using the expectation-maximization (EM) algorithm by Helu et al. ^[40]. Wei et al. ^[41] studied Bayes estimation of the Lomax distribution parameters in the composite LINEX loss of symmetry. Furthermore, based on the progressively type-I hybrid censored sample, Asl et al. ^[42] applied the classical and Bayesian inferential procedures for the Lomax distribution. Chandra and Khan ^[43] have derived a step-stress ALTs model for the Lomax failure time under modified progressive type-I censoring. Recently, under generalized progressive hybrid censoring, Mahto et al. ^[44] analyzed a partially observed competing risk model for the Lomax model.

Meanwhile, no articles appeared on estimating parameters of the Lomax distribution based on a constant-stress PALT under progressive first failure type-II censored data. Our aim in this article is to address this research point. Two approaches will be used to estimate the parameters of the Lomax distribution, the maximum likelihood (ML) and Bayesian methods. When calculating the maximum likelihood estimates (MLEs) and their associated confidence intervals, the Newton-Raphson (NR) and EM algorithms, are derived and discussed in detail. On the other hand, when the Bayes procedure is performed, there is difficulty in obtaining the posterior distributions because this process often requires the computation of integrals, which is often difficult to compute, whether using complex high or low-dimensional models. In such a case, we can rely on the MCMC to estimate the unknown parameters. Assuming independent gamma priors for the parameters, the Bayes estimates will be performed under different loss functions. Based on various techniques of the MCMC, the Gibbs within Metropolis samples can be applied to generate samples from the posterior distributions. The performance of the Bayes estimators can be compared with classical MLEs through extensive computer studies. Approximate 95% confidence intervals of the unknown parameters can be calculated based on the MLEs via NR, EM, and MCMC algorithms. Also, we will compare them through their average lengths and coverage probability.

As for the form and organization of this paper as well as its sequence, it is as follows: Model description and basic assumptions are presented in Section 2. Frequentest estimations including MLEs and approximate confidence intervals of unknown parameters, are provided in Section 3. During this section, we also developed the EM algorithm for calculating the MLEs and the observed Fisher information matrix. Bayes point estimates via (symmetric-asymmetric) loss functions and their corresponding credible intervals, are constructed in Section 4. A real-life example as well as simulation studies are presented and discussed in Section 5. Finally, a few concluding remarks are reported in Section 6.

2.   Model descriptions and assumptions

2.1.   Model description

As explained above, in constant-stress PALT, $N_{1}$ number of test units are randomly selected from all available units ($N$), and are operated under normal conditions, while the remaining $N_{2} = N-N_{1}$ number of units are operated under accelerated conditions. Assume the lifetime of the test unit follows the Lomax distribution, therefore, the PDF, CDF, and failure rate function at normal conditions are given, respectively by

and

where$\ \beta > 0$ and$\ \theta > 0.\ $The hazard rate function for any item tested at the accelerated condition is given as follows

where $\lambda$ ($\lambda > 1$) is the acceleration factor. Under the accelerated condition, the PDF and its corresponding CDF are obtained, respectively by

and

Here, the progressive first-failure censoring is combined with constant-stress PALT. Hence, the total of $N$ test units will be divided into two groups, namely, group I and group II. Elements of the first group ($N_{1} = n_{1}k_{1}$) are assigned to normal conditions whereas the items of the second group ($N_{2} = n_{2}k_{2}$) are allocated to stress conditions. Each group is divided under normal or accelerated conditions into several groups containing the number of test units $k_{d}$. In this scheme, $R_{1i}$ and $R_{2i}$ are the progressive censoring schemes for normal and accelerated tests, respectively. This scheme runs until $m_{d}$; $d = 1, 2$ failures are observed in each test condition, respectively, as follows.

For group I (items dedicated to normal condition) or group II (items are dedicated to stress condition), suppose that a random sample $X_{d} = (X_{d1}, X_{d2}, ..., X_{n_{d}k_{d}})$ of size $N_{d} = n_{d}k_{d},$; $d = 1, 2$ is put on a life test experiment where the lifetimes $X_{di}^{^{\prime}}$s are independent and identically distributed (IID) random variables with PDF $f_{X_{d}}(x_{d}; \Theta)$ and CDF $F_{X_{d}}(x_{d}; \Theta)$ where $\Theta$ is a vector of unknown parameters, and $N_{d}$ test units are divided into $n_{d}$ independent groups with $k_{d}$ items in each group. Now, when the first failure $X_{d1}$ is appeared, then the group with observed failure and also $R_{d1}$ of the remaining ($n_{d}-1$) groups are randomly removed from the experiment. Next, when the second failure $X_{d2}$ is observed, then the group with the second observation and also $R_{d2}$ of the remaining ($n_{d}-R_{d1}-1$) groups are randomly removed. The experiment continues until the $m_{d}$th failure with lifetime $X_{dm_{d}}$ is observed and the group with this observed failure and the remaining $R_{dm_{d}}$ groups are all removed from the experiment. In this way, a progressive first-failure censored sample $X_{d} = (X_{d1}, X_{d2}, ..., X_{dm_{d}})$ of size $m_{d}$ can be obtained from a parent sample of size $N_{d}$ with censoring scheme $\Re_{d} = (R_{d1}, R_{d2}, ..., R_{dm_{d}})$ where ${ \sum\limits_{i = 1}^{m_{d}}} R_{di} = n_{d}-m_{d}$. In our study, $R_{di}, i = 1, 2, ..., m_{d}$ are fixed prior and $m_{d} < n_{d}$. With a progressively first-failure type-II censoring scheme under constant-stress PALT, the likelihood function of the observed sample $X_{d}$ can be written as

2.2.   Basic assumptions

Regarding the proposed PALT approach, the following assumptions are made:

$(\textbf{A1})$ The total number of units under test is $N = n_{1} k_{1}+n_{2}k_{2}$.

$(\textbf{A2})$ The lifetime of all units tested under various normal or accelerated conditions follow the Lomax distribution.

$(\textbf{A3})$ The lifetimes of test items are IID random variables.

$(\textbf{A4})$ The lifetimes $X_{1i}$, $i = 1, 2, \ldots, m_{1}$ of items assigned to the normal condition, while the lifetimes $X_{2i}$, $i = 1, 2, \ldots, m_{2}$ of items assigned to the accelerated condition are mutually independent.

$(\textbf{A5})$ The lifetime of any item at accelerated condition is $X_{2} = \lambda^{-1}X_{1}$.

Next based on a progressive first failure censored sample from the Lomax distribution, the MLEs via both NR and EM algorithms are derived for the unknown parameters $\theta$, $\beta, \ $and$\ \lambda$.

3.   Maximum likelihood estimation and confidence interval

3.1.   Maximum likelihood estimation via Newton-Raphson method

Consistency, asymptotic efficiency, asymptotic unbiased, and asymptotic normality are some of the characteristics of the MLE, which made it one of the most important and popular methods for fitting the statistical model. In this section, the MLEs of $\theta$, $\beta, \ $and$\ \lambda$ can be derived from the likelihood function (LF) presented in (7). The substitution of (1)–(6) into (7) yields

where

The MLEs of $\theta$, $\beta\ $and$\ \lambda$ are the values that maximize the LF in (8). Maximizing the LF is difficult, so it is preferable to maximize the logarithm of this function, where the log-LF, say $\ell = \log L\left(\theta, \beta, \lambda|y\right)$, can be formulated as

By calculating the first derivatives of (10) with respect to $\theta$, $\beta$, and $\lambda$ and then setting them equal to zero, the resulting simultaneous equations are expressed as follows

and

The preceding equations mathematically represent a system of three nonlinear equations in three unknowns $\theta$, $\beta$, $\ $and$\ \lambda$. Theoretically, there is difficulty in providing closed-form solutions to the previous nonlinear equations. Thus, the numerical NR technique will be applied to solve these simultaneous equations to get the MLEs ($\hat{\theta}_{ML}$, $\hat{\beta}_{ML}$, $\hat{\lambda}_{ML})$ of ($\theta$, $\beta$, $\lambda)$.

Theorem 1. The MLEs of the parameters $\theta > 0$, $\beta > 0$ and $\lambda > 1$ are exist and unique.

Proof. To prove this theorem, three cases should be discussed, namely, (i) for given $\beta > 0$ and $\lambda > 1$, the MLE of $\theta$ exists and is unique; (ii) for given $\theta > 0$ and $\lambda > 1$, the MLE of $\beta$ exists and is unique; and (iii) for given $\theta > 0$ and $\beta > 0$, the MLE of $\lambda$ exists and is unique.

For the first case, when $\beta$ and $\lambda$ are known, the MLE of $\theta$ can be obtained from (11) as follows

where

Using (10) and $\hat{\theta}_{ML},$ we can write

where

Using inequality $\log\frac{\theta}{\hat{\theta}_{ML}}\leq\frac{\theta} {\hat{\theta}_{ML}}-1$, we have

Hence,

Equality holds if and only if $\theta = \hat{\theta}_{ML}$ or equivalent, $\ell(\theta, \beta, \lambda) = \ell(\hat{\theta}_{ML}, \beta, \lambda)$ if and only if $\theta = \hat{\theta}_{ML}.$ This proves that the function $\ell (\theta, \beta, \lambda)$ reaches its maximum value at the point $\ell (\hat{\theta}, \beta, \lambda)$ where$\ \beta$ and $\lambda$ are given. Using the similarity steps to prove (ii) and (iii).

For more results about the existence and uniqueness theorem of MLEs of the Lomax parameters, see Cramer and Schmiedt ^[38], and recently Qin and Gui ^[45]. The second derivatives are very important for distributions that are expressed using over one parameter for several reasons. For one of those reasons, they are going to confirm that maxima have been identified. In our case, the second partial derivatives of the LF can be expressed as

Arranging the second partial derivatives (14) in matrix form yields the Hessian matrix, which when multiplied by -1 in a statistical context, is known as a Fisher information matrix. In an optimization context, a sufficient condition for a stationary point to be a maximum is that the Hessian matrix is negative semi-definite. In the statistical literature, the inverse form of the Fisher information matrix is noted by the variance-covariance matrix. Thus, for the MLEs of the Lomax parameters, the Hessian matrix can be proposed as

The corresponding Fisher information matrix is

and the asymptotic variance-covariance matrix is

Thus, an approximate $(1-\alpha)100\%$ confidence intervals for $\theta$, $\beta,$ and$\ \lambda$ are obtained from the asymptotic normality of the ML results of intervals calculated according to

where $z_{\alpha/2}$ is defined as the percentile of the standard normal model with $\alpha/2$ right-tail probability. In the next two sections, the EM algorithm is applied to calculate MLEs of unknown parameters and their corresponding confidence intervals.

3.2.   Maximum likelihood estimation via EM algorithm

In this segment, our focus will be on the application of the EM algorithm in estimating the Lomax parameters, as it represents a general algorithm to obtain the MLEs, especially, in the problems of incomplete or missing data. The standard NR algorithm does not converge in some cases (see Pradhan and Kundu ^[46]), while the convergence is reliably an advantage of the EM method over other methods. Thus, the EM algorithm is utilized in this subsection as a good alternative to the NR approach, in order to calculate the MLEs of the unknown parameters of the Lomax model based on the fact that the available sample (constant partially PALT under progressive first-failure type-II censored samples) is considered incomplete data. This algorithm was proposed by Dempster et al. ^[47] and later was thoroughly discussed alongside its extensions in the book by McLachlan and Krishnan ^[48]. Now, let us symbolize the observed data by $X_{d} = (X_{d1}, X_{d2}, ..., X_{dm_{d}})$ whereas the censored data by $Z_{d} = (Z_{d1}, Z_{d2}, ..., Z_{dm_{d}})$, where each $Z_{di}$ is $1\times(k_{d}R_{di} +k_{d}-1)$ vector with $Z_{di} = (Z_{di1}, Z_{di2}, ..., Z_{di(k_{d}R_{di} +k_{d}-1)})$; $i = 1, 2, ..., m_{d}; \ d = 1, 2$ and they are not observable. Let $W = (X_{d}, Z_{d})$ represents the complete data set, hence the complete data LF will be in the following form

After ignoring the constants in the previous equation, the log-LF can be written as

The EM algorithm runs in two main stages: The expectation stage (E-stage), followed by the maximization stage (M-stage), which are repeated until achieving at least one convergence criterion. At each iteration, the missing data are filled with expected data and consequently, the estimations of the parameters are updated. About the requirements of the E-stage, one needs to compute the pseudo-log-LF, which can be obtained from $\ell_{c}\left(W; \theta, \beta, \lambda\right)$ by replacing any function of $z_{dij}$, say $g(z_{dij})$, with an accompanying conditional expectation "$E\left[g(z_{dij})|z_{dij} > X_{di}\right]$". Therefore,

For simplicity, $x_{di}$ has been used instead of $x_{di:m_{d}:n_{d}:k_{d}}$; $d = 1, 2$. Then, the M-stage begins, in which the pseudo-log-LF in (19) is maximized relative to $\theta$, $\beta, \ $and $\lambda$. Thus, at $s^{th}$ stage, the estimate of ($\theta, \beta, \lambda$) is ($\theta^{(s)}, \beta ^{(s)}, \lambda^{(s)}$), then ($\theta^{(s+1)}, \beta^{(s+1)}, \lambda^{(s+1)}$) can be obtained by maximizing $L_{s}(\theta, \beta, \lambda)$ with respect to ($\theta, \beta, \lambda$). Hence, by equating the partial derivative of $L_{s}(\Theta)$ to zero, a nonlinear system of equations can be generated as follows

and

The conditional distribution of $z_{di}$ follows the truncated Lomax distribution with left truncation at $x_{di}$ (see Ng et al. ^[49]). That is,

Hence,

Then, the conditional expectation $E\left(\ln(\beta+z_{di})|z_{di} > x_{di}\right)$ can be obtained as follows

using $u = \ln(\beta+t)$ and $dv = \left(\beta+t\right) ^{-\left(S_{d}\left(\lambda\right) \theta+1\right) }dt,$

and

3.3.   Fisher information matrix

Depending on the method presented by Louis ^[50], we can have an observed Fisher information matrix. Let $\bf{X}$ indicate observed data, $W$ denotes the complete data, $I_{X}(\Theta)$ shows the observed information, while $I_{W}(\Theta)$ and $I_{W|X}(\Theta)$ stand for the complete and missing information, then according to the Louis approach we have

where $I_{W}(\Theta)$ is obtained as

The complete data LF of the Lomax distribution in constant-stress PALT can be expressed as

Hence, the log-LF based on $X_{d}$ is given by

hence,

where

and

The Fisher information matrix for a single observation which censored at $i^{th}$ time failure is given by

where $f_{z_{di}}\left(z_{di}|x_{di}, \Theta\right)$ is given in (23) and the expected values of the second partial of the log-LF of $Z$ given $X$ are calculated as

consequently, the total missing information is obtained as

where $I_{W|X}^{i}(\Theta)|_{\Theta = (\theta, \beta, \lambda)}$ is given by (35). Finally, the asymptotic variance-covariance matrix of the MLE of $\Theta$ is then obtained as

Once the relevant variance estimates have been obtained from MLEs for $\theta$, $\beta, \ $and$\ \lambda$, asymptotic confidence intervals of $100(1-\gamma)\%$ can be easily generated using the normality property of MLEs.

4.   Bayesian estimation

In this section, the Bayes estimates for the unknown parameters $\theta$, $\beta$, and $\lambda$ are derived and discussed. Here, we consider independent gamma priors for $\theta$ and $\beta$, with PDFs given, respectively by

and

where ($a_{1}, b_{1}$) and ($a_{2}, b_{2}$) are hyper-parameters chosen to reflect prior knowledge about $\theta$ and $\beta$. The non-informative distribution with $a_{i} = 0$, $\ b_{i} = 0$; $i = 1, 2$ can be chosen if no prior knowledge is available. Moreover, a vague prior is selected for the acceleration factor $\lambda$ with the following PDF

The joint prior density function of $\theta$, $\beta$, and$\ \lambda$ is then obtained by

The joint posterior distribution of $\theta$, $\beta$, and\ $\lambda$ is reported using Bayes theorem by combining the LF of ($\theta, \beta, \lambda$) given in (8) and the joint prior distribution $\pi(\theta, \beta, \lambda)$ in (42) as

where $A$ is a normalizing constant and is given by

In Bayes' approach, to reach the best estimator, one should choose a loss function corresponding to each of the possible estimators. Here, the estimates are derived concerning two different types of loss functions, namely, symmetric loss functions and asymmetric loss functions. Squared error loss function (SELF) will be a representative of the first type, while both general entropy loss function (GELF) and LINEX loss function (LLF) will represent examples of choosing the second type. When an overestimation or underestimation occurs, the SELF is not appropriate. In this case, GELF and LLF can be taken as an alternative selection to estimate the parameters. The LLF is useful when overestimation is more serious than underestimation and vice-versa. Let $\psi$ be an estimator for the unknown parameter $\phi$, the loss functions under study can be presented as follows:

● Symmetric loss function, SELF: Form "$(\psi-\phi)^{2}$" and Bayes estimate "$E_{\phi}\left(\phi|\underline{t}\right)$".

● Asymmetric loss function:

- GELF: Form "$(\frac{\psi}{\phi})^{q}-q\ln(\frac{\psi}{\phi})-1$, $\ q\neq0$" and Bayes estimate "$\left(E_{\phi}\left(\phi^{-q} |\underline{t}\right) \right) ^{-\frac{1}{q}}$".

- LLF: Form "$\exp\left[c(\psi-\phi)\right] -c(\psi-\phi)-1$, $c\neq0$" and Bayes estimate "$-\frac{1}{c}\ln\left(E_{\phi}\left(\exp (-c\phi)|\underline{t}\right) \right)$".

It is not difficult to notice that, when $q = -1$, the Bayes estimation for any parameter $\phi$ based on GELF and SELF are the same (see Calabria and Pulcini ^[51]). Referring to (43), we can notice the difficulty of calculating integrals, and then the inability to obtain the joint posterior in a closed form that enables us to calculate Bayes estimations of the unknown parameters $\theta$, $\beta$, and $\lambda$. Therefore, in order to obtain these estimates, we will rely on the MCMC approach, which enables us to obtain simulated samples from the posterior distributions of the parameters. These generated samples will be used for calculating the point and interval estimation of unknown parameters. As for the mechanism of this method, it depends on the calculation of conditional posterior functions where the conditional distribution of $\theta$ given $\beta$ and\ $\lambda$ can be expressed as

where

Similarly, the conditional distribution of $\beta$ given $\theta$, $\lambda,$ and data can be reported as

Further, the conditional distribution of $\lambda$ given $\theta$, $\beta,$ and data, can be listed as

It can be seen that $\pi_{1}^{\ast}(\theta|\beta, \lambda, \underline {\text{x}} )$ and $\pi_{3}^{\ast}(\lambda|\theta, \beta, \underline {\text{x}} )$ are gamma densities. Thus, samples of $\theta$ and $\lambda$ can be generated using a gamma generator. Moreover, $\pi_{2}^{\ast}(\beta|\theta, \lambda, \underline {\text{x}} )$ cannot be reduced for drawing samples directly by standard methods. In such a case, to obtain Bayes' estimate for $\beta$, we can use one of the well-known algorithms in the MCMC method, which is Metropolis-Hastings (MH) algorithm model, which was presented in the literature by Metropolis et al. ^[52]. To apply this algorithm, we need to assume a proposal function to sample from it. In this algorithm, we may choose either symmetric or non-symmetric proposal distribution to decrease the rejection rate as much as possible. Since the marginal distribution of $\beta$ is not well-known, then the normal distribution is listed as a symmetric proposal distribution.

The Metropolis-Hastings steps are included in the Gibbs sampler to update $\beta$, while $\theta$ as well as $\lambda$ is updated directly from its full conditional. Below is a hybrid algorithm with Gibbs sampling steps for updating the parameters $\theta$ and $\lambda$ with MH steps for updating $\beta$:

Step 1: Choose an initial guess of ($\theta, \beta, \lambda$), say ($\theta^{(0)}, \beta^{(0)}, \lambda^{(0)}$), and set $l = 1.$

Step 2:Generate $\beta^{(l)}$ according to the following steps:

(i) Generate $\beta^{\ast}\ $from normal $N(\beta^{(l-1)}, Var(\hat{\beta}_{ML}))$ distribution where $Var(\hat{\beta}_{ML})$ denotes the variance of $\beta$;

(ii) Compute $r = min\left\{ 1, \dfrac{\pi_{2}^{\ast}(\beta^{\ast }|\theta^{(l-1)}, \lambda^{(l-1)}, \underline {\text{x}} )}{\pi_{2}^{\ast}(\beta ^{(l-1)}|\theta^{(l-1)}, \lambda^{(l-1)}, \underline {\text{x}} )}\right\};$

(iii) Generate a sample $u$ from the $U(0;1)$ distribution;

(iv) If $u\leq r,$ then set $\beta^{(l)} = \beta^{\ast}$; otherwise $\beta^{(l)} = \beta^{(l-1)}$.

Step 3:Generate $\theta^{(l)}$ from Gamma($m_{1}+m_{2}+a_{1}, m_{1}+m_{2}+a_{1}, b_{1}+T\left(x_{d}; S_{d}\left(\lambda^{(l-1)}\right), \beta\right)$).

Step 4:Generate $\lambda^{(l)}$ from Gamma($m_{2}, { \sum\limits_{i = 1}^{m_{2}}} \theta^{(l)}k_{2}\left(R_{2i}+1\right) \ln\left(\beta^{(l)} +x_{2i}\right) -\theta^{(l)}k_{2}n_{2}\ln\left(\beta^{(l)}\right)$)$.$

Step 5:Set $l = l+1.$

Step 6:Repeat steps ($2-5)$ $M$ times to obtain the desired number of samples.

We discard the initial $M_{0}$ number of burn-in samples and obtain estimates using the remaining $M-M_{0}$ samples. Hence, the Bayes estimate of $\phi = (\theta, \; \beta,$ or $\lambda)$ under SELF can be considered as the mean of the generated samples from the posterior densities as follows

Further, the approximate Bayes estimate for $\phi$, under LLF and GELF are then given, respectively by

and

To construct the highest posterior density (HPD) credible interval of $\phi = (\theta$, $\beta$, or $\lambda)$ using generated MCMC sampling procedure, we first refer to the ordered random sample generated in the previous algorithm in the form $\phi^{_{^{(1)}}} < \phi^{_{^{(2)}}} < ... < \phi^{_{^{(M)}}}$, then the $100(1-\gamma)\%$ confidence intervals for the parameter $\phi$ can be listed as

where $\left[.\right]$ refers to the greatest integer function.

5.   Simulation study and real data analysis

In this section, we introduce a numerical investigation of the estimation methods discussed in previous sections for the Lomax distribution using simulated data and a realistic data set.

5.1.   Censoring sample data

Here, the estimation procedure described in the previous sections is applied to the set of simulated progressive first-failure type-II censoring sample data under the constant-stress PALT. A data set of system lifetime is generated from the Lomax model with $\theta = 0.2$, $\beta = 2$ and $\lambda = 1.2$, respectively. Based on $N = 200$ ($n_{1} = n_{2} = 50$, $k_{1} = k_{2} = 2$), using the algorithm described in Balakrishnan and Sandhu ^[53], we simulate two samples of size $m_{1} = 20$ and $m_{2} = 30$ from the Lomax($\theta, \beta$) and Lomax($\theta\lambda, \beta$), respectively under the following two progressive censoring schemes

The simulated data are given in Table 1.

For the generated data, we first compute the MLEs of $\theta$, $\beta$, and $\lambda$ using the NR and EM methods, the results are listed in Table 2. Of the tabulated values, the estimates obtained using the NR method are far more distant from the true parameter value than those obtained under the EM algorithm. Then, next, we begin to calculate Bayes estimates, depending on the proposed loss functions and using informative priors ($a_{1} = 20$, $b_{1} = 100$, $a_{2} = 200$, $b_{2} = 100$). For LLF, $c = -10$ and $10$ were used as two options for the constant $c.$ These choices give more weight to underestimation and exaggeration, respectively. Similarly, two choices for $q$ such as $q = -10$ and $10$ are considered under GELF. The results are reported also in Table 2. From Table 2, it can be seen that, the estimates under the SELF, LINEX, and the GELF are close to the actual values of the parameters. From (15) and (34), the observed variance-covariance matrix based on NR and EM\ algorithms are given, respectively by

and

Table 2 also contains all interval estimations of the unknown parameters $\theta$, $\beta$, and $\lambda$ using the various methods mentioned. The Bayesian estimator gives narrower credible intervals as compared to the classical MLEs.

The numerical example shows that the means of the Bayes estimates and MLEs via NR or EM algorithm of the unknown parameters $\theta$, $\beta$ and $\lambda$ are close to the true values. It is seen that the Bayes estimates have better performances than MLEs where the Bayes estimates feature smaller bias than MLEs and that the interval lengths of HPD credible intervals are shorter than approximate confidence intervals, which also shows that HPD credible intervals are superior to approximate confidence intervals in terms of interval lengths.

5.2.   Engineering data analysis and interpretation

In this segment, we would like to know how the estimation techniques suggested in the above sections work in practice for the accelerated data set. The data are from the accelerated life test of oil breakdown times of insulating fluid under various levels of high voltage. Table 3 shows this data under two high test voltages stress levels (in kilovolt, Kv), see Nelson ^[12]. Here, the stress level (32 Kv) was assumed as normal conditions while the other stress level (36 Kv) was considered to be the accelerated conditions.

The initial shape is reported using the non-parametric approaches like histogram, kernel density, box, violin, TTT, and normal Q-Q plots in Figures 1 and 2. It is noted that the data is asymmetric as well as some outliers observations were founded. Moreover, the failure rate is decreasing.

Before further proceeding, we first check if the Lomax model can be used as a proper one to fit this data through the Kolmogorov-Smirnov (K-S) statistic. K-S distances and its associated p-values were calculated for the data set with the two stress levels, and the calculated values were: $0.1671$ ($0.7964$) and $0.1435$ ($0.9169$), respectively. Thus, we ensured that the Lomax distribution is considered an appropriate distribution for this data. The MLEs are listed in Table 4. By moving to Bayes estimates, since no prior information is available about the unknown population parameters, the noninformative (or vague) gamma priors are appropriate in this case. In this case the hyperparameters take zero values ($a_{i} = b_{i} = 0$, $i = 1, 2$). As noted earlier, Metropolis within the Gibbs algorithm was relied upon to generate $10000$ MCMC samples based on the MLEs of $\theta$, $\beta$, and $\lambda$ as initial values at the beginning of the algorithm. Trace plots for the first 10000 MCMC outputs of $\theta$, $\beta$, and $\lambda$ are shown in Figure 3. It is showing good convergence of MCMC procedure. Moreover, the histogram plots of generated samples of $\theta$, $\beta$, and $\lambda$ were presented in Figure 4. It is observed that the histograms of the generated samples quite well with the theoretical posterior density functions.

For Bayesian estimations under asymmetric LLF, it is known in the literature that, $c < 0$ implies that underestimation results in more penalty than overestimation and the reverse is true for $c > 0$. When $c$ close to zero LLF becomes symmetric, and it behaves roughly as the SELF. The Bayes estimates of $\theta$, $\beta$, and $\lambda$ are computed and reported in Table 4 with $c = -2$ and$\ 2$. The second asymmetric loss used here is GELF. In this loss function $q > 0$ means overestimation is more serious than underestimation and the opposite is true when $q < 0.$ The symmetrical case in which this loss function approximately corresponds to SELF occurs when $q = -1$. Also, the Bayes estimates are computed and recorded in Table 4 with $q = -2$ and$\ 2$. This table also displays the $95\%$ ML and Bayes confidence intervals with the corresponding lengths for all estimates based on the complete sample data.

Go ahead, based on the complete failure data in Table 3, in order to reduce cost and time, we assume that the number of items put on a life test is equal to $n_{d}\times k_{d}$, $d = 1, 2$, where $n_{d}$ denotes the number of groups and $k_{d}$ the number of items in each group. The previous data were randomly grouped into $n_{1} = n_{2} = 5$ sets, with $k_{1} = k_{2} = 3$ observations in each. Here, It is assumed that we observed only $m_{1} = m_{2} = 4$ data with arbitrarily chosen censoring schemes $\Re_{1} = \Re_{2} = (1, 0, 0, 0)$. Thus, the progressive first failure censored samples under the two stress are obtained as follows: $X_{1} = \{0.27, 2.75, 3.91, 82.85\}$ and $X_{2} = \{0.35, 1.69, 2.07, 2.71\}$. The MLEs relative to both NR, EM techniques and Bayesian MCMC method of unknown parameters $\theta$, $\beta$, and $\lambda$ are computed and listed in Table 5 for Nelson data under schema V: $n_{1} = n_{2} = 5$, $k_{1} = k_{2} = 3,$ $m_{1} = m_{2} = 3$. Moreover, the results of 95% approximate CI, credible intervals, and HPD intervals of $\theta$, $\beta$, and $\lambda$ are given in Table 5.

Finally, we can conclude that the estimated Lomax distribution provides an excellent good fit for the given data and the Bayes estimates fits the data better than MLEs. Also, we have decreased the number of units used here and yet we have obtained results close to the complete sample.

5.3.   The performance of different estimators

In this numerical section, we compare the performance of various estimators and confidence intervals mentioned above by conducting a simulated study. First, under normal and accelerated conditions from the Lomax($\theta$, $\beta$) and Lomax($\theta\lambda$, $\beta$) distributions, respectively, the algorithm which is given by Balakrishnan and Sandhu [53] can be utilized to generate progressive first-failure censored samples. This simulation was earned by considering different values of ($n_{1}, n_{2}, m_{1}, m_{2}, k_{1}, k_{2}$), and by choosing $(\theta, \beta, \lambda) = (0.4, 1.5, 1.3)$ in all cases. Further, we are giving results only for two different group sizes $(k_{1}, k_{2} ) = (1, 2)\ $and $(2, 1)$, $\ $and the censoring schemes used in this simulation are:

In order to improve the appearance of the numerical tables, the various censoring schemes have been represented by short notations as $(10, 0^{\ast4})$ denotes $(10, 0, 0, 0, 0)$. In each case, estimates of $\theta$, $\beta$, and $\lambda$ were calculated based on both the ML and Bayes methods. The NR and EM algorithm methods were used for ML estimation computation, while the MH algorithm with the Gibbs sampling algorithm is used to compute Bayes estimators. The informative prior distributions for both $\theta$ and $\beta$, with hyperparameters $a_{1}$, $b_{1}$, $a_{2}$, and $b_{2}$ were utilized and the associated hyperparameters are defined through the following discussion. According to the idea of Arabi Belaghi et al. ^[54], let $M_{1}$ random samples from Lomax($\theta = 0.4, \beta = 1.5$) distribution are available, and that ($\hat{\theta}_{i}, \hat{\beta}_{i}$); $i = 1, 2, \ldots, M_{1}$ are the MLEs of ($\theta$, $\beta$). Suppose that a parameter $\phi$ is a prior distributed as gamma with density proportional to $\phi^{h_{1}-1}e^{-h_{2}\phi}$, where $\phi$ can be either $\theta\ $or $\beta$. Then prior mean and variance are obtained by $\frac{h_{1}}{h_{2}}$and $\frac{h_{1}}{h_{2}^{2}}$, respectively. Now, equating the sample mean and associated variance of $\hat{\phi}_{i}$ with the mean and variance of the prior distribution, we get

From these equations, we get

Thus, if the unknown parameter is $\theta$, then the corresponding hyperparameters are estimated as $h_{1} = a_{1}, h_{2} = b_{1}$ and likewise, for $\beta$, we have $h_{1} = a_{2}$ and $h_{2} = b_{2}$. Following this, we compute hyperparameters by taking $M_{1} = 1000$ and for each $M_{1}$ a sample of size $50$ is taken into consideration to obtain desired estimates. Consequently, hyperparameters are assigned values as $a_{1} = 2.021$, $b_{1} = 4.9488$, $a_{2} = 5.827$ and $b_{2} = 3.2019$. Considering different loss functions, all the Bayes estimators are derived, and the corresponding approximate estimates are provided by applying Gibbs sampling under the MH sampling technique (see Section 4) approach. The different values of loss parameters are arbitrarily taken as $c = \pm1.5$ and $q = \pm1.5$, which produces estimators under both symmetric and asymmetric loss functions. The MH with the Gibbs algorithm is implemented using the normal distribution as the proposal distribution. Markov chains of size $5000$ are generated and the first $1000$ of the observations are removed to eliminate the effect of the starting distribution. While utilizing this algorithm, we take the corresponding MLEs as an initial guess for the respective unknown parameters and the associated variance-covariance matrix as the variance-covariance matrix of the normal proposal density. All the simulations were performed using the Mathematica 9 software and the results were based on $N_{S} = 1000$ Monte Carlo runs. We compute average values (Avg) and corresponding mean square error (MSE) of different estimates. Here,

where $\hat{\phi}$ is the estimate of $\phi$. Also, we obtain the average length (AL) of $95\%$ confidence/HPD credible intervals and $95\%$ CP of the parameters based on the simulation. The results of the Monte Carlo simulation study are presented in Tables A1–A4 (see appendix). From these tables the following conclusions are made:

$(1)$ The MSEs of all estimates decrease as the sample sizes increase in all cases, as expected.

$(2)$ The MLEs of $\theta$, $\beta$, and $\lambda$ using the EM algorithm have smaller MSEs than the MLEs using the NR algorithm. Hence, the MLEs via the EM algorithm perform better than those obtained by the NR method.

$(3)$ In most simulations, the Bayes estimates outperform the MLEs for the estimation of $\theta$, $\beta$, and $\lambda$. So, in general, we would recommend using the Bayes estimate of the unknown parameters of Lomax distribution based on a first-failure progressively censoring scheme under constant-stress PALT. Also, from the tabulated values, we notice that the Bayes estimates of $\theta$, $\beta$, and $\lambda$ under asymmetric loss function (LLF, GELF) are sensitive to the value of the parameters $c$ and $q$. It is observed that, as expected, the positive values of $c$ and $q$ respectively under LINEX and GELF lead to smaller estimates as compared to the negative value of $c$ and $q$. Moreover, the choice $c = 1.5$ for the LINEX loss seems to be reasonable. For the GELF, $q = 1.5$ is a better value in computing the Bayes estimate based on the LINEX or the SELE. Hence, it can be concluded that the asymmetric loss function makes the Bayes estimates attractive for use, in reality, the scale parameters $c$ and $q$ of the LINEX and GELF make one estimate the unknown parameters with more flexibility.

$(4)$ Finally, in most cases, the MSEs and the AL of the intervals for ($\theta, \beta, \lambda)$ are maximized when the censoring is at the end. Hence, when we compare progressive type-II censoring and type-II censoring plans, the progressive type-II plan provides a better result.

6.   Conclusions

In this work, different parameter estimation methods for the two parameters of the Lomax distribution are discussed based on constant stress partially accelerated life test using the progressive first failure type-II censored sample data. We derived the MLEs and Bayes estimates of the parameters as well as the acceleration factor and the corresponding confidence intervals. We have also exploited the EM algorithm to obtain the MLEs for the unknown parameters. In addition, the supplemented EM algorithm was used for finding the asymptotic covariance matrix of the ML estimators. Meanwhile, because the explicit expressions of MLEs of part of parameters cannot be derived, we look to the Bayesian method for support. Based on the square error, LINEX, and general entropy loss functions, the Bayes estimates are obtained on the premise of independent gamma priors. The posterior distributions of unknown parameters indicated that some parameters do not follow a well-known distribution. Consequently, we utilized MH sampling within the Gibbs sampling steps algorithm to calculate the Bayes estimates with associated HPD credible intervals. In order to compare the performance of all mentioned methods directly and to investigate the effect of different values of $n$, $m$, and $k$, the simulation study was implemented. The average, MSEs, the average widths of the interval estimators, and the coverage percentages of each estimator are illustrated. The simulation results show that for fixed values of $m$ and $k$, the performances of both estimation methods improve with $n$, and for fixed $n$ and $k$, their performances improve with $m$, also for fixed $n$ and $m$, their performances improve with $k$.

Further, the numerical results demonstrated that the Bayes method based on informative priors outperformed the ML method under both NR and EM methods. Also, it is perceived that the estimates obtained under the EM method generally perform well compared to those obtained under the NR method in terms of yielding relatively low values of MSEs and average widths of the interval estimators. Further, the Bayes estimates based on LINEX, and general entropy loss functions are more efficient than the squared error loss function under informative and non-informative priors. Practically, the study has demonstrated that the Lomax model has provided good flexibility for modeling oil breakdown times of insulating fluid data sets. The novelty in this research is that if a progressively first failure type-II censored sample is used, different sample sizes ($k_{1}\neq k_{2}$) can be considered in each group, and this is fully consistent with real examples when doing life tests. Finally, in this work, although we have mainly considered progressively first failure type-II censoring and Lomax distribution, the same method can be extended for other distribution and censoring schemes also. There are numerous further works to be done here, for example, the design of optimal censoring schemes, inference of competing risks model with more failure factors, and making statistical prediction of the further order statistics based on the PALTs from Lomax distribution, these topics can be investigated in the future. Finally, we recommend the use of MCMC, and EM procedures based on partially accelerated life test using the progressive first failure type-II censored on data related to life testing, reliability modeling as well as biological analysis.

Acknowledgements

The authors extend their appreciation to the Deputyship for Research & Innovation, Ministry of Education, Saudi Arabia for funding this research work through the project number (QU-IF-2-5-3-25108). The authors also thank to Qassim University for technical support.

Conflict of interest

The authors declare no conflicts of interest.

Appendix