A new least squares method for estimation and prediction based on the cumulative Hazard function

1.   Introduction

The main objective of statistical inference is to draw conclusions about a certain population based on a random sample of that population. Conclusions may be about the functional form of the distribution for the population of interest, estimates of its characteristics, testing statistical hypotheses regarding parameters or making predictions about future events based on current knowledge.

One of the essential elements of reliability analysis is testing the product failure times under normal operating conditions. Products are frequently highly dependable, and extensive applied testing is prohibitively costly. The use of complete samples is not a good choice in statistical analysis when the product is highly reliable and applied testing is prohibitively expensive. In such situations, statisticians have recommended many types of censoring to save time and money. Type Ⅱ-right censorship is one of these censoring types and has been more popular over the past few decades. It is expected for this type that a few lower-order statistics are observed and that subsequent inferences are then required, such as determining the distribution from which the data is collected, estimating the distribution parameters and making predictions about unobserved events. Here, the emphasis is on estimate for censored samples as well as prediction problems.

There are many contributions to the estimation of complete and censored data. Excellent references for this subject, among others, include Balakrishnan and Cramer ^[1], Casella and Berger ^[2], Lehmann and Casella ^[3] and the references therein.

Prediction of future observations is a major concern in many real-world problems. The prediction for ordered random variables (RVs) is frequently used in industrial applications and survival studies to estimate the future prevalence of defective products. In lifetimes testing experiments, the interval and point predictions are also helpful in determining the best censoring strategy. To put it more precisely, we can place identical items in a given life-test experiment and wait until a manageable number of failed items are obtained in relation to cost and time. We can therefore predict the failure times of the survivor items using these observed failure times. This prediction enables the experimenter to choose an appropriate censoring scheme. Such censoring schemes may be Type Ⅰ or Type Ⅱ censoring, progressive Type Ⅱ censoring, or hybrid censoring. The prediction has additional benefits that help us decide whether or not the life test needs to be accelerated and when a future failed item should cause the life test to end. Many authors have considered the prediction issue in statistical literature from both a theoretical and an applied standpoint. Ahsanullah ^[4] studied the linear prediction of record values of the two-parameter exponential distribution. Bayesian predictions of both ordinary order statistics (OOSs) and generalized order statistics (GOSs) were studied in AL-Hussaini ^[5] and AL-Hussaini and Ahmed ^[6], respectively. In Aly ^[7], two-point predictors of the fractional kth upper record values from the exponential distribution are given. Some predictive and reconstructive results of dual GOSs from the inverse Weibull distribution are obtained in Aly ^[8] via the pivotal quantities approach. The problem of predicting future lifetimes from the Weibull probability model for a simple step-stress plan under the Khamis-Higgins model is studied by Amleh and Raqab ^[9]. Barakat et al. ^[10] suggested a new method for predicting future order statistics based on some of their distributional properties. Moreover, by applying the cumulative hazard transformation, Barakat et al. ^[11] constructed nonparametric prediction intervals for GOSs, and their exact coverage probabilities were determined. For all distributions that are attracted to the Weibull distribution, El-Adll et al. ^[12] have constructed different asymptotic prediction intervals for future lower extreme order statistics. The maximum likelihood predictor was first introduced by Kaminsky and Rhodin ^[13]. In life-testing experiments with fixed sample sizes, Lawless ^[14] suggested a pivotal quantity for constructing a prediction interval for an order statistic. Lingappaiah ^[15] suggested a different pivotal quantity for the same purpose. Raqab ^[16] obtained an optimal prediction interval for $m$-GOSs based on a one-parameter exponential distribution. Concise reviews that are related to this topic and include several prominent prediction results may be found in the books by Aitcheson and Dunsmore ^[17], David and Nagaraja ^[18], Geisser ^[19] and Kaminsky and Nelson ^[20].

The main contribution of the paper is presented in Section 3. More specifically, based on Rényi representation and its extensions, a new least squares method for estimation and prediction is introduced, relying on the cumulative hazard function. The proposed method is formulated in a general manner so that it can be applied to any model of ordered random variables. Furthermore, the method is used for complete or type Ⅱ right-censored data and simplifies the computations significantly for those distributions with an explicit form for the cumulative hazard function.

The rest of this paper is organized as follows: In Section 2, some necessary preliminary results that are related to the present work are presented. A new least-squares method for estimation and prediction is presented in Section 3. The methodology presented in Section 3 is validated in Section 4 through comprehensive simulation studies. In Section 5, two real data sets are analyzed for illustrative and comparison purposes.

2.   Preliminary results

As a unified model of various ascending models of ordered RVs, Kamps ^[22] established the GOSs, based on any continuous distribution function (CDF) $F$, via their joint probability density function (PDF)

on the cone $F^{-1}(0+) < x_1\leq \cdots\leq x_n < F^{-1}(1),$ where $C_{r} = \prod \limits_{i = 1}^{r} \gamma_{i}, \; r = 1, 2, ..., n,$ $F^{-1}$ denotes the quantile function of $F.$ The model parameters are defined by the vector $\underline\gamma = (\gamma_{1}, ..., \gamma_{n}),$ where $\gamma_{r} = k+\sum \limits_{i = r}^{n-1} (m_i+1) > 0, \; i = 1, 2, ..., n, \; \mbox{with}\; \gamma_{n} = k > 0\; \mbox{and}\; (m_1, ..., m_{n-1}) \in \mathbb{R}^{n-1}.$ Particular models, such as OOSs, sequential order statistics (SOSs), progressive type Ⅱ censored order statistics, standard $k$th record values and Pfeifer's record, result from specific choices of the model parameters, $\gamma_1, ..., \gamma_n$. When $m_i = m$ for all $i = 1, 2, ... n,$ extensive asymptotic results for bootstrapping $m$-GOSs are obtained by Barakat et al. ^[21].

According to the restriction $\gamma_{i} \neq \gamma_{j}$ for $i \neq j$, a wide subclass of GOSs-model, excluding record values, have discussed by Kamps and Cramer ^[23], where some important distributional properties have been considered. In particular, the PDF of the $r$th GOS $X^{*}_{r:n}$ is given by

and the CDF of $X^{*}_{r:n}$ is

where

In view of the probability integral transform, the RVs, $U^{*}_{r:n} = F^{-1}(X^{*}_{r:n}), \; r = 1, 2, ..., n,$ are uniform GOSs. Consequently, it can be shown that

and

For the uniform OOSs, $U_{r:n}, \; r = 1, 2, ..., n,$ we have

and

The ordinary least squares method of estimation for complete samples was originally proposed by Swain et al. ^[24]. The method is based on minimizing the function

with respect to the unknown parameter vector $\Theta = (\theta_1, \theta_2, ..., \theta_\ell),$ where ${\tilde x}_n = (x_1, x_2, ..., x_n)$ is an observed ordered sample. The weighted least squares estimates can be accomplished by minimizing the function

with respect to $\Theta.$ Several authors used the least squares method and the weighted least squares method. For estimating the parameters from different distributions; among them are Gupta and Kundu ^[25] and Kundu and Raqab ^[26].

El-Adll and Aly ^[27] have extended the above results to type Ⅱ censored samples. Namely, based on the first $r$ observed OOSs, $x_1, x_2, ..., x_r,$ approximate least square estimates of $\Theta$ can be obtained by minimizing the least squares function

For the GOSs model, we can exptend the ordinary least squares and weighted least squares functions, respectively, to take the formulas

and

where

${\tilde x}_r = (x_1, x_2, ..., x_r)$ are observed values of the GOSs $X_{1:n}^{*}, X_{2:n}^{*}, ..., X_{r:n}^{*},$ $E[U_{i:n}^{*}]$ and $Var(U_{i:n}^{*})$ are given by (2.4) and (2.5), respectively.

Kamps ^[22] extended Rényi's representation ^[29], to GOSs model (see also Barakat et al. ^[28]). This representation enables us to express the $r$th GOS based on the exponential distribution as a linear combination of independent and identically distributed (iid) RVs from the EXP(1). Namely,

where $Z_1, ..., Z_n$ are iid RVs from the EXP(1) and "$U{ \begin{array}{c} \stackrel{d}{ = }\\ {}\end{array}}V$" means that the RVs $U$ and $V$ have the same CDF.

3.   A new least squares method for estimation and prediction

Minimizing the least squares and weighted least squares functions in (2.11) and (2.12) is not always a simple problem. In this section, a novel and efficient method for solving such problems is proposed.

3.1.   Least squares estimation via the cumulative hazard function

The cumulative hazard function of any RV $X$ with CDF $F$ is defined by

Evidently, the function $H(x)$ is a nonnegative and nondecreasing function on $x.$ Suppose now that $X^{*}_{i:n}: = X(i, n, \tilde{m}, k)$ denotes the $i$th GOS from a continuous CDF $F$. Therefore, the RV's

are GOSs from the standard exponential distribution (denoted by EXP(1). Clearly, if the CDF $F$ depends on an unknown vector of parameters $\Theta = (\theta_1, \theta_2, ..., \theta_\ell),$ with $\ell \geq 1,$ then $H$ depends on the same vector of parameters. Hence, we can develop a least square method based on $H$. In view of Kamps ^[22] and Barakat et al. ^[28], we have

Consequently, the parameters can be estimated by minimizing the least squares function

where $-\infty < x_{1} < x_{2} < ... < x_{n} < \infty$ are observed values of the GOSs $X^{*}_{1:n}, X^{*}_{2:n}, ..., X^{*}_{n:n}.$ Moreover, the weighted least squares estimates based on the cumulative hazard function can be obtained by minimizing the weighted least squares function.

According to the extended Rényi's representation (2.13), we have

where $z_1, ..., z_n$ is an observed random sample from the EXP(1). Clearly, $E\left[\sum_{j = r+1}^{i} \frac{Z_j-1}{\gamma_j}\right] = 0.$ Therefore, we can approximate the sum of the last term by its mean. Consequently, we can suggest the modified least squares function.

An approximate least squares estimate of $\Theta$ based on the first $r$ observed GOSs, ${\tilde x}_r: = (x_{1}, x_{2}, ..., x_{r})$ for $r \leq n,$ can be obtained by minimizing $\mathcal{L}^{*}_{H, r}(\Theta|{\tilde x}_r)$ with respect to $\Theta.$ Similarly, minimizing the function

produces a modified weighted least squares estimate of $\Theta$ based on $x_{1}, x_{2}, ..., x_{r}.$

3.2.   A least squares predictor relying on GOSs

The results of El-Adll and Aly ^[27] can be extended to GOSs through minimizing the predictive least squares function

with respect to $\Theta$ and $x_s.$ Similarly, a weighted least squares predictor can be derived by choosing the weights $w_i$ defined by (2.12). As we proceeded in (3.5) and (3.6), an approximate point predictor of the unobserved $s$th GOS can be obtained via the minimization of the proposed predictive least squares function

By the same manner, approximate weighted least squares estimates of $\Theta$ and $x_s$ based on $x_{1}, x_{2}, ..., x_{r}$ are derived via minimizing the predictive weighed least squares function

Remark 3.1.

One advantage of choosing the cumulative hazard function transformation is that it always follows the standard exponential distribution. In addition, its mean and variance can simply be computed for any model of ordered random variables, and they do not depend on unknown parameters. Moreover, as a quick comparison between the proposed method and one of the most well-known methods of estimation, we find that the Maximum likelihood estimation (MLE) for the distribution parameters may be difficult to obtain in certain cases. Particularly when the support of the distribution is unknown. Also, the MLEs may not be robust enough to depart from the assumed distribution. These considerations motivated the least-squares approach to be used.

4.   Numerical simulation experiments

In this section, simulation experiments are carried out to compare the proposed method with different estimation and prediction techniques. We are mainly interested in some important probability distributions that can be widely applied in survival analysis, reliability theory and life-testing experiments. For brevity, we compare only three methods, namely the maximum likelihood, the ordinary weighted least squares and the weighted least squares via the cumulative hazard function. The following assumes that the first $r$ ($r\leq n$) elements of the GOSs based on a continuous distribution $F$ are observed and used to estimate the unknown distribution parameters and predict some future observations. In this simulation, 10,000 independent random samples are generated from $F$ and then used in the proposed estimation and prediction methods.

4.1.   Estimation

In view of Kamps ^[22], the likelihood function based on the first $r$ elements of the GOSs, which is the joint PDF of $X_{1:n}^{*}, X_{2:n}^{*}, ..., X_{r:n}^{*},$ is given by

The maximum likelihood estimates (MLEs) of the unknown parameters $\Theta = (\theta_1, \theta_2, ..., \theta_\ell)$ can be accomplished via maximizing $\mathcal{L}^{*}(\Theta|{\tilde x}_r)$ by solving the nonlinear equations $\frac{\partial \log (\mathcal{L}^{*}(\Theta|{\tilde x}_r))}{\partial \theta_j} = 0,$ $\; j = 1, 2, ..., \ell,$ numerically. The ordinary weighted least squares estimates (WLSEs) are obtained by minimizing the function $\mathcal{{WL}}_{F, r}^{*}(\Theta|{\tilde x}_r)$ given in (2.12), through solving the nonlinear equations $\frac{\partial \mathcal{{WL}}_{F, r}^{*}(\Theta|{\tilde x}_r)}{\partial \theta_j} = 0, \; j = 1, 2, ..., \ell,$ numerically. Finally, the modified weighted least squares estimates (MWLSEs) can be obtained by minimizing the function $\mathcal{{WL}}_{H, r}^{*}(\Theta|{\tilde x}_r)$ in (3.7), via solving the nonlinear equations $\frac{\partial \mathcal{{WL}}_{H, r}^{*}(\Theta|{\tilde x}_r)}{\partial \theta_j} = 0, \; j = 1, 2, ..., \ell,$ numerically. The OOSs model, where $\gamma_i = n-i+1,$ is primarily the focus of the next sections of the study.

We use the root mean square error to compare different estimators or predictors. As a result of its sensitivity to extreme values, the root-mean-square error may not even exist. Pitman's measure of closeness (Pitman ^[30]) is also used for comparing estimators or predictors. According to Keating et al. ^[31], who provided several inspiring instances and examples, Pitman's measure of closeness is an efficient criterion for selecting among estimators. Pitman's measure of closeness is widely applied to assess estimators and predictors by several authors. Pitman's measure is used by Balakrishnan et al. ^[32] for comparing different point predictors for type Ⅱ censored data that follows an exponential distribution in one sample and two sample cases. Nagaraja ^[33] used the mean square error and Pitman's measure of closeness to compare the best linear predictor with the best linear invariant predictor for the record value and order statistic. Raqab et al. ^[34] compared different point predictors of progressively censored units using Pitman's measure of closeness.

Let $\hat{\theta}_1$ and $\hat{\theta}_2$ be two estimators of the same parameter $\theta$. Then, $\hat{\theta}_1$ is Pitman closer than $\hat{\theta}_2$ if

for all values of $\theta$, with strict inequality for at least one value of $\theta$.

Pitman's measure of closeness is applied to explore which of two estimators is the closest (in probability) to the true value of a parameter.

4.1.1.   Two parameter-exponential distribution

If the CDF of the RV $X$ is given by

it is said to have a two-parameter exponential distribution, with a location parameter $\theta$ and a scale parameter $\sigma$. The MLEs of $\theta$ and $\sigma$ based on the first $r$ elements of GOSs are given in El-Adll ^[35] by

The WLSEs, $\tilde{\theta}$ and $\tilde{\sigma},$ are obtained by minimizing the function

subject to the constrains $\sigma > 0$ and $\theta-a > 0,$ for some real constant $a.$ We get the MWLSEs, $\theta^*$ and $\sigma^*$ by minimizing the function,

subject to $\sigma > 0$ and $\theta-a > 0,$ where $a$ is a suitable real number chosen, provided that it is less than the minimum of the data. Minimization can be accomplished by equating the first partial derivatives with respect to $\theta$ and $\sigma$ with zero and then solving the resulting equations, numerically. The results are presented in Tables 1 and 2.

Remark 4.1.

The main reason for choosing the two-parameter exponential distribution in the simulation study is not only its theoretical importance but also because it has explicit forms for the maximum likelihood estimator of its parameters, and consequently, it can be considered a reference in comparison with the proposed method.

4.1.2.   Three-parameter Weibull distribution

We consider the three-parameter Weibull distribution with CDF

where $\alpha, \; \sigma$ and $\theta$ are the shape, scale and location parameters, respectively. As in the exponential distribution, the MLEs, WLSEs and MWLSEs of the parameters are obtained numerically and shown in Tables 3 and 4.

4.1.3.   The two-parameter modified Kies exponential distribution

Al-Babtain et al. ^[36] introduced a new family of continuous probability distributions, which they called the "new modified Kies family." They discussed the two-parameter MKi-exponential (MKiExp) distribution as a special case in detail. In the same paper, they demonstrated the practical importance and flexibility of fitting several types of real data. The CDF of MKiExp distribution is given by

where $\alpha$ and $\sigma$ are the shape and scale parameters, respectively. Abd El-Raheem et al. ^[37] consider the estimation problem of multiple constant-stress tests for progressive type-Ⅱ censored MKiExp data with binomial removals.

The problem of determining numerical estimates in this work is reduced to a suitable minimization problem subject to the following constraints: $\alpha > 0$ and $\sigma > 0.$ The results are given in Tables 5 and 6.

4.2.   Prediction

According to Kaminsky and Rhodin ^[13], the predictive likelihood function (PLF) based on the first $r$ observed order statistics is given by

The above PL function was extended to the GOSs model in Barakat et al. ^[28] for fixed and random sample sizes.

4.2.1.   Prediction of a future order statistic from the Weibull distribution

In view of (4.3), the log of the PLF for the three-parameter Weibull distribution, Weibull($\theta, \sigma, \alpha$), is as follows:

$x_1 < x_2 < ... < x_r < x_s.$ The predictive WLS function takes the form

where $w_i = \cfrac{r (n-i+1)}{(n+2) (n+1)^2}, \; i = 1, 2, ..., n.$ Similarly, the modified predictive WLS function based on the cumulative hazard function is given by

For OOSs we have $\mu^{*}_{i:n} = \sum_{j = 1}^{i} (n-j+1)^{-1}$ and $w^{*}_i = \left[\sum_{j = 1}^{i} (n-j+1)^{-2} \right]^{-1}, \; i = 1, 2, ..., n.$

Remark 4.2.

By setting $\alpha = 1$, the above results are specialized to the two-parameter exponential distribution. The prediction problem of the two-parameter exponential distribution is discussed in details in Barakat et al. ^[38]. Moreover, $\theta = 0$ yields the two-parameter Weibull distribution.

4.2.2.   Predicting future observations from the modified Kies exponential distribution

As in the Weibull distribution, the prediction of a future order statistic from the two-parameter modified Kies exponential distribution can be accomplished by minimizing the following three functions:

and

The minimization problem is subject to the constraints $\sigma > 0, \; \alpha > 0$ and $x_s > x_r.$

Remark 4.3.

The RMSE and PMC are obtained numerically via a simulation, and all computations are performed through Mathematica 13.1.

In view of the simulation studies given above, the following comments are extracted:

1. In all cases, for the scale and shape parameters as well as the point predictors, the RMSEs decreased as $r$ increased.

2. In most cases, the RMSE of the MLEs is smaller than the RMSEs of both the WLSEs and MWLSEs of the parameters.

3. For the Exp$(\theta, \sigma),$ according to PMC shown in Table 2, the MLE is the best followed by the MWLSE, which is followed by the WLSE for estimating the location parameter $\theta,$ while for estimating the scale parameter $\sigma$, the MWLSE is the best followed by the WLSE, which is followed by the MLE.

4. For the Weibull distribution, from Table 3, it is noted that

5. $RMSE(\hat{\theta}) < RMSE(\tilde{\theta}) < RMSE(\theta^{*})$, $RMSE(\sigma^{*}) \leq RMSE(\tilde{\sigma}) < RMSE(\hat{\sigma})$ and $RMSE(\alpha^{*}) < RMSE(\hat{\alpha})$ and $RMSE(\tilde{\alpha}) < RMSE(\hat{\alpha}).$

6. In view of PMC, Table 4 reveals that the WLSE is the best for the location parameter, while the MWLSEs of the scale and shape parameters are the best whenever $n-r$ is small, but the WLSEs are the best whenever $n-r$ is large. Moreover, for complete samples (i.e., $r = n$), the MWLSEs of the location, scale and shape parameters are the best.

7. According to the results presented in Tables 5 and 6, there is no preference for one method over the others when estimating the MKiExp distribution parameters.

8. In most cases, it is noted that the MWLSP and WLSP are better than the MLP, according to both the RMSE and PMC (see Tables 7–10).

5.   Applications

5.1.   The aircraft windscreen failure times

The windscreen on a large aircraft is a complex piece of equipment comprised of several layers. Failures of these items typically involve damage to or delamination of the heating system's nonstructural outer ply. These failures do not result in damage to the aircraft but do require replacement. Data of this type is incomplete in that all failure times have not yet been observed and may include failures to date of a particular model or combination of models. Murthy et al. ^[39] reported failure and service times for a specific model windscreen in Table 16.11 on page 297. The data represents 84 observed failure times for a specific windscreen device. Al-Babtain et al. ^[36] have shown that the MKiExp is an appropriate model for fitting this data.

We assume that the first $r$ failure times have been observed, and we apply our prediction methods to predict the remaining failure times in two different scenarios. Three-point predictors, along with their related estimates and their relative errors, are obtained in the first scenario. Alternatively, to avoid complications in computations, in the second scenario, we first estimate the parameters based on the first $r$ observed failure times and then compute three different point predictors for the future $s$th failure. In order to assess the prediction results, we compute the relative error (RE) for each point predictor. Recall that the RE is defined by $RE = 100 \times \frac{|\tilde{X}_{s:n}-X_{s:n}|}{X_{s:n}}$, where $\tilde{X}_{s:n}$ denotes the point predictor and $X_{s:n}$ is the exact value of the quantity to be predicted. The results are presented in Table 11.

5.2.   Application to reticulum cell sarcoma

The second data set is reported and analyzed by Hoel ^[40] and Abu El Azm ^[41] among others. According to these data, male mice received a radiation dose of 300 roentgen at an age of 5–6 weeks. Each mouse's cause of death was identified by autopsy as either thymic lymphoma, reticulum cell sarcoma or other causes. Reticulum cell sarcoma is designated as cause 1 in this instance, and the other two causes of death are merged to make cause 2. There were 77 observations in this set of data, 39 of which are attributable to the second cause of death, while 38 are attributed to the first. For analysis purposes, we consider the following observations that are due to the first cause of death: 317,318,399,495,525,536,549,552,554,557,558,571,586,594,596,605,612,621,628,631,636,643,647,648,649,661,663,666,670,695,697,700,705,712,713,738,748,753.

It has been shown that the Makeham-Gompertz distribution is quite adequate for fitting reticulum cell sarcomas (e.g., Hoel ^[40]). Arguing as in the first data set, the prediction results are shown in Table 12.

The real data analysis given the the two preceding applications revealed that the maximum RE appears for the WLS method, followed by the MWLS method, followed by the ML method.

6.   Conclusions

In this article, by using the cumulative hazard function, a new least squares method for estimation and prediction has been proposed. The method is presented in a general setup so that it can be applied to any model of GOSs. A simulation study and numerical comparisons based on the RMSEs and PMCs have been performed through three important probability distributions. For applicability, two examples of real data sets are provided to illustrate the prescribed method. The comparisons reveal that the method is comparable with the ML and WLS methods, in the sense that there is no obvious preference for one method over the others for all estimation and prediction situations. Moreover, analyzing the real data revealed that the second scenario, in which we first estimate the unknown distribution parameters using Type Ⅱ right censoring and then predict future unobserved failures, performs better than the first scenario for the three methods.

Use of AI tools declaration

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

Acknowledgments

The authors would like to express their sincere gratitude to the anonymous reviewers for their insightful comments and recommendations that raised the caliber of this paper.

This study is supported via funding from Prince Sattam Bin Abdulaziz University project number (PSAU/2023/R/1444).

Conflict of interest

The authors declare no conflict of interest.