Statistical inference for the Nadarajah-Haghighi distribution based on ranked set sampling with applications

1.   Introduction

McIntyre ^[12] developed a sampling strategy to ascertain the mean pasture yield in Australia, which was later known as ranked set sampling (RSS). This sampling technique is useful when it is more convenient and affordable to rank the sample units without utilizing their precise values. In order to generate RSS with set size $n$, we first draw $n$ simple random samples (SRSs) of size $n$ from the population of interest. Then, we arrange each sample of size $n$ in ascending order of magnitude. When ranking is done at this step, it is not required to actually measure the sample units because it can be done using a subjective method like eye inspection or personal judgement. The sample unit from the $i$th sample with judgement rank $i$, $i = 1, 2, ..., n$, has been chosen for actual quantification. It is important to note that a limited set size is suggested to facilitate informal rankings. RSS-based methods are frequently more effective than their SRS equivalents; this issue has been investigated for a number of common issues in the RSS literature, including, Dell and Clutter ^[2], Huang ^[5], Kvam and Samaniego ^[8], MacEachern et al. ^[11], Mohie El-Din et al. ^[14,15], Stokes ^[18], Stokes and Sager ^[19], Takahasi and Wakimoto ^[22], Wang et al. ^[23] and Zamanzade and Mahdizadeh ^[24].

The Nadarajah-Haghighi distribution (NHD) was developed by Nadarajah and Haghighi ^[16] as an alternative to the gamma, Weibull, and generalized-exponential distributions. Numerous NHD features are discussed by Nadarajah and Haghighi ^[16]. Let's assume that the lifetime of a testing unit $X$ follows a two-parameter NHD$(\theta, \mu)$; $\theta$ is the shape parameter and $\mu$ is the scale parameter. The probability density function (PDF), cumulative distribution function (CDF), survival function (SF) and hazard rate function (HRF) are all presented in accordance with the mission time $t$ by

Statistical inferences for the NHD have been studied in considerable detail. Singh et al. ^[20,21] have addressed maximum likelihood (MLEs) and Bayes estimations (BEs) based on progressively type Ⅱ censoring samples. Ashour et al. ^[1] calculated MLEs using the Newton-Raphson technique and BEs using the Markov chain Monte Carlo method (MCMC) based on progressive first-failure censoring samples. They employed highest posterior density intervals (HPD) and asymptotic confidence intervals (ACIs) to estimate intervals. Other sources, such as ^[6,7] have also called attention to the problems with an NHD's order statistics. They both worked on moment recurrence equations in order statistics.

Statistical inference for NHD unknown lifetime parameters have not yet been studied using RSS. This study's main objective was to obtain MLEs and BEs based on lifetime data collected under RSS to produce point and interval estimates of the unknown NHD parameters as well as some lifetime-parameters, such as the SF and HRF. Independent conjugate gamma priors of the unknown parameters are taken into account under the squared error loss function (SEL). Additionally, we construct ACIs and HPD for $\theta$ and $\mu$, as well as any function of them. In order to compute the BEs and produce the associated HPD, we recommend using Lindley-approximation and Metropolis-Hastings (MH) algorithm. A Monte-Carlo simulation (MC-simulation) analysis was carried out to assess the outcomes of various estimations through their mean-square error (MSE) and absolute-relative bias (RA). The following portions of the study are organized as follows: We derive the MLEs and associated two-sided ACIs in Section 2. In Section 3, we develop the SEL-based Bayesian inference. The results of MC-simulation are presented in Section 4. Section 5 describes inferences based on life-time data analysis as an illustration. Section 6 presents a final conclusion.

2.   MLE based on RSS

Consider $x_{j}$, $j = 1, 2, ..., n$, as an RSS from NHD($\theta, \mu$). Then, the likelihood-function (LF) based on $x_{j}$ is

When Eqs (1.1)–(1.3) are substituted into Eq (2.1), we obtain

where $\zeta(x_{j}, \mu) = (1+\mu x_{j})$ and $\psi(x_{j}, \theta, \mu) = 1-(1+\mu x_{j})^{\theta}$; the related log-LF, $\ell(.)\propto \log L(.)$, of Eq (2.2) may be expressed as

Differentiating Eq (2.3) partially with regard to $\theta$ and $\mu$, we find MLEs $\tilde{\theta}$, $\tilde{\mu}$, correspondingly, as

We have a system of two nonlinear formulae with unknown parameters, and it is clear from Eqs (2.4) and (2.5) that it is difficult to obtain a closed-form solution. Therefore, in order to acquire the required MLEs of the two undetermined parameters, we must use an appropriate iterative method, such as the Newton-Raphson approach. MLEs $\tilde{\kappa}(t)$ and $\tilde{\xi}(t)$ of $\kappa(t)$ and $\xi(t)$, as in Eqs (1.3) and (1.4), respectively, for a given time $t$, are further obtained via the invariance-property of MLEs $\tilde{\theta}$ and $\tilde{\mu}$. The asymptotic variance-covariance (VC) matrix of the MLEs $\tilde{\vartheta} = (\tilde{\theta}, \tilde{\mu})^{T}$ can be got by inverting the Fisher-information matrix, $I_{ij}(\vartheta) = E[-\partial^{2} \ell(\vartheta| data)/\partial \vartheta^{2}], \; i, j = 1, 2$, which is difficult to obtain by dropping the expectation operator $E$ and replacing $\vartheta$ by $\tilde{\vartheta}$. Consequently, we determine the MLEs' approximately VC-matrix via

as a result, the Fisher components will be

where

where

where

It is feasible to employ the asymptotic normality of MLEs $\tilde{\vartheta} \sim N(\tilde{\vartheta}, I^{-1}_{0}(\tilde{\vartheta}))$ to establish $100(1-a)%$ two-sided ACIs for the unknown parameters $\theta$ and $\mu$, consequently, as $\tilde{\theta} \pm Z_{a/2} \sqrt{Var(\tilde{\theta})}$ and $\tilde{\mu} \pm Z_{a/2} \sqrt{Var(\tilde{\mu})}$, where $Z_{a/2}$ represents the percentile of the standard normal distribution. According to Greene ^[3], we employ the delta approach to get approximations of the variances of $\tilde{\kappa}(t)$ and $\tilde{\xi}(t)$. As a result, the variance of $\tilde{\kappa}(t)$ and $\tilde{\xi}(t)$ may be roughly calculated by using $Var(\tilde{\kappa}(t)) = A^{t}I^{-1}_{0}(\tilde{\vartheta})A$ and $Var(\tilde{\xi}(t)) = B^{t}I^{-1}_{0}(\tilde{\vartheta})B$, where $A = \left[ \frac{{\partial} \kappa(t)}{\partial {\theta}}, \frac{{\partial} \kappa(t)}{\partial {\mu}}\right]_{(\theta = \tilde{\theta}, \; \mu = \tilde{\mu})}$ and $B = \left[ \frac{{\partial} \xi(t)}{\partial {\theta}}, \frac{{\partial} \xi(t)}{\partial {\mu}}\right]_{(\theta = \tilde{\theta}, \; \mu = \tilde{\mu})}.$ Thus, the $100(1-a)%$ two-sided ACIs of $\kappa(t)$ and $\xi(t)$, respectively, provided by $\tilde{\kappa}(t) \pm Z_{a/2} \sqrt{Var(\tilde{\kappa}(t))}$ and $\tilde{\xi}(t) \pm Z_{a/2} \sqrt{Var(\tilde{\xi}(t))}$.

3.   Bayesian inference

We use the assumption that the unknown parameters $\theta$ and $\mu$ are both independent random variables, and have conjugate gamma priors, i.e., $\theta \sim Gamma(c_{1}, d_{1})$ and $\mu \sim Gamma (c_{2}, d_{2})$, respectively, in order to derive the BEs. In order to reflect prior information regarding the two unknown parameters $\theta$ and $\mu$, the hyper parameters $c_{j}, d_{j}$, $j = 1, 2$ were chosen. As a result, $\pi(\theta, \mu) \propto \theta^{c_{1}-1} \mu^{c_{2}-1} e^{ -(d_{1}\theta +d_{2}\mu)}, \; \theta, \mu > 0\; \mathrm{and} \; c_{j}, d_{j} > 0, \; j = 1, 2$ can be applied to express the joint prior density of $\theta$ and $\mu$. The joint posterior PDF of $\theta$ and $\mu$ in the continuous Bayes' theorem can be defined by $\pi(\theta, \mu| {\bf{\underline{x}}}) \propto L(\theta, \mu | {\bf{\underline{x}}}) \pi(\theta, \mu)$. The joint posterior distribution is

As a consequence, the BEs of $\delta(\vartheta)$, $\vartheta = (\theta, \mu)$ under SEL, is the posterior expectation of $\delta(\vartheta)$ and is provided by

Nevertheless, the BEs are computed as two-dimensional integrals based on Eq (3.3), for which there is no closed-form solution. As a result, we advise adopting Lindley's approximation and MCMC methodology as two techniques for approximating Eq (3.3).

3.1.   Lindley-approximation based on RSS

Considering that the ratio form of posterior distribution contains an integration-denominator and couldn't be simplified to closed-form, so Lindley ^[10] introduced a way to assess the posterior expectation. Consequently, we approximate BEs using Lindley's technique. The posterior-expectation of $\phi(\alpha)$ is shown by

where $\phi(\alpha) = \ell(\alpha) +g(\alpha)$, $\ell(\alpha)$ represents log-LF and $g(\alpha) = \mathrm{log} (\pi(\alpha))$ represents log prior-function. Approximate BE, $\tilde{\phi}(\alpha_{1}, \alpha_{2})$, for the two-parameter case $\phi(\alpha_{1}, \alpha_{2})$ can be defined as

where

and $\sigma_{ij}$ represents the $(i, j)$th element of VC-matrix $I_{0}^{-1}(\tilde{\alpha_{1}}, \tilde{\alpha_{2}})$.

We can write Eq (3.4) as

where $\varphi_{1} = \sum^{2}_{i, j = 1} L_{ij1} \sigma_{ij}$ and $\varphi_{2} = \sum^{2}_{i, j = 1} L_{ij2} \sigma_{ij}$. Approximate BEs $\tilde{\theta}_{L}$ and $\tilde{\mu}_{L}$ of $\theta$ and $\mu$ based on SEL are supplied as

and

where $g_{1} = \frac{\partial{g} }{\partial{\theta}} = \frac{c_{1} -1}{\theta} -d_{1}$, $g_{2} = \frac{\partial{g} }{\partial{\mu}} = \frac{c_{2} -1}{\mu} -d_{2}$. So, to calculate the BEs using Lindley's approach, we need to have

where $D^{(\theta)}(x_{j}, \theta, \mu) = \frac{\partial D(x_{j}, \theta, \mu) }{\partial{\theta}}$.

where

The following equation provides an estimate of the BE $\tilde{\kappa}(t, \theta, \mu)$ of $\kappa(t, \theta, \mu)$:

where

Similarly, approximate BE of $\tilde{\xi}_{L}(t, \theta, \mu)$ of $\xi(t, \theta, \mu)$ is

where

However, we are unable to create the HPD credible intervals using Lindley's approximation. The Metropolis-Hastings (MH) method may be used to accomplish this goal by generating samples from the posterior distribution using MCMC, then computing BEs and their related HPD.

3.2.   MH-algorithm based on RSS

To implement MCMC approach and produce samples from the joint posterior distribution, we utilize MH sampler algorithm (Metropolis et al. ^[13] and Hastings ^[4]). According to Eq (3.1), the conditional posterior distributions of $\theta$ and $\mu$ are provided by

and

Analytical reduction to well-known distributions of the conditional posterior distributions of Eqs $(3.27)$ and $(3.28)$ of unknown parameters $\theta$ and $\mu$ is not possible. Hence, in order to get the BEs as well as accompanying HPD credible interval, MH algorithm with normal proposed distribution has been utilized to get random samples from such conditional posterior distributions. MH algorithm's stages are carried out as follows:

Step 1: Select initial values $\theta^{(0)}$ and $\mu^{(0)}$.

Step 2: Start with $J = 1$.

Step 3: Obtain $\theta ^{(J)}$ and $\mu ^{(J)}$ from Eqs $(3.27)$ and $(3.28)$ with normal proposal distributions $N(\theta^{(J-1)}, \sigma^{2}_{\theta})$ and $N(\mu^{(J-1)}, \sigma^{2}_{\mu})$ as

(a) Generate $\theta^{(*)}$ from $N(\theta^{(J-1)}, \sigma^{2}_{\theta})$, and $\mu^{(*)}$ from $N(\mu^{(J-1)}, \sigma^{2}_{\mu})$.

(b) Calculate the acceptance probability using

and

(c) Generate $u_{1}$ and $u_{2}$ from uniform distribution.

(d) If $u_{1}\le \phi_{\theta}$, don't reject proposal, then put $\theta^{(J)} = \theta^{*}$, otherwise put $\theta^{(J)} = \theta^{(J-1)}$.

(e) If $u_{2}\le \phi_{\mu}$, don't reject proposal, then put $\mu^{(J)} = \mu^{*}$, otherwise put $\mu^{(J)} = \mu^{(J-1)}$.

Step 4: At a specific time $t$, BEs of the SF and HRF are provided by

and

Step 5: Put $J = J+1$.

Step 6: Steps 3–5 should be repeated $M$-times to obtain

Step 7: To determine the HPD of $\rho = (\theta, \mu, \kappa(t), \xi(t))$, arrange the MCMC sample of $\rho^{(j)}$, $j = 1, 2, ..., M$ into $(\rho_{(1)}, \rho_{(2)}, ..., \rho_{(M)})$.

The $100(1-\gamma) %$ HPD of $\rho$ is given by $(\rho_{(J^{*})}, \rho_{(J^{*}+(1-\gamma) M)})$, where $J^{*}$ is selected such that

So, we can obtain the shortest length of HPD interval for $\rho$. The initial simulated iterations of the algorithm, $M_{0}$, are frequently discarded at the start of the analytical implementation (burn-in period), to guarantee convergence. The drawn samples, $\rho_{(j)}$, $j = M_{0}+1, ..., M$, are then large enough to construct the BEs. As a result, any function of $\theta$ and $\mu$ under SEL has MCMC BEs of $\rho(\theta, \mu)$ that are provided by

4.   MC-simulation analysis

This section compares the estimated values for $\theta$, $\mu$, $\kappa(t)$ and $\xi(t)$ using MC-simulation analysis. We calculated MLEs as well as BEs for the reliability characteristics and unknown-parameters. For Bayesian computing, Lindley's approximation as well as MCMC techniques using SEL were utilized. With regard to their MSE and RA values, various estimators have been compared. A $95%$ ACI/HPD intervals' performances have also been compared using average lengths (ALs). We determined our mean estimates, or MSEs in addition RAs, which are provided, respectively, by $\bar{\tilde{\vartheta}} = \frac{1}{\varsigma}\sum^{\varsigma}_{i = 1}\tilde{\vartheta}_{i}$, $\mathrm{MSE}(\tilde{\vartheta}) = \frac{1}{\varsigma}\sum^{\varsigma}_{i = 1}(\tilde{\vartheta}_{i}-\vartheta)^{2}$ and $\mathrm{RA}(\tilde{\vartheta}) = \frac{1}{\varsigma}\sum^{\varsigma}_{i = 1}\frac{|\tilde{\vartheta}_{i}-\vartheta|}{\vartheta}$, where $\varsigma = 1000$ and $\tilde{\vartheta}$ denotes the estimates of the parametric function. We have used the informative priors of $\theta$ and $\mu$ with $c_{1} = c_{2} = 1$ and $d_{1} = 0.2$, $d_{2} = 0.3$ when $(\theta, \mu) = (0.2, 0.3)$, such that the hyper-parameter values selected to satisfy the prior-mean are transformed into the predicted value of the associated parameter. For specified time $t = 5$, the true values of $\kappa(t)$ and $\xi(t)$ were further approximated to be $\kappa(5) = 0.8178$ and $\xi(5) = 0.0288$.

Based on true values of $\theta$ and $\mu$, we generated $n$ sets of random samples, each one of size $n$ from NHD$(\theta, \mu)$ distribution, by using the transformation $X_{i} = \frac{1}{\mu}\left(\left(1-\ln\left(1-U_{i}\right)\right)^{\frac{1}{\theta}}-1\right),$ $i = 1, \cdots, n$, $U_{i}$ was taken from $U(0, 1)$. Then, by using an RSS scheme (SCH-Ⅰ) and SRS scheme (SCH-Ⅱ), samples of size $n = 4, 7, 10$ were obtained.

The various BEs were produced by using MH sampler technique suggested in the theoretical section, with $12,000$ MCMC samples; we dropped the first $2000$ samples as burn-in. We calculated ALs of the $95%$ ACI/HPD of $\theta$, $\mu$, $\kappa(t)$ and $\xi(t)$. Statistical computer language Mathematica was used to carry out extensive computations. Table 1 displays mean of MLEs and BEs for $\theta$, $\mu$, $\kappa(t)$ and $\xi(t)$. Additionally, Table 2 includes the ALs of $\theta$, $\mu$, $\kappa(t)$ and $\xi(t)$.

Table 1 displays that the reliability features of the NHD in terms of the MSEs and RAs based on MLEs and BEs are extremely superb. MSEs and RAs of all estimates decline as $n$ increases, as predicted. Additionally, BEs utilizing gamma informative priors are superior to MLEs in terms of the MSEs and RAs because they involve prior knowledge. Based on the MSEs and RAs, MCMC technique employing the MH algorithm is superior to Lindley's approximation method. HPD outperform ACIs in regard to ALs for interval estimation. Additionally, by increasing $n$, the ALs of ACIs and HPD credible intervals decrease. We recommend applying point and interval BEs using the MH method. Furthermore, MSEs and RAs of MLEs and BEs for parameters were higher for SCH-Ⅱ than SCH-Ⅰ in most cases when comparing SCH-Ⅰ and SCH-Ⅱ. Given that the experiment's observations for SCH-Ⅰ were more random than for SCH-Ⅱ, therefore, it would be expected that the data collected by SCH-Ⅰ would reveal more details regarding reliability-parameters than the sample collected by SCH-Ⅱ.

5.   Real-life data analysis

In order to show how the suggested technique may be used to study real-life events, a real data set was examined. The analyzed data-set, which depicts failure times for 33 small electronic applications with known failure causes in an autonomous life test, is shown in ^[9]. There were totals of 17 and 16 observed failures attributed to causes 1 and 2, respectively, from the entire failure times. For computational ease, each failure time point in the original data set was split by a thousand. Salem et al. ^[17] analyzed and addressed transformed failure times of electronic applications. Before sampling, we fitted the NHD to the whole set of data and compared it with the fitting of three-lifetime distributions, specifically, the gamma distribution (G-D), Weibull distribution (W-D) and generalized-exponential distribution (GE-D).

The data's compatibility with the G-D, W-D, GE-D and NHD is one issue that has to be addressed. The reliability of the proposed model has been evaluated using the goodness of fit test (GFT) statistics with related pvalues from the Kolmogorov Smirnov (KS) and Anderson Darling (AD) tests. Additionally, we evaluated the following:

● negative-log-LF$= (-\widetilde{l})$ at MLEs,

● Akaike-criteria, $AC = -2\widetilde{l}+2w$,

● Akaike-criteria-correction, $ACc = AC+\frac{2nw}{n-w-1}$,

● Akaike-criteria-consistent, $CAC = -2\widetilde{l}+\frac{2nw}{n-w-1}$,

● Bayesian-criterion, $BC = -2\widetilde{l}+w\log (n)$,

● Hannan-Quinn-criteria, $HQC = -2\widetilde{l}+2w\log(\log(n))$, $n$ is sample size and $w$ is the number of model parameters are taken into account, see Ashour et al. ^[1].

The optimal distribution often has the greatest pvalues and lowest values for the $-\widetilde{l}$, AC, ACc, CAC, BC, HQC, KS, and AD statistics. Table 3 lists the values for MLEs of the distribution parameters and the accompanying GFT metrics. Table 4 provides the statistical results of the KS and AD GFTs, and their pvalues. Additionally, we assessed the distributions' goodness of fit by using a graphing method. For the G-D, W-D, GE-D and NHD, we created quantile-quantile (Q-Q) plots, which are shown in Figure 1. The points $\left(F^{-1}\left(\frac{(i-0.5)}{n}; \widetilde{\theta}\right), x_{(i)}\right)$, $i = 1, 2, ..., n$ are shown on a Q-Q plot; $\widetilde{\theta}$ is the MLE of $\theta$ and $x_{(i)}$ is ordered data. We have also included two graphs that were produced using the estimated parameters for a more precise explanation. The SFs of the G-D, W-D, GE-D and NHD have been fitted by using the empirical CDF plot, which is the first plot. As shown in Figure 2, the subsequent plot represents the histogram of the real-data together with PDFs for the G-D, W-D, GE-D and NHD. The fact that it is frequently impossible to demonstrate the existence and uniqueness of the generated estimators is one of the main problems with ML inference. The contour plot of the log-LF based on NHD parameters is presented in Figure 3 to help address this issue using the observed failure data. It shows the existence and uniqueness of the MLEs $\widetilde{\theta}$ and $\widetilde{\mu}$. Therefore, to begin the computational iteration, we suggest using these estimates as initial hypotheses.

The NHD is the best option among the comparable G-D, W-D and GE-D that have been studied to fit lifetime data as it had the lowest GFTs and greatest pvalues, as can be observed from Tables 3 and 4. The Q-Q plots also lend credence to the aforementioned conclusions.

We can obtain at random six sets with $n = 6$ items in each set to develop RSS to demonstrate the inferential techniques used in the previous sections (see Table 5). The reliable features $\kappa(t)$ and $\xi(t)$ at the specified time $t = 2$, in addition to MLEs and BEs of $\theta$ and $\mu$, were calculated and are listed in Table 6. Under a non-informative prior, Lindley-approximation and MCMC method were used to develop BEs. The first 10,000 iterations were removed from the sequence after we generated 50,000 MCMC samples. The MLEs for $\theta$, $\mu$, $\kappa(t)$ and $\xi(t)$ were used as initial-values for the MCMC sampler procedure. Additionally, Table 6 lists a $95%$ ACIs and HPD intervals. Forty-thousand outcome of $\theta$, $\mu$, $\kappa(t)$ and $\xi(t)$ are shown in Figure 4 together with their sample means (solid line (—)) and $95%$ two sided intervals (dashed lines (- - -)). Additionally, this plot demonstrates that MCMC algorithm has achieved good convergence and burn-in-period sample size is appropriate to eliminate initial proposals effects.

6.   Conclusions

The issue of determining the NHD's life-time parameters when the observed data are drawn by RSS is addressed in this article based on MLEs and BEs. These estimations can be computed numerically but cannot be achieved in closed forms. The BEs were got by utilizing the Lindley-approximation and MCMC techniques with gamma-priors under SEL. Asymptotic-normality MLEs as well as the delta technique were used to construct $95%$-ACIs of life-time parameters. We used MH technique to get point estimations and related HPD credible intervals because Lindley approximation approach is unable to develop the HPD credible intervals.

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 extend their appreciation to the Deputyship for Research & Innovation, Ministry of Education in Saudi Arabia for funding this research work through the project (No. IFKSUOR3–216–1).

Conflict of interest

The authors declare that they have no conflicts of interest.