Lie symmetry analysis of fractional ordinary differential equation with neutral delay

1.   Introduction

In the analysis of survival and reliability data, different engineering and medical applications arise. Due to the limitations of cost and time, life-testing experiments are frequently performed under censoring schemes. In the literature contexts, Type-I and Type-II censoring schemes are initially used. In Type-I censoring, the life testing experiment is terminated at a pre-determined time T and the total number of failures is considered a random variable, while in Type-II censoring, the total number of failures is pre-fixed and in this case, the termination time is considered a random variable. The weakness of the traditional Type–I and Type-II censoring schemes is that they do not have the flexibility to remove some of the surviving units at different stages of the life-testing experiment. Recently, a generalization of Type-II censoring scheme called progressive Type-II censoring is the most commonly used. Suppose that n units are put on the life test at time zero for failure and the experimenter initially specified the number m of units that will be failed during the life test with the proposed censored scheme: $R = (R_1, R_2, ..., R_m)$ of random removals. The mechanism to perform the life test under the progressive Type-II censoring scheme is carried out as follows: At the time of the first failure, $x_{1:m:n}$, $R_1$ of the remaining $n-1$ surviving units are randomly removed from the test. Continuing on, at the second failure time $x_{2:m:n}$, $R_2$ of the remaining $n-R_1-1$ surviving units are randomly removed from the test and so on. This experiment is terminated at the time $x_{m:m:n}$ of the $m^{th}$ failure and at this time, the remaining $n-m-\sum_{i = m+1}^{n} R_i$ surviving units are all removed.

Balakrishnan ^[1] and Aggarwala and Balakrishnan ^[2] have introduced the concept of Type-II progressive censoring and an algorithm to simulate progressively Type-II censored samples from any continuous distribution. Recent statistical inferences for many lifetime distributions are also included in ^{[3,4,5,6,7,8,9,10]}.

In many real practical settings, experimenters need to control the termination time of the life test and to change the number of removed units during the test. For these reasons, Ng et al. ^[11] have proposed the adaptive progressive Type-II censoring scheme (APT-II). In this scheme, the expected value of the ideal time $T > \theta$ for terminating the life test is initially determined. Now, if ($x_{i:m:n}, i = 1, 2, \ldots, m$) are the $m$ completely observed failure times and($R_i, i = 1, 2, \ldots, m$) are the removed units at the times $x_{i:m:n}, i = 1, 2, \ldots, m$, respectively. APT-II censoring scheme is carried out as follows: if $x_{m:m:n} < T$, the test is terminated at $x_{m:m:n}$ and we will have the standard progressive Type-II censoring, in this case: $R_m = n-m-\sum_{i = 1}^{m}R_i$, otherwise, if $x_{j: m: n} < T < x_{j+1: m: n}$, where $j+1 < m$ and $x_{j:m:n}$ is the $j^{th}$ failure before the ideal time $T$, then any surviving units will be removed from the test by placing $R_{j+1} = R_{j+2} = \ldots = R_m = 0$, and at the time $x_{m:m:n}$ all the remaining units are removed, in this case: $R_m = n-m-\sum_{i = 1}^{j}R_i$. This scheme provides a design in which the experimenter is assured of getting m observed failure times for efficiency of statistical inference and the total test time will be closest to the ideal time T. Figure 1 illustrates a schematic presentation of the life testing experiment based on the APT-II censoring scheme.

Many authors have considered classical and Bayesian estimations for the parameters of different distributions based on APT-II. Examples of many others, under the APT-II censoring scheme, are: Chen and Gui ^[12] discussed the estimation problem of a lifetime distribution with a bathtub-shaped failure rate function, Abu-Moussa et al. ^[13] discussed the maximum likelihood estimation for the Gompertz distribution. Mahmoud et al. ^[14] studied the Bayes estimation of the Pareto distribution. Asadi et al. ^[15] investigated the accelerated life test in Gompertz distribution for the evaluation of virus-containing micro-droplets data. AL Alotaibi et al. ^[16] considered the estimation of parameters of the XLindley distribution, Almalki et al. ^[17] have explored the partially constant-stress accelerated life tests model for estimation of the parameters of the Kumaraswamy distribution and Ahmad et al. ^[18] adopted both Bayesian and non-Bayesian inference for the exponentiated power Lindley distribution.

Using the power transformation method (see Shukala and Shanker ^[19]), Bhat and Ahmad ^[20] introduced a new form of Rayleigh distribution called the Power Rayleigh distribution (PRD) with probability density (PDF), cumulative distribution (CDF) and the reliability functions given, respectively, by:

The hazard rate function (HRF) of the PRD is given by

As a generalization of the standard Rayleigh distribution, which can be used in modeling the lifetimes of many pieces of equipment like resistors, transformers and capacitors (see Kundu and Raqab ^[21] and Ateeq et al. ^[22]), the hazard rate function of PRD is constant when $\beta = 0.5$, increasing when $\beta < 0.5$ and decreasing when $\beta > 0.5$. For this flexibility, PRD can be utilized for more reliability and survival data analysis. Particularly, as shown by Bhat and Ahmad ^[20], this new distribution has superiority for modeling the stress and strength data. The PDF and the HRF of the PRD with different combinations of $(\theta, \beta)$ are presented in Figure 2.

Under the APT-II censoring scheme. Many authors have adopted statistical inference for the Rayleigh distribution and some of its generalizations. Recent studies were completed by Mahmoud and Ghazal ^[23] for the Exponentiated Rayleigh distribution, Almetwally et al. ^[24] for the Generalized Rayleigh distribution, Panahi and Moradi ^[25], Gao and Gao ^[26] and Fan and Gui ^[27] for the Inverted Exponentiated Rayleigh distribution. To the best of our knowledge, there exists no statistical inference about the parameters of the PRD based on any censoring scheme. The aim of this paper is to focus on classical and Bayesian point and interval estimation of the unknown parameters of the PRD using the APT-II censored data. This will extend the application of PRD as a life model to include more reliability and survival data.

The rest of this paper is organized as follows: Section 2 is devoted to point and interval estimation using the likelihood approach. Bayesian estimation under the squared error (SE) and the general entropy (GE) loss functions is explored in Section 3. In Section 4, both point and interval estimates are presented and the Bootstrap sampling method is also implemented. To evaluate the performance of the obtained point and interval estimates, a Monte Carlo simulation study is conducted in Section 5. Section 6 illustrates the application of the developed estimation procedure to a set of real-life data. Finally, a conclusion about this work is summarized in Section 7.

2.   Maximum likelihood estimation

Based on an APT-II observed censored sampled data from a continuous distribution with PDF $f\left(x; \varphi\right)$ and CDF $F\left(x; \varphi\right)$, Ng et al. ^[11] derived the likelihood function for a vector of parameter $\varphi$ in the form:

where: $x_i = x_{i, m, n}\, \ c_j = n-m-\sum_{i = 1}^{j}R_i$ and $C$ is a constant that is independent of the parameters $\varphi$. Substituting the corresponding formulas for $f\left(x_i; \theta, \beta\right)$, $F\left(x_i; \theta, \beta\right)$ of PRD from Eqs (1.1) and (1.2) in Eq (2.1). This implies that the corresponding likelihood function of the parameters $\theta$ and $\beta$ using APT-II data is given by:

Therefore, the natural logarithmic function of the likelihood function is

Differentiating both sides of Eq (2.2) with respect to $\theta$ and $\beta$ and equating the derived equations to zero, we get the likelihood equations:

When the shape parameter $\beta$ is known, solving Eq (2.4) with respect to $\theta$, the maximum likelihood estimate (MLE) of $\theta$ exists in an explicit form given by

When the two parameters are both unknown, numerical methods should be used to solve the system of the nonlinear equations in Eqs (2.4) and (2.5). In this paper, a multivariate Newton-Raphson numerical method is implemented to obtain the approximated values for the MLEs ($\hat{\theta}$, $\hat{\beta}$) of $\theta$ and $\beta$. Utilizing the invariance property of the MLEs, the MLE of the reliability function at a given time t is given by:

Under the APT-II censoring scheme, another desired statistical inference for the PRD is to find confidence interval estimates for the unknown parameters $\theta$ and $\beta$. It is mathematically difficult to construct exact confidence intervals because the distributions of the MLEs of $\theta$ and $\beta$ are not in closed forms. The observed Fisher information matrix of ($\theta$, $\beta$) can also be used to find approximated confidence intervals for the two unknown parameters (see Lawless ^[28]) which are defined as:

where:

Hence, the approximated variance-covariance matrix of ($\hat{\theta}$, $\hat{\beta}$) is simplified by

Consequently, based on the large sample theory of the MLEs, the asymptotic $100(1-\alpha)\%$ confidence intervals (ACIs) for $\theta$ and $\beta$ are constructed respectively as:

where $z_\frac{\alpha}{2}$ is the $(1-\alpha)\%$ quantile of the standard Normal distribution, and ${\hat{v}}_{11}$, ${\hat{v}}_{22}$ are the diagonal elements of $\hat{v}$. Accordingly, the approximate $100(1-\alpha)\%$ confidence intervals for the reliability function at a given time t is

3.   Bayesian estimation

In statistical inference, the Bayesian estimation approach is a complement part to the conventional likelihood estimation approach. To apply the estimation procedures using the Bayesian approach, prior information of the parameters $\theta$ and $\beta$ should be determined first. When no previous information is available about these parameters, it is convenient to specify their priors in such a way that their joint prior is proportional to their likelihood function. In this research, the joint prior distribution of the parameters $\theta$ and $\theta$ is proposed to take the form:

where $\pi_1(\theta)$ is the Inverse Rayleigh distribution with scale parameter $\gamma$ and $\pi_2(\beta|\theta)$ is the Exponential distribution with scale parameter $(\frac{\lambda}{\theta^2})$ with prior densities given, respectively, by:

Substituting for $\pi_1(\theta)$ and $\pi_2(\beta|\theta)$ from Eqs (3.2) and (3.3) in Eq (3.1), the joint prior of $\theta$ and $\beta$ has the form

Now, combining Eq (2.1) and (3.3), the posterior functions of $\theta$ and $\beta$ is:

where: $u = \sum_{i = 1}^{m}x_i^{2\beta}+\sum_{i = 1}^{j}{R_ix}_i^{2\beta}+c_jx_m^{2\beta}$. Clearly, the joint prior given in Eq (3.4) is proportional to the posterior function of $(\theta, \beta)$ in Eq (3.5) and thus, it can be considered a conjugate prior. The conditional posteriors of $\theta$ and $\beta$ can be obtained, respectively, as:

3.1.   Bayesian estimation when the shape parameter is known

When the shape parameter $\beta$ is known (As in the case of the standard Rayleigh distribution). The posterior function $\pi^\ast\left(\theta, \beta\middle| x\right)$ given in Eq (3.4) becomes a function of $\theta$ only given by

Consequently, the Bayes point estimate of $\theta$ under the squared error loss function (SE) is obtained as the mean of the posterior, which is given by:

Similarly, the Bayes point estimate of $\theta$ under the general entropy loss function (GE) indexed with parameter $\alpha > 0$ is obtained as:

When the shape parameter is known, another common Bayesian inference is to find a credible interval in the form $\left[l_l, l_u\right]$ for the parameters $\theta$. The $\left(1-\alpha\right)\%$ credible interval for the parameter $\theta$ can be obtained as

Substituting for $\pi^\ast\left(\theta\middle|\beta, x\right)$ from Eq (3.5) in Eq (3.10) and simplifying the integral in the above equation, the lower and the upper bounds of the $\left(1-\alpha\right)\%$ credible interval for $\theta$ can be obtained from the following equation:

where: $z = ({\frac{u}{2}+\gamma}^2+\lambda\beta), and I_g(x, y) = \frac{1}{\Gamma(y)}\int_{0}^{x}{t^{y-1}e^{-t}\partial t}$ is the incomplete Gamma function.

In addition to the constraint Eq (3.12), a $\left(1-\alpha\right)\%$ highest posterior credible intervals (HPD) credible interval for $\theta$ can also be obtained by adding the constraint:

3.2.   Bayesian estimation when both parameters are unknown

When the parameters $\theta$ and $\beta$ are both unknown, it is difficult to obtain explicit forms of their Bayes point estimates under SE or GE loss functions because, in this case, the corresponding integrals:

cannot be found in closed forms for any number $\alpha\in R$.

To overcome this problem, the Markov Chain Monte Carlo (MCMC) method can be used to approximate the above integrals. In this paper, we use the common type of MCMC, which is called the Metropolis-Hasting (M-H) developed by Metropolis et al. ^[29] and extended by Hasting ^[30] with Gibbs sampling to find Bayes point and interval estimates for the parameters $\theta$ and $\beta$. Assuming normal is the target distribution, the procedure of implementing the MCMC is performed through the following steps:

Step 1. Set the initial values of $\theta$ and $\beta$, $(\theta^{(0)}, \, \beta^{(0)})$.

Step 2. Put j = 1.

Step 3. Simulate $\theta^\ast$, $\beta^\ast$ from the normal distributions N$(\theta^{(j-1)}, \hat{v}_{11})$, N$(\beta^{(j-1)}, \hat{v}_{22})$ respectively, where ${\hat{v}}_{11}$ and ${\hat{v}}_{22}$ are the diagonal elements of the variance-covariance matrix computed given the values of $\theta^{(j-1)}, \, \beta^{(j-1)}$.

Step 4. Compute the acceptance ratios:

Step 5. Generate $u_1$, $u_2$ from the uniform distribution $U(0, 1)$.

Step 6. If $u_1\le r_{\theta^\ast}$, $\ u_2\le r_{\beta^\ast}$, set: $\theta^{(j)} = \theta^\ast, \beta^{(j)} = \beta^\ast$, else: $\theta^{(j)} = \theta^{(j-1)}, \beta^{(j)} = \beta^{(j-1)}$.

Step 7. Use the generated sample to obtain an estimate of the reliability function: $R_{(t)}^{(j)}$.

Step 8. Set $j = j+1$.

Step 9. Repeat Step 3- Step 8 $N$ times to have $\left(\theta^{\left(1\right)}, \theta^{\left(2\right)}, \ldots, \theta^{\left(N\right)}\right), \, \left(\beta^{\left(1\right)}, \beta^{\left(2\right)}, \ldots, \beta^{\left(N\right)}\right)$ Then the Bayesian estimates of $\theta$ and $\beta$ under SE are given, respectively, by

And under the GE, the Bayesian estimates of $\theta$ and $\beta$ are given, respectively, by

where $\alpha$ is the parameter of the GE function and $M$ is the number of iterations before the stationary distribution is achieved (the burn in-period) from $N$ (the total number of iterations).

Based on the (M-H) algorithm, the Bayesian estimations of the reliability function under SE and GE at a given time $t$ are:

Ascending order the values of $\theta^{\left(i\right)}$, $\beta^{\left(i\right)}$ and $R_{(t)}^{(i)}$, $i = M+1, \ldots, N$ generated by the M-H algorithm. The $100\left(1-\alpha\right)\%$ symmetric two-sided credible intervals of $\theta$, $\beta$ and the reliability function $R(t)$ are given, respectively, by:

Using the method of Chen and Shao ^[31]. The symmetric $100\left(1-\alpha\right)\%$ two sided HPD credible intervals of $\theta$, $\beta$ and the reliability function $R(t)$ are given, respectively, by:

where $R^\ast\left(t\right)$ is the estimated $R(t)$ by the (M-H) algorithm and $i^\ast$ is chosen such that:

where $[y]$ denotes the greatest integer less than or equal $y$.

4.   Estimation using the Bootstrap technique

One common method to estimate the parameters of the underlying distribution is the Bootstrap resampling technique. This method gives better estimates of the parameters than the approximate Bayes and likelihood estimates when the effective sample $m$ is small enough. Details about the Bootstrap estimation method are explored in Hall ^[32] and Efron and Tibshirani ^[33]. The procedure to obtain point and interval estimates for the parameters $\theta$ and $\beta$ by this method is described in the following steps:

Step 1. Based on the specified values of $n$, $m$, $T$, $R$, ncompute ${\hat{\theta}}_{mle}$ and ${\hat{\beta}}_{mle}$.

Step 2. Using ${\hat{\theta}}_{mle}$ and ${\hat{\beta}}_{mle}$ generate a Bootstrap resample.

Step 3. Calculate the Bootstrap estimates ${\hat{\theta}}_p^1$ and ${\hat{\beta}}_p^1$.

Step 4. Repeat the steps 2 and 3 up to $k$ times to have ${\hat{\theta}}_p^1, {\hat{\theta}}_p^2, \ldots, {\hat{\theta}}_p^k, \, {\hat{\beta}}_p^1, {\hat{\beta}}_p^2, \ldots, {\hat{\beta}}_p^k$.

Step 5. Use the generated sample to obtain an estimate of the reliability function: $R(t)_{p}^{(1)}, R(t)_{p}^{(2)}, \ldots, \ R(t)_{p}^{(k)}$.

Step 6. Compute the Bootstrap point estimates

Step 7. Rearrange the Bootstrap estimates in ascending order to have $100(1-\alpha)\%$ percentile confidence intervals for $\theta$, $\beta$ and the reliability function $R(t)$, respectively, as:

5.   Monte Carlo simulation study

We evaluate the performance of the point and interval estimates obtained in the previous sections. Thus, Monte Carlo simulations are performed by generating $1000$ APT-II samples from the PRD with different settings of sample size, affective sample size, true parameter values and ideal termination time with the following censoring schemes:

SC.I presents left censoring at the beginning of the experiment, SC.II presents ordinary Type-II censoring and SC.III presents progressive APT-II at the expected mid-termination time of the experiment.

The simulation process is implemented through the following steps:

$1)$ Determine the values of $n$, $m$, $R$, $T$ and set the true parameter values ($\theta$, $\beta$).

$2)$ Generate a random sample of size $n$ from the uniform distribution and apply the Balakrishnan and Sandho ^[34] algorithm to generate a random sample of size $m$ from the PRD using the relation:

$3)$ Determine the value of $j$ and remove all observations greater than $x_j$.

$4)$ Generate $(m\ -\ j\ -\ 1)$ Type-II censored sample from $[\frac{f(x_j)}{1\ -\ F(x_{j+1})}]$ where f(x) and F(x) are the PDF and the CDF of the PRD, and stop the experiment at $x_m$.

$5)$ Obtain point and interval estimates of the parameters $\theta$, $\beta$ and the reliability function at time $t = 1.5$.

$6)$ Repeat steps 2–5 and compute the average absolute biases (ABs), the mean squared errors (MSEs) of the point estimates, the average lengths (ALs) and the coverage probabilities (CPs) of the interval estimates.

For the case where the shape parameter $\beta$ is assumed to be known, the true parameter values are taken to be ($\theta = 1$, $\beta = 0.8$), the values of $n$ and $m$ are determined such that:

with the ideal terminations $T = 1$ and $T = 2$.

When both parameters are assumed to be unknown, the true parameter values are considered to be $\left\{(\theta, \beta) = (2, 0.4)\text{ and }(2, 0.8) \right\}$ with ideal termination times ($T = 200 \text{ for }(\theta = 2, \beta = 0.4)$) and ($T = 20\text{ for }(\theta = 2, \beta = 0.8)$), the sample sizes are proposed to be $n = 80$, n = 50 and the effective sample size $m$ is indicated such that the proportion of observed failures $r = (\frac{m}{n})$ are: $0.80$, $0.60$ and $0.40$.

In the Bayesian estimation, the values of the hyperparameters ($\gamma$, $\lambda$) are chosen in a way such that the prior mean is equal to the expected value of the unknown parameter; see Kundu ^[35]. The shape parameter of the GE loss is fixed to be equal to $\alpha = 1.1$. Bayesian point estimates and associated credible intervals are computed by generating $12000$ samples using the M-H algorithm and discarding the first $2000$ corresponding values of ($\theta$, $\beta$) as burn-in periods. Simulation results are listed in Tables 1–7. Tables 1–7 manifest the following conclusions:

$1)$ All the obtained (point and interval) estimates have satisfactory behavior and gave an efficient inference about the unknown parameters of the PRD and the reliability function.

$2)$ All estimates perform better with an increase in sample size $n$, the proportion of effective sample size $r$ and the ideal time $T$.

$3)$ The point and interval estimates obtained by the Bootstrap sampling method are clearly better than their corresponding estimates obtained by the conventional likelihood method.

$4)$ Under the same settings, the Bayes estimates obtained by the MCMC method are better than estimates obtained by the Likelihood and Bootstrap classical approaches.

$5)$ We notice that the obtained estimates become better with a decrease in the value of the shape parameter.

$6)$ It is noted that the CPs of the obtained interval estimates are almost close at a $95\%$ confidence level, but when the lengths of (Asymptotic, Bootstrap or Credible) interval estimates become relatively small enough, the CPs of the interval estimates are strictly less than their nominal $95\%$ confidence level.

$7)$ In comparing schemes 1–3, we conclude that:

$(i)$ Point estimates obtained under the progressive SC.II scheme are better than point estimates obtained under progressive (SC.I and SC.III) censoring schemes.

$(ii)$ Interval estimates obtained under SC.III have smaller lengths than interval estimates obtained under progressive (SC.I and SC.II) schemes.

6.   Application to real data

We illustrate the application of the developed estimation methods. The inferential procedure is investigated in a real data set from Al-Aqtash et al. ^[36] This data set consists of $66$ observations representing the breaking stress of carbon fibers of $50$ mm length (GPa). Bhat and Ahmad ^[20] have analyzed this data to confirm the priority fit of the PRD in real-life data as compared to other generalizations of Rayleigh and Weibull distributions. Assuming the PRD modeling for this data, the MLEs of the parameters $\theta$ and $\beta$ using the complete data set (see Bhat and Ahmad ^[20]) are ${\hat{\theta}}_{com} = 4.850$, ${\hat{\beta}}_{com} = 1.721$, the asymptotic 95% confidence intervals for $\theta$ and $\beta$ are $(2.818, 6.883)$ and $(1.396, 2.045)$. The Kolmogorov-Smirnov (K-S), the Akaike's information and the negative log-likelihood goodness of fit test statistics are $0.08$, $172.14$ and $176.33$, respectively, with a p-value = $0.768$. Figures 3–6 present the empirical and the theoretical CDF, the histogram plot, the empirical survival function and the empirical cumulative hazard function of these data.

To confirm the validity of the PRD for the APT-II censored data, we generate APT-II censored data from the original data with $m = 40$ and $T = 3.5$ using two different censoring schemes as presented in Table 8. Based on the APT-II censored data, the performance of point estimates $\hat{\theta}$ and $\hat{\beta}$ of $\theta$ and $\beta$, respectively, are measured in terms of their relative efficiencies defined as:

The relative efficiencies of the point estimates and the ABLs of the interval estimates are listed in Table 9. Based on ($2, 0^{\ast20}, 2, 0^{\ast20}$) and ($8, 0^{\ast25}, 8, 0^{\ast15}$) APT-II censoring schemes, Figures 7 and 8 present the trace and the histogram plots for point estimates of ($\theta$, $\beta$) using 12000 samples generated by the MCMC method and assuming M$= 2000$ is the burn-in period. As clearly observed from Table 9, all point estimates have high relative efficiency and The Bayes point estimates under GE are the best among others as they have the highest relative efficiency. All estimates became better with a decrease in the sum of random removals and the interval estimates are also very close to the approximate intervals for ($\theta$, $\beta$) obtained by the asymptotic distributions using the complete data set. It is clearly observed from Figures 7 and 8 that Bayesian estimates obtained by the MCMC converge to their means and thus the MCMC is recommended.

7.   Conclusions

In this paper, based on the Type-II censoring, we have considered point and interval estimation for the parameters and the reliability function of PRD via likelihood, Bayesian and Bootstrap estimation methods. The Bayesian estimation approach is carried out under the SE and GE loss functions. When the shape parameter is known, the ML and Bayesian point estimates of the scale parameter, as well as the corresponding asymptotic confidence, credible and HPD credible intervals, are obtained in explicit forms. When both parameters are unknown, the MLEs and the asymptotic confidence intervals of the unknown parameters are obtained using the Newton-Raphson numerical method, while Bayesian point estimators and associated credible and HPD credible intervals are obtained based on the MCMC method.

To evaluate the performance of the obtained point and interval estimates, a simulation study is conducted in different combinations of the life testing experiment, and the estimation procedure is also applied to a real-life data set. Concluding results indicated that both point and interval estimates are efficient and have good characteristics for estimating the model parameters with a specific preference for the Bayesian point and interval estimates obtained under the GE loss function.

As a generalization of the standard Rayleigh distribution, the main conclusions conducted from this work confirmed that the PRD is an appropriate lifetime model and can be effectively used with Type-II censored data. Future researchers may extend this work by employing PRD with other censoring schemes.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The authors declare no conflict of interest.