Numerical analysis of the fractional nonlinear waves of fifth-order KdV and Kawahara equations under Caputo operator

1.   Introduction

Because of advancements in manufacturing and technical development, products and devices are becoming increasingly reliable, making typical life tests under normal conditions difficult, if not impossible. For the industrial markets, such tests are too time-consuming and costly to get sufficient information about a product's lifetime distribution or even a prototype. As a result, the accelerated life test (ALT) is becoming more popular as well as important, as it gives information on a highly reliable product's lifetime in a short amount of time, see ^[1,2]. More failures can be collected quickly by conducting the life test at higher stress levels than normal operating conditions. A suitable stress-response regression model is then used to estimate the lifetime distribution at the usage stress. As a particular class of ALT, the progressive-stress test implements a special stress-loading scheme where the stress increases in time, see ^[3,4,5,6].

Censoring often occurs when some lifetimes of products are missing or when some experimental design purposes are being implemented. Type-I and type-II censoring are the most prevalent schemes. These two types of censoring do not have the flexibility to allow the experimenter to withdraw units from a life test at different stages during the experiment. Because of this lack of flexibility, progressive censoring is proposed as a more general censoring technique. This technique enables the experimenter to withdraw units from a life test at different stages, other than the final point, through the test experiment.

One of various risk factors could cause the units to fail in reliability, engineering, biomedicine, physical studies, and other fields. If the risks are unclear, an issue may occur regarding the factor that caused the unit to fail, and as a result, the lifetime associated with this particular risk cannot be determined. In such cases, the maximum (minimum) lifetime value with all risks can only be observed. This is referred to as a "series (parallel) system" in the literature. These two systems can be joined to form a new system known as a "series-parallel system", see its description in Figure 1. In a series system, components are connected in such a way that the failure of a single component leads to the failure of the system. On the other hand, a system in which failure of all components leads to the failure of the system is called a parallel system. Finally, a series-parallel system is a system in which $m$ subsystems are connected in series and each subsystem consists of $k$ components connected in parallel, see Figure 1. For such a system, failure is observed if a subsystem fails. The technique of compounding of distribution functions may be used to construct the lifetime distribution of the series-parallel system.

Abdel-Hamid and Hashem ^[7] introduced a new distribution for the series-parallel system by compounding two Poisson distributions (truncated at zero) with an exponential distribution. They used six estimation methods to estimate the included parameters. In ^[8], they obtained a new distribution by compounding two discrete distributions with a mixture of continuous distributions based on a parallel-series system. Nadarajah et al. ^[9] applied the progressive-stress ALT technique using type-I progressively hybrid censored data with binomial removal to components connected as a parallel-series structure. Hu et al. ^[10], suggested an ideally distributed series-parallel system, as well as two analytical reliability assessment approaches to analyze the reliability of a distributed power system. Hashem and Alyami ^[11] introduced a new distribution based on a parallel-series system.

In this article, we introduce a new distribution that describes the lifetime of a series-parallel system when the lifetimes of the included components are subject to a mixture of exponential and gamma distributions. The distribution will be called Poisson-geometric-Lomax distribution (PGLD) which can be obtained by compounding truncated Poisson, geometric, and Lomax distributions

The following is how the rest of the article is put together: The distribution of a series-parallel system based on a finite mixture of distributions is constructed in Section 2. In Section 3, the model is described based on progressive-stress ALT and progressive type-II censoring. In Section 4, estimation of the parameters is obtained using the maximum likelihood and Bayes methods. In Section 5, the parameters are estimated numerically using maximum likelihood and Bayes estimation methods. An illustrative example, based on two real data sets, is studied in Section 6. A simulation study is worked done in Section 7. Concluding remarks followed by certain features and motivations of the PGLD as well as future work are presented in Section 8.

2.   Finite mixture of distributions and distribution of a series-parallel system

A mixture of distributions is considered a combination of statistical distributions, that occur when sampling from inhomogeneous populations (or mixed populations) with a different probability density function (PDF) in each component.

Let $F(y \mid \omega)$ be the cumulative distribution function (CDF) of a continuous random variable (RV) $Y$, with realization $y$, that is dependent on a continuous RV $\Omega$, with realization $\omega$, and $\Pi(\omega)$ be the CDF of $\Omega$. The marginal CDF $F(y)$ is given by

which is called a mixture (according to Teicher ^[12]) of the CDFs $F(y \mid \omega)$ and $\Pi(\omega)$. Fisher ^[13] called $F(y)$ "compound distribution".

The corresponding PDF is given by

If the RV $\Omega$ assumes only a finite number of points $\{ {\omega_j}, \, j = 1, \dots, \kappa\}$, $\Pi(\omega)$ is then a mass function and assigns positive probabilities to only ${\omega_j}$. The integral in (2.1) is then replaced by a sum to give a { finite mixture} of the form

Suppose, in (2.3), $\Pi({\omega_j}) = p_j$, $j = 1, \dots, \kappa$, and $F(y; {\omega_j}) = F_j(y)$. Then (2.3) takes the form

A corresponding finite mixture of PDFs is given by

In (2.4) and (2.5) the masses $p_j$ are called mixing proportions with the conditions:

The functions $F_j$ and $f_j$ are called $j^{th}$ components in the finite mixture of CDFs (2.4) and PDFs (2.5), respectively. For more details on finite mixture of distributions, see ^[14,15]. Several distributions may arise due to compound distribution (2.2) and mixture distribution (2.5), see ^[16], such as

(i) If in (2.2), $f(y\mid \omega)$ has $\text{exp}(y\mid\omega)$ distribution and $\pi(\omega)$ has $\text{gamma}(\alpha, \gamma)$ distribution, then we get Lomax$(\alpha, \gamma)$ distribution (LD) as follows:

Therefore, the corresponding CDF is given by

(ii) If in (2.5), $\kappa = 2, p_1 = \frac{\theta}{\theta+1}$, $f_1(y)$ has $\text{exp}(\theta)$ distribution and $f_{2}(y)$ has $\text{gamma} (2, \theta)$ distribution, then we get Lindley distribution as follows:

2.1.   Mixture of two components with mixing proportion as a random variable

If in (2.2), $f(y\mid \Omega = \omega) = \omega f_1(y)+(1-\omega) f_2(y)$ and the RV $\Omega$ subjects to a beta distribution, $B(a, b)$, then

2.2.   Distribution of a series-parallel system

Several new distributions of the series-parallel system may emerge by choosing different distributions for the number of subsystems as well as those describing the lifetimes of the included components, as explained in ^[7,8], but in exchange, several questions may arise, for example:

● Does the presented distribution have any interesting properties?

● Does it arise naturally in some observable process like natural phenomena?

● Does it have any motivations?

There is no need for new distributions if there are no positive answers to the above questions. A response to these questions will be illustrated in Section 8.

According to ^[7,8], the following theorem gives the CDF and PDF of the PGLD which represents the failure time distribution of the series-parallel system when the number of series subsystems and the number of their components, that are connected as a parallel structure, are RVs subject to truncated Poisson and geometric distributions, respectively. At the same time, the failure times of the components have a finite mixture of distribution functions each of which may be responsible for a different cause of failure.

Theorem 2.1. Suppose that, for $i = 1, 2, \dots, k_j$, $j = 1, 2, \dots, m$, a mixed system of series-parallel structure type has the lifetime $X = \min_{ j }(\max _{i}Y_{ij})$, where $Y_{ij}$ are IID RVs, see Figure 1, with PDF(CDF) $f_{Y}(y)$($F_{Y}(y)$) of Lomax distribution given by $(2.6)$ ($(2.7)$). Consider $M$ and $K_j$ are two discrete RVs subject to truncated Poisson and geometric distributions with PMFs $P(M = m) = \frac{e^{-\lambda} \lambda^{m}}{m! (1-e^{-\lambda})}, \, m = 1, 2, \dots, \, (\lambda > 0)$ and $P(K_j = k_j) = (1-\theta)\, \theta^{k_j-1}, \, k_j = 1, 2, \dots, (0 < \theta < 1)$, respectively. Then the distribution of $X$ has the PGLD with CDF and PDF given, respectively, by

where

Proof. The proof is similarly as in ^[7].

Remark 2.1. It is worth noting that while $\gamma$ is a scale parameter for LD with CDF $(2.7)$, it is also a scale parameter for PGLD with CDF $(2.9)$.

The PDF and hazard rate function (HRF), $h(x) = f(x)/(1-F(x))$, of the PGLD are plotted in Figure 2, in which one can observe that the PDF exhibits decreasing and unimodal shapes while the HRF exhibits a unimodal shape and may be sudden fluctuation at its end. These fluctuations usually imply that the product's performance has degraded over time. Non-stationary data can exhibit this characteristic, and the PGLD can help to represent such data. The non-stationary nature of failure times may help the researcher forecast how some items will behave in the environment.

The features and motivations of the PGLD with CDF and PDF (2.9) and (2.10) are summarized in the last section.

In the following section, some important properties of the PGLD are given.

3.   Properties of the PGLD

In this section, some important properties of the PGLD, such as the $q$-th quantile, mode, $r$-th moment, mean residual lifetime, Bonferroni and Lorenz curves, Rényi and Shannon's entropies, PDF and CDF of the $i$-th order statistic, are given.

Theorem 3.1. The q-th quantile $x_{q}$ of the PGLD with CDF $(2.9)$ can be obtained as

Proof. By solving the equation $F_X(x_q) = q$ with respect to $x_q$, the proof can be achieved immediately.

Remark 3.1. As a particular case, the median of PGLD with CDF $(2.9)$ can be obtained by putting $q = 1/2$ in Eq $(3.1)$.

Theorem 3.2. Let $X$ be a RV subject to the PGLD with PDF $(2.10)$. Then the mode is given by

where

Proof. The mode can be directly obtained by solving $\frac{{\text d}\ln[f_{X}(x)]}{{\text d}x} = 0$ with respect to $x$.

Theorem 3.3. Let $X$ be a RV subject to PGLD with PDF $(2.10)$. Then, for $r = 1, 2, \dots$, the $r$-th moment of $X$ is given by

where $y_{\imath}, \, \nu_{\imath}$ are the zeros and the corresponding Christoffel numbers of the Legendre-Gauss quadrature formula on the interval (-1, 1), see Canuto et al. ^[17].

Proof. The $r$-th moment of $X$ is given by

where $y^{*} = \frac{1+y}{1-y}$, and

The last integral can be approximated, by using Legendre-Gauss quadrature formula, as

where

and $L_{N}(.)$ is the Legendre polynomial of degree $N$.

Theorem 3.4. The mean residual lifetime of the PGLD is given by

where $\nu_{\imath}$ is given by $(3.5)$.

Proof. The mean residual lifetime of the PGLD is given by

where $S(x) = 1-F(x)$ is the survival function.

The last integral can be approximated, by using Legendre-Gauss quadrature formula, as

where $\nu_{\imath}$ is given by (3.5).

Theorem 3.5. Let $X$ be a RV subject to the PGLD with PDF $(2.10)$. Then, the Bonferroni curve (${BC}$) and Lorenz curve (${LC}$) are given, respectively, by

where $\nu_{\imath}$ is given by $(3.5)$, $0 < \eta < 1$ and

and $\Phi^{(1)}(.)$ and $\mathfrak{m}^{(1)}$ are given, respectively, by $(3.3)$ and $(3.4)$ at $r = 1.$

Proof. The Bonferroni curve of PGLD is given by

The Lorenz curve of PGLD is given by

The Bonferroni and Lorenz curves are plotted in Figure 3.

Theorem 3.6. Let $X$ be a RV subject to the PGLD with PDF $(2.10)$. Then, the Rényi and Shannon's entropies of $X$ are given, respectively, by

where $\nu_{\imath}$ is given by $(3.5)$ and

Proof. The Rényi entropy of $X$ is given by

where $\ell > 0$ and $\ell\neq 1$.

Based on PDF (2.10), we obtain

Then,

The Shannon's entropy of $T$ is given by

Theorem 3.7. Let $X_{1}, \dots, X_{n}$ is a random sample from the PGLD with CDF $(2.9)$ and PDF $(2.10)$. Then, the PDF and CDF of the $i$-th order statistic, say $X_{i:n}$, are given, respectively, by

Proof. The PDF $f_{j:n}(t)$ of the $j$-th order statistic, see ^[18,19], is given by

where $F_{X}(x)$ and $f_{X}(x)$ are given by Eqs (2.9) and (2.10), respectively.

Therefore,

The CDF $F_{j:n}(x)$, corresponding to PDF (3.10), is given by

4.   Model description

Based on units connected in a series-parallel structure, we discuss, in this section, the application of a progressive-stress model to units whose lifetime distribution subjects to the PGLD with CDF (2.9). We assume that the units are subject to progressive type-II censoring and that the number of surviving units eliminated follows a binomial distribution.

It is well-known in the previous literature regarding the progressive-stress model that the stress is considered an increasing linear function of time. Now, in the current article, the stress is supposed to be represented by an increasing nonlinear function of time.

4.1.   Basic assumptions

$(1)$ The lifetime of units under normal conditions is governed by CDF (2.9) of PGLD with scale parameter $\gamma$.

$(2)$ According to the progressive-stress model, it is assumed that the stress $s$ is a function of time $t$ and affects the scale parameter $\gamma$ of CDF (2.7) which is also a scale parameter of CDF (2.9).

$(3)$ The parameter $\gamma$ follows the inverse power law with two parameters $c$ and $d$. This means that $\gamma \equiv \gamma(t) = \frac{1}{c(s(t))^{d}}$. For the sake of simplicity, we'll assume that the parameter $d$ assigns the value 1 from now on.

$(4)$ The progressive-stress $s(t)$ assigns an increasing nonlinear function of time, $s(t) = \text{sinh}(v\, t)$, which is also continuous and differentiable for $t > 0$.

$(5)$ To start the testing procedure, the total $N$ units to be tested are split into $\hbar(\geq 2)$ groups, each of them consists of $n_i$ units under progressive-stress $s_i(t) = \text{sinh}(v_i t), \; i = 1, 2, \dots, \hbar,$ such that the stress rates satisfy $0 < v_1 < v_2 < \dots < v_\hbar.$

$(6)$ At whatever stress level $v_i, \; i = 1, 2, \dots, \hbar$, a unit's failure mechanisms are the same.

$(7)$ The cumulative exposure model holds for the effect of changing stress, see ^[1].

According to Steps 1, 3, and 5, and using Assumption 7, we can write

where $c$ is a parameter that should be estimated.

Therefore, if $G_{iY}(.)$ denotes the CDF for a unit in group $i$ under progressive-stress model, then using Assumption 7, it takes the form

where $F_{iY}(.)$ is the CDF (2.7) of LD, included in (2.9), under group $i$, with scale parameter value equal to 1 and

where $\Phi_i(x)$ is given by (4.1).

The corresponding PDF of (4.2) is given by

Then, by replacing $F_Y(y)$ and $f_Y(y)$ with CDF (4.2) and PDF (4.4), respectively, CDF (2.9) under progressive-stress ALT becomes

The corresponding PDF of (4.5) is given by

where

and $\Psi_i(x)$ is given by (4.3).

4.2.   The model under progressive type-II censoring

There are a variety of censored tests available. Progressive type-II censoring is one of the most widely used censored tests. It is implemented under progressive-stress model as follows: Suppose that, in the $i$-th group, $i = 1, 2, \dots, \hbar$, $n_i$ units are put through a life test and the experimenter specifies beforehand quantity $w_i$, the number of failure units that will be observed. Now, when the first failure occurs, $r_{i1}$ of the remaining $n_i-1$ surviving units are eliminated from the experiment at random. Continuing, when the second failure occurs, $r_{i2}$ of the remaining ${n_i}-r_{i1}-2$ surviving units are eliminated from the experiment at random. Finally, at the time of occurring the $w_i$-th failure, all the remaining $r_{iw_i} = n_i-w_i-r_{i1}-r_{i2}-\dots-r_{i\, w_i-1}$ surviving units are eliminated from the experiment at random. Note that, in this scheme, $r_{i1}, r_{i2}, \dots, r_{iw_i}$ are all prefixed. However, these numbers may occur at random, in some practical situations. For example, in some reliability experiments, an experimenter may assess that it is unsuitable or too risky to perform the test on some of the tested units even though these units have not failed. In such situations, the pattern of elimination at each failure is random. Such a situation leads to progressive censoring with random removals, see ^[20,21].

5.   Estimations under progressive-stress ALT and progressive type-II censoring with binomial removals

We discuss, in this section, two estimation methods (maximum likelihood estimation (MLE) and Bayesian estimation (BE)) for the parameters included in CDF (4.5) under progressive type-II censoring with binomial removals.

5.1.   Maximum likelihood estimation

For group $i$, $i = 1, 2, \dots, \hbar$, suppose $x_{i1}, x_{i2}, \dots, x_{in_i}$ are $n_i$ lifetimes from a population with CDF (4.5) and PDF (4.6), and ${\bf x} = \left(x_{i1:w_i:n_i}^{(r_{i1}\dots, r_{iw_i})}, x_{i2:w_i:n_i}^{(r_{i1}\dots, r_{iw_i})}\right.$, $\left.\dots, x_{iw_i:w_i:n_i}^{(r_{i1}\dots, r_{iw_i})}\right)$ is a vector of $w_i$ progressively type-II ordered lifetimes out of $n_i$ with progressive censoring scheme $(r_{i1}\dots, r_{iw_i})$. Suppose also $w_i$ and $r_{ij}$ are all predetermined before the test. From now on, we will write $x_{ij}$ instead of $x_{ij:w_i:n_i}^{(r_{i1}\dots, r_{iw_i})}$, $j = 1, 2, \dots, w_i$, for simplicity. Then, based on Eqs (4.2)–(4.7), the conditional likelihood function is given by

The numbers $r_{ij}$ may occur at random in some practical scenarios as a result of unanticipated dropout experimental units. Therefore, we assume that $R_{ij} \; (i = 1, 2, \dots, \hbar, \; j = 1, 2, \dots, w_i-1)$ are RVs, with realizations $r_{ij}$, subject to the following binomial distributions

whereas,

Then

where, for fixed $i$ and $j$, $r_{ij}$ is an integer number satisfying $0\leq r_{ij} \leq n_i-w_i-(r_{i1}+r_{i2}+\dots, r_{i(j-1)}), \; i = 1, 2, \dots, \hbar, \, j = 2, 3, \dots, w_i-1,$ and

Furthermore, suppose that $R_{ij}$ is independent of $X_{ij}$ for $i = 1, \dots, \hbar$, $j = 1, \dots, w_i$. Then

Since $P({\bf R} = {\bf r})$ does not depend on the parameters $(\alpha, \lambda, \theta, c)$, hence the MLE of them can be derived by maximizing (5.1). Similarly, (5.1) does not involve the binomial parameter $p$, the MLE of $p$ can be found by maximizing (5.3) directly. Thus

Therefore,

The local Fisher information matrix, ${\bf{I}}$, for $(\hat{\alpha}, \hat{\lambda}, \hat{\theta}, \hat{c})$ is the $4\times 4$ symmetric matrix of negative second partial derivatives of $£ = \text{ln} L_1$ with respect to $\alpha$, $\lambda$, $\theta$, and $c$, see ^[1]. Let $(\varphi_1 = \alpha, \, \varphi_2 = \lambda, \, \varphi_3 = \theta, \, \varphi_4 = c)$. Therefore, matrix ${\bf{I}}$ is given by

where the caret $\hat{}$ denotes that the derivative is computed at $(\hat{\varphi_1} = \hat{\alpha}, \, \hat{\varphi_2} = \hat{\lambda}, \, \hat{\varphi_3} = \hat{\theta}, \, \hat{\varphi_4} = \hat{c})$. It is easy to get the matrix's elements.

The local estimate ${\bf{V}}$ of the asymptotic variance-covariance matrix of $(\hat{\alpha}, \hat{\lambda}, \hat{\theta}, \hat{c})$ can be obtained by inverting matrix ${\bf{I}}$. Therefore,

The sampling distribution of $\frac{\hat{\varphi_i}-\varphi_i}{\sqrt{\text{var}(\hat{\varphi_i})}}, \; i = 1, \dots, 4,$ follows the general asymptotic theory of MLEs and hence it may be approximated by a standard normal distribution which can be used to create confidence intervals (CIs) for unknown parameters.

Therefore, a two-sided $(1-\eta^{\star})100\%$ asymptotic CIs for the parameters $\varphi_i, \; i = 1, \dots, 4$, can then be created as follows:

where $z_{\eta^{\star}/2}$ is the value of a standard normal RV that leaves an area $\eta^{\star}/2$ to the right and $\sqrt{\text{var}(\hat{\varphi_i})}$ can be determined from (5.5).

It's commonly beneficial in practical applications to have a concept of how long a test should last. Because the time it takes to complete the test is proportional to the expense. For $i = 1, \dots, \hbar$, the expected termination time under progressive type-II censoring scheme with binomial removals is given by, see ^[22].

where ${{\bf R}_i} = (R_{i1} = r_{i1}, \dots, R_{i(w_i-1)} = r_{i(w_i-1)})$, ($F_X(x), f_X(x)$) is given by ((4.5), (4.6)), $\mathcal{A}_i = \sum_{k = 1}^{w_i}\ell_k$, $h(\ell_k) = \ell_1+\dots+ \ell_k+k$, and $C({{\bf r}_i}) = n_i(n_i-r_{i1}-1)(n_i-r_{i1}-r_{i2}-2)\dots(n_i-\sum_{k = 1}^{w_i-1}(r_{ik}+1))$. Therefore,

where $q(r_{i1}) = n_i-w_i, \, q(r_{ij}) = n_i-w_i-r_{i1}-r_{i2}-\dots-r_{i(j-1)}, \, i = 1, 2, \dots, \hbar, \, j = 2, 3, \dots, w_i-1,$ and $P({\bf R} = {\bf r})$ is given by Eq (5.3). Then compute the ratio of expected experiment time (REET) as follows:

5.2.   Bayesian estimation

Based on two asymmetric loss functions (general entropy (GE) and linear exponential (LINEX)), we discuss the Bayesian estimation of the parameters of CDF (4.5). Because they offer overestimation or underestimation of the parameters, symmetric loss functions may be unsuitable in many real-life situations. Overestimating the parameters can have worse or worse repercussions than underestimating them or vice versa. As a result, asymmetric loss functions have been the subject of research, see ^[23,24].

The LINEX loss function was suggested by Varian ^[25]. It is given by

where $\xi = \tilde{\Theta}-\Theta$ and $\tilde{\Theta}$ is the LINEX estimator of $\Theta$.

The Bayes estimate under LINEX (BEL) loss function of $\Theta$ is given by

The GE loss function was introduced by Calabria and Pulcini ^[26]. It is given by

The Bayes estimate under GE (BEG) loss function of $\Theta$ is given by

5.2.1.   Prior and posterior functions

Assume that the experimenter's prior belief is measured by a function $\vartheta(\alpha, \lambda, \theta, c)$, with all parameters being independent and having log-normal distributions except $\theta$, which has a beta distribution. Therefore, if $(\varphi_1 = \alpha, \, \varphi_2 = \lambda, \varphi_3 = c, \varphi_4 = \theta)$, then the prior function of $\varphi_i$, $i = 1, 2, 3$, is given by

and the prior function of $\varphi_4$ is given by

The joint prior density function is then calculated as follows:

From (5.1) and (5.11), the joint posterior density function is then calculated as follows:

where

It can be noted that $K$ involves quad integral and it is not reducible in a closed form, and hence generating samples directly from the joint posterior density is not possible. The MCMC algorithm, presented in the following subsection, can be implemented in this case in which we need the following conditional posterior distributions of the parameters $\alpha$, $\lambda$, $\theta$ and $c$,

5.2.2.   MCMC algorithm

The Metropolis-Hasting technique can be used if the conditional posterior distribution isn't one of the known parametric distributions. Then, it can be used to generate samples from conditional posterior distributions using the MCMC algorithm, see ^[27].

The following procedure can be used to compute BELs and BEGs of $\alpha$, $\lambda$, $\theta$, and $c$,

$(1)$ Assign some initial values of $\alpha$, $\lambda$, $\theta$ and $c$ say $\alpha_{0}$, $\lambda_{0}$, $\theta_{0}$ and $c_{0}$.

$(2)$ For $i = 1$, using Metropolis-Hastings technique, generate $\alpha_{i}$, $\lambda_{i}$, $\theta_{i}$ and $c_{i}$ from the conditional posterior distributions presented in (5.13).

$(3)$ Repeat Step 2, $\mathbb{M}$ times.

$(4)$ Calculate the BELs of $\alpha$, $\lambda$, $\theta$ and $c$, using Eq (5.7), as

where $\mathbb{W}$ is the burn-in period.

$(5)$ Calculate the BEGs of $\alpha$, $\lambda$, $\theta$ and $c$, using Eq (5.8), as follows:

$(6)$ The highest posterior density (HPD) credible interval of parameter $\alpha$ can be computed by arranging in an ascending order $\alpha_{\mathbb{W}+1}$, \dots, $\alpha_{\mathbb{M}}$ as $\alpha_{[1]} < \dots < \alpha_{[\mathbb{M}-\mathbb{W}]}$, then compute the lower and upper bounds of all $(1-\eta^{\star})100\%$ credible intervals of $\alpha$ as follows:

where $[x]$ denotes the largest integer number less than or equal to $x$. Then the HPD credible interval of $\alpha$ is the one with the shortest length. The HPD credible interval of $\lambda$, $\theta$ and $c$ can also be computed in the same way.

$(7)$ The symmetric credible intervals of $\alpha$ can be computed as follows:

The symmetric credible intervals of $\lambda$, $\theta$ and $c$ can also be computed in the same way.

6.   Illustrative example based on real data

Two real data sets are proposed, in this section, for fitting and comparing the PGLD, Poisson-Lomax distribution (PLD) ^[28], geometric-Lomax distribution (GLD), exponentiated LD (ELD) ^[29] and LD with CDF (2.7).

The CDFs of PLD, GLD, and ELD are given, respectively, by

It can be noticed that the GLD is a special case of the PGLD as $\lambda\rightarrow 0^+$.

● The first data set is taken from ^[30]. It consists of 72 exceedances for the years from 1958 to 1984, rounded to one decimal place. The data are given as follows:

1.7, 2.2, 14.4, 1.1, 0.4, 20.6, 5.3, 0.7, 1.9, 13.0, 12.0, 9.3, 1.4, 18.7, 8.5, 25.5, 11.6, 14.1, 22.1, 1.1, 2.5, 14.4, 1.7, 37.6, 0.6, 2.2, 39.0, 0.3, 15.0, 11.0, 7.3, 22.9, 1.7, 0.1, 1.1, 0.6, 9.0, 1.7, 7.0, 20.1, 0.4, 2.8, 14.1, 9.9, 10.4, 10.7, 30.0, 3.6, 5.6, 30.8, 13.3, 4.2, 25.5, 3.4, 11.9, 21.5, 27.6, 36.4, 2.7, 64.0, 1.5, 2.5, 27.4, 1.0, 27.1, 20.2, 16.8, 5.3, 9.7, 27.5, 2.5, 27.0.

● The second data set represents the marks of 48 slow space students in Mathematics in the final examination of the Indian Institute of Technology, Kanpur in year 2003, ^[31]. The data are given as follows:

29, 25, 50, 15, 13, 27, 15, 18, 7, 7, 8, 19, 12, 18, 5, 21, 15, 86, 21, 15, 14, 39, 15, 14, 70, 44, 6, 23, 58, 19, 50, 23, 11, 6, 34, 18, 28, 34, 12, 37, 4, 60, 20, 23, 40, 65, 19, 31.

The Kolmogorov-Smirnov (K-S) statistic and its corresponding p-value are used to check the validity of the PGLD, PLD, GLD, ELD, and LD to fit the above two data sets. The results are shown in Table 1. Also, the results of a comparison among these five distributions using some criteria such as the Akaike information criterion (AIC), consistent AIC (CAIC), and Bayesian information criterion (BIC), are shown in Table 1, where

where $£(\hat{\beta})$ stands for the log-likelihood function calculated at the MLE $\hat{\beta}$ of $\beta$, $r$ and $w$ denote the number of parameters and the sample size, respectively.

According to the values of K-S statistic and its corresponding p-value, presented in Table 1, it's worth noting that the PGLD has the smallest(largest) K-S (p-)values than those for the PLD, GLD, ELD, and LD. Therefore, the PGLD is a better fit for the data than the other four distributions. Since the PGLD has the smallest values of AIC, CAIC, and BIC, then this is considered another indicator of the superiority of the PGLD. Figure 4 shows the comparison graphically by plotting the empirical CDF against the CDF of PGLD, PLD, GLD, ELD, and LD. In Table 1, we note that the PLD, ELD, and LD (GLD, ELD, and LD) do not fit well the first (second) data set (based on the p-value ($< 0.05$)) but we use them for comparison purposes.

7.   Simulation study

In this section, a simulation study is performed to evaluate the performance of the estimation methods presented in Section 5. The MLEs, BELs and BEGs of the parameters $\alpha$, $\lambda$, $\theta$ and $c$ are computed and compared, through their mean squared errors (MSEs) and relative absolute biases (RABs), via a Monte Carlo simulation. The $95\%$ CIs, symmetric and HPD credible intervals are also computed and compared through their average interval lengths (AILs). The following algorithm is used to perform a simulation study:

$(1)$ For $i = 1, \dots, \hbar$, generate the values of progressive censoring with binomial removals, $r_{ij}$, such that $R_{ij}\sim$ Binomial$(n_{i}-w_{i}-\sum_{k = 1}^{j-1}R_{ik}, p)$ distribution, where $p$ is the removal probability, $j = 1, \dots, w_{i}-1$ and $R_{i w_{i}} \; = \; n_{i} \; - \; w_{i} \; - \; \sum_{j = 1}^{w_{i}-1} \; R_{ij}$.

$(2)$ For given values of the prior parameters ($\mu_{1}$, $\mu_{2}$, $\mu_{3}$, $\sigma_{1}$, $\sigma_{2}$, $\sigma_{3}$, $a_{1}$, $a_{2}$), generate values for the parameters ($\alpha$, $\lambda$, $\theta, c$).

$(3)$ For $i = 1, \dots, \hbar$, generate a progressively type-II censored sample of size $w_{i}$ from PGLD with CDF (2.9), according to the algorithm given in ^[32].

$(4)$ The MLEs, BELs and BEGs of the parameters $\alpha$, $\lambda$, $\theta$ and $c$ are computed as shown in Section 5. The BELs and BEGs of the parameters $\alpha$, $\lambda$, $\theta$ and $c$ are computed based on $\mathbb{M}( = 5500)$ MCMC samples and discard the first 500 values as burn-in period.

$(5)$ Repeat the above steps $\mathbb{N}( = 1,000)$ times.

$(6)$ If $\hat{\Theta}$ is an estimate of $\Theta$, then the average of estimates, MSE and RAB of $\hat{\Theta}$ over $\mathbb{N}$ samples are given, respectively, by

$(7)$ Calculate the average of estimates of the parameters $\alpha$, $\lambda$, $\theta$ and $c$ and their MSEs and RABs as shown in Step 6. Calculate also the average of the MSEs (MMSE) and the average of the RABs (MRAB) over the four parameter estimates.

$(8)$ Calculate the $95\%$ CIs (as shown in Subsection 5.1), symmetric and HPD credible intervals (as shown in Subsection 5.3) of the parameters and then calculate the average interval lengths (AILs) of them. Calculate also the average of the AILs (MAIL) over the four parameter estimates.

Tables 2–4 provide the computational results, which take into account the prior parameter values: $\mu_1 = 0.154$, $\mu_2 = 0.8$, $\mu_3 = -1.65$, $\sigma_1 = 0.8$, $\sigma_2 = 0.07$, $\sigma_3 = 0.3$, $a_{1} = 2.3$ and $a_{2} = 0.25$ to generate population parameter values: $\alpha = 1.6$, $\lambda = 2.3$, $\theta = 0.9$ and $c = 0.2$. Three distinct removal probabilities $p = 0.15, 0.55,$ and $0.95$ are considered.

The total number of observations $N$ is divided into two groups, $\hbar = 2$, and then into four groups, $\hbar = 4$, for comparison among MLEs, BELs, and BEGs.

● When there are two groups $(\hbar = 2)$,

● When there are four groups $(\hbar = 4)$,

The MLEs (BELs and BEGs) of the parameters $\alpha$, $\lambda$, $\theta$ and $c$ with their MSEs, RABs, MMSEs, MRABs are displayed in Table 2 (Table 3). The HPD credible intervals, symmetric credible intervals, and CIs with their AILs and MAILs are presented in Table 4.

We can observe from Tables 2–4 that:

$(1)$ Through the MMSEs and MRABs, the BELs and BEGs outperform the MLEs.

$(2)$ Through the MAILs, the HPD and symmetric credible intervals are superior to CIs. This confirms that BELs and BEGs are superior to MLEs.

$(3)$ The HPD is better than symmetric credible intervals via the MAILs.

$(4)$ For fixed $p$, $\hbar$, $N$, and $n_{i}$ (or $w_{i}$), $i = 1, \dots, \hbar$, by increasing $w_{i}$ (or $n_{i}$) the MSEs, RABs, MMSEs, MRABs, AILs and MAILs decrease. This assures that the more data we collect, the more accurate our results will be.

$(5)$ The BELs and BEGs are better at $\nu = 0.5$ than those at $\nu = -0.5$ through comparing the MMSE.

$(6)$ The BELs are better than BEGs at $\nu = 0.5$ but the converse is true at $\nu = -0.5$ through comparing the MMSE.

$(7)$ The MLEs are better at $p = 0.15$ than those at $p = 0.55, 0.95$ through comparing the MMSE.

$(8)$ Better results have been obtained at $\hbar = 2$ than those obtained at $\hbar = 4$ since the number of observations in the subgroups became smaller than those at $\hbar = 2$.

Except in a few rare cases, the above results are satisfactory, which could be related to data fluctuations.

Furthermore, if the hyperparameters are unknown, we can estimate them using past samples following the empirical Bayes approach, see ^[33]. Alternatively, the hierarchical Bayes technique, which uses a suitable prior for the hyperparameters, could be utilized, see ^[34].

8.   Conclusions

Due to the progress in manufacturing devices in the last decades, physicists and engineers may face a problem in assessing the lifetime distribution of these devices if they are connected in a certain mixed system such as a series-parallel system. To overcome this problem, we have introduced a new distribution called PGLD. The mentioned distribution can describe the lifetime distribution of series-parallel systems when the number of series subsystems, as well as the number of their parallel components, are random variables. The PGLD may arise by compounding two discrete (truncated Poisson and geometric) distributions with a mixture of continuous distributions, LD. Some important properties of the PGLD have been investigated, such as the $q$-th quantile, mode, $r$-th moment, mean residual lifetime, Bonferroni and Lorenz curves, Rényi and Shannon's entropies, PDF and CDF of the $i$-th order statistic. The progressive-stress model has been applied to units connected in a series-parallel structure. The lifetimes of these units under normal stress conditions are assumed to follow LD which is considered a mixture of exponential and gamma distributions. The progressive-stress model was used when the stress is increasing nonlinearly over time, and the inverse power law model had established a relationship between the stress and the proposed distribution's scale parameter. Based on progressively type-II censored data with binomial eliminations, two estimation methods were performed to estimate the unknown parameters. The Bayesian estimation was performed using two asymmetric (LINEX and GE) loss functions. CIs, symmetric, and HPD credible intervals for the unknown parameters were established. The numerical results showed that the Bayes estimates were performed well than the MLEs. An illustrative example, based on two real data sets, demonstrated the superiority of the proposed distribution over some other four distributions.

8.1.   Certain features and motivations to the PGLD

Finally, the features and motivations to the PGLD can be summarized as follows:

$(1)$ The CDF of the PGLD has a closed-form. This feature simplifies its use.

$(2)$ The four parameters included in the CDF of PGLD give it the ability to fit several data.

$(3)$ The CDF of PGLD includes the CDF of GLD, PLD, and LD as special cases.

$(4)$ The PGLD can represent the failure time of a series-parallel system. This is considered a motivation to introduce the PGLD.

$(5)$ The HRF of PGLD has a unimodal shape. This feature gives it more flexibility to fit and analyze several data arising from increasing and decreasing hazard rates.

$(6)$ The PGLD can represent the non-stationary data. This feature may be useful for the experimenter to forecast how some products would behave in different environments.

$(7)$ The PGLD is better to fit the data than some other distributions such as PLD, GLD, ELD, and LD.

$(8)$ Some other distributions may be emerged from Theorem 2.1 by choosing some other continuous distributions rather than LD.

Based on the above features and motivations, we hope that the PGLD will attract more attention from physicists and engineers in the near future.

8.2.   Future work

The following points may be investigated in future work:

$(1)$ The inference procedure may be implemented based on a general progressive censoring scheme such as hybrid progressive censoring.

$(2)$ Some other estimation methods such as the moments and probability weighted moments may be discussed.

$(3)$ Some other types of ALTs such as the step-stress ALT may be considered.

$(4)$ Prediction of future order statistics based on the PGLD may be investigated.

Acknowledgments

The authors would like to thank the Deanship of Scientific Research at Umm Al-Qura University for supporting this work by grant code (19-SCI-1-03-0009). They are also grateful to the editor and reviewers for their valuable comments and suggestions, which improved the paper.

Conflict of interest

The authors declare that they have no conflicts of interest regarding the publication of this article.