Bayesian and Classical Estimation of Stress-Strength Reliability for Inverse Weibull Lifetime Models

Abstract

In this paper, we consider the problem of estimating stress-strength reliability for inverse Weibull lifetime models having the same shape parameters but different scale parameters. We obtain the maximum likelihood estimator and its asymptotic distribution. Since the classical estimator doesn’t hold explicit forms, we propose an approximate maximum likelihood estimator. The asymptotic confidence interval and two bootstrap intervals are obtained. Using the Gibbs sampling technique, Bayesian estimator and the corresponding credible interval are obtained. The Metropolis-Hastings algorithm is used to generate random variates. Monte Carlo simulations are conducted to compare the proposed methods. Analysis of a real dataset is performed.

1. Introduction

The inference of stress-strength reliability in statistics is an important topic of interest. It has many applications in practical areas. In the stress-strength modeling, $R=P(Y<X)$ is a measure of component reliability. If $X\le Y$, the component fails or the component may malfunction, where X is subject to Y. R can also be considered in electrical and electronic systems. Many authors have studied its properties for many statistical models including double exponential, Weibull, generalized Pareto and Lindley distribution (see [1,2,3,4]). Classical and Bayesian estimation of reliability in a multicomponent stress-strength model under a general class of inverse exponentiated distributions was researched by Ref. [5]. Ref. [6] studied classical and Bayesian estimation of reliability in multicomponent stress-strength model under Weibull distribution. Furthermore, Refs. [7,8,9] also considered the problem of the stress-strength reliability.

The inverse Weibull (IW) distribution has attracted much attention recently. If T denotes the random variable from Weibull model, we define X as follows $X=1/T$, then the random variable X is said to have inverse Weibull distribution. It is a lifetime probability distribution that can be used in the reliability engineering discipline. The inverse Weibull distribution has the ability to model failure rates, which is quite common in reliability and biological studies. The inverse Weibull model was referred to with many different names like “Frechet-type” ([10]) and “Complementary Weibull” ([11]). Ref. [11] discussed a graphical plotting technique to settle the suitability of the model. Ref. [12] presented the IW distribution for modeling reliability data, this model was further discussed by researching the failures of mechanical components subject to degradation. Ref. [13] proposed a discrete inverse Weibull distribution and its parameters were estimated. The mixture model of two IW distributions and its identifiability properties were studied by [14]. For the theoretical analysis of IW distribution, we can refer to [15]. Ref. [16] proposed the generalized IW distribution and several properties of this model. For more details on the inverse Weibull distribution, see [17].

In this paper, we focus on the estimation of the stress-strength reliability $R=P(Y<X)$, where X and Y follow the inverse Weibull distribution. As far as we know, this model has not been previously studied, although, we believe it plays an important role in reliability analysis.

We obtain the maximum likelihood estimator (MLE), approximate maximum likelihood estimator (AMLE) and the asymptotic distribution of the estimator. The asymptotic distribution is used to construct an asymptotic confidence interval. We also present two bootstrap confidence intervals of R. By using Gibbs sampling technique, we obtain Bayes estimator of R and its corresponding credible interval. Finally, we present a real data example to illustrate the performance of different methods.

The layout of this paper is organized as follows: in Section 2, we introduce the distribution of the inverse Weibull. In Section 3, we obtain the MLE of R. In Section 4, we derive the estimator of R by approximating maximum likelihood equations. Different confidence intervals are presented in Section 5. In Section 6, Bayesian solutions are introduced. In Section 7, we compare different proposed methods using Monte Carlo simulation. A numerical example is also provided. Finally, in Section 8, we conclude the paper.

2. Inverse Weibull Distribution

The probability density function of the known Weibull distribution is given by

$$\begin{array}{c}\hfill f(t;\alpha ,\theta )=\frac{\alpha }{\theta }{t}^{\alpha -1}{e}^{-\frac{t^{\alpha }}{\theta }},t>0,\end{array}$$

(1)

where $\alpha >0$ is the shape parameter and $\theta >0$ is the scale parameter.

Let T denote the random variable from Weibull model, namely, $W(\alpha ,\theta )$. Define X as follows:

$$\begin{array}{c}\hfill X=\frac{1}{T}.\end{array}$$

(2)

The random variable X is said to have inverse Weibull distribution, and its probability density function (pdf) is given by

$$\begin{array}{c}\hfill f(x;\alpha ,\theta )=\frac{\alpha }{\theta }{x}^{-\alpha -1}{e}^{-\frac{x^{-\alpha }}{\theta }},x>0.\end{array}$$

(3)

The cumulative distribution function(cdf) is given by

$$\begin{array}{c}\hfill F(x)=e^{-\frac{x^{-\alpha }}{\theta }},x>0,\end{array}$$

(4)

where $\alpha >0$ and $\theta >0$. The inverse Weibull distribution will be denoted by $IW(\alpha ,\theta )$.

3. Maximum Likelihood Estimator of $\mathit{R}$

In this section, we consider the problem of estimating $R=P(Y<X)$ under the assumption that $X\sim IW(\alpha ,{\theta }_1)$ and $Y\sim IW(\alpha ,{\theta }_2)$. Then, it can be easily calculated that

$$\begin{array}{c}\hfill R=P(Y<X)=\frac{{\theta }_2}{{\theta }_1+{\theta }_2}.\end{array}$$

(5)

To computer the MLE of R, first we obtain the MLEs of ${\theta }_1$ and ${\theta }_2$. Suppose $X_1,X_2,...,X_n$ is a random sample from $IW(\alpha ,{\theta }_1)$ and $Y_1,Y_2,...,Y_m$ is a random sample from $IW(\alpha ,{\theta }_2)$. The joint likelihood function is:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill l(\alpha ,{\theta }_1,{\theta }_2)=\frac{{\alpha }^{n+m}}{{\theta }_1^n{\theta }_2^m}(\prod _{i=1}^n{x}_i^{-\alpha -1})(\prod _{j=1}^m{y}_j^{-\alpha -1})e^{-({\sum }_{i=1}^n\frac{x_i^{-}\alpha }{{\theta }_1}+{\sum }_{j=1}^m\frac{y_j^{-}\alpha }{{\theta }_2})}.\end{array}\end{array}$$

Then, the log-likelihood function is

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill L(\alpha ,{\theta }_1,{\theta }_2)=& (m+n)\ln \alpha -n\ln {\theta }_1-m\ln {\theta }_2-(\alpha +1)\left[\sum _{i=1}^n\ln x_i+\sum _{j=1}^m\ln y_j\right]\hfill \\ & -\frac{1}{{\theta }_1}\sum _{i=1}^n{x}_i^{-\alpha }-\frac{1}{{\theta }_2}\sum _{j=1}^m{y}_j^{-\alpha }.\hfill \end{array}\end{array}$$

(6)

$\widehat{\alpha }$, $\widehat{{\theta }_1}$ and $\widehat{{\theta }_2}$, the MLEs of the parameters $\alpha$, ${\theta }_1$ and ${\theta }_2$, can be numerically obtained by solving the following equations:

$$\begin{array}{c}\frac{\partial L}{\partial \alpha }=\frac{m+n}{\alpha }-\sum _{i=1}^n\ln x_i-\sum _{j=1}^m\ln y_j+\frac{1}{{\theta }_1}\sum _{i=1}^n{x}_i^{-\alpha }\ln x_i+\frac{1}{{\theta }_2}\sum _{j=1}^m{y}_j^{-\alpha }\ln y_j=0,\hfill \end{array}$$

(7)

$$\begin{array}{c}\frac{\partial L}{\partial {\theta }_1}=-\frac{n}{{\theta }_1}+\frac{1}{{\theta }_1^2}\sum _{i=1}^n{x}_i^{-\alpha }=0,\hfill \end{array}$$

(8)

$$\begin{array}{c}\frac{\partial L}{\partial {\theta }_2}=-\frac{m}{{\theta }_2}+\frac{1}{{\theta }_2^2}\sum _{j=1}^m{y}_j^{-\alpha }=0.\hfill \end{array}$$

(9)

From (8) and (9), we obtain

$$\begin{array}{c}\hfill \widehat{{\theta }_1}(\alpha )=\frac{1}{n}\sum _{i=1}^n{x}_i^{-\alpha }\mathrm{and}\widehat{{\theta }_2}(\alpha )=\frac{1}{m}\sum _{j=1}^m{y}_j^{-\alpha }.\end{array}$$

(10)

Putting the expressions of $\widehat{{\theta }_1}(\alpha )$ and $\widehat{{\theta }_2}(\alpha )$ into (7), we obtain

$$\begin{array}{ccc}\hfill \frac{m+n}{\alpha }-\frac{1}{\alpha }\left[\sum _{i=1}^n\ln x_i^{\alpha }+\sum _{j=1}^m\ln x_j^{\alpha }\right]+\frac{{\sum }_{i=1}^n{x}_i^{-\alpha }\ln x_i}{\frac{1}{n}{\sum }_{i=1}^n{x}_i^{-\alpha }}+\frac{{\sum }_{j=1}^m{y}_j^{-\alpha }\ln y_j}{\frac{1}{m}{\sum }_{j=1}^m{y}_j^{-\alpha }}& =& 0.\hfill \end{array}$$

(11)

Therefore, $\widehat{\alpha }$ can be obtained as a fixed point solution of the non-linear equation of the form

$$\begin{array}{ccc}\hfill h(\alpha )& =& \alpha ,\hfill \end{array}$$

(12)

where

$$\begin{array}{ccc}\hfill h(\alpha )& =& \frac{-m-n+{\sum }_{i=1}^n\ln x_i^{\alpha }+{\sum }_{j=1}^m\ln y_j^{\alpha }}{\frac{{\sum }_{i=1}^n{x}_i^{-\alpha }\ln x_i}{\frac{1}{n}{\sum }_{i=1}^n{x}_i^{-\alpha }}+\frac{{\sum }_{j=1}^m{y}_j^{-\alpha }\ln y_j}{\frac{1}{m}{\sum }_{j=1}^m{y}_j^{-\alpha }}}.\hfill \end{array}$$

(13)

Using a simple iterative procedure $h({\alpha }_{(j)})={\alpha }_{(j+1)}$, where ${\alpha }_{(j)}$ is the j-th iterate of $\widehat{\alpha }$, we stop the iterative procedure when $|{\alpha }_{(j)}-{\alpha }_{(j+1)}|$ is adequately small less than a specified level. Once we obtain $\widehat{\alpha }$, then $\widehat{{\theta }_1}$ and $\widehat{{\theta }_2}$ can be calculated from (10). Therefore, we obtain the MLE of $R=P(Y<X)$ as

$$\begin{array}{ccc}\hfill \widehat{R}& =& \frac{\frac{1}{m}{\sum }_{j=1}^m{y}_j^{-\widehat{\alpha }}}{\frac{1}{n}{\sum }_{i=1}^n{x}_i^{-\widehat{\alpha }}+\frac{1}{m}{\sum }_{j=1}^m{y}_j^{-\widehat{\alpha }}}.\hfill \end{array}$$

(14)

4. Approximate Maximum Likelihood Estimator of $\mathit{R}$

The MLEs do not take explicit forms; therefore, we approximate the likelihood equation and derive explicit estimators of the parameters.

Since the random variable X follows $IW(\alpha ,\theta )$, then $V=\ln X$ has the extreme value distribution with pdf as

$$\begin{array}{c}\hfill f(v;\mu ,\sigma )=-\frac{1}{\sigma }{e}^{\frac{v-\mu }{\sigma }-e^{\frac{v-\mu }{\sigma }}},-\infty <v<+\infty ,\end{array}$$

(15)

where $\mu =-\frac{1}{\alpha }\ln \theta$ and $\sigma =-\frac{1}{\alpha }$. The $\mu$ and $\sigma$ are location and scale parameters, respectively. The pdf and cdf of the standard extreme value distribution can be obtained as

$$\begin{array}{c}\hfill g(v)=e^{-v-e^{-v}},G(v)=e^{-e^{-v}}.\end{array}$$

(16)

Suppose $X_{(1)}<X_{(2)}<...<X_{(n)}$ and $Y_{(1)}<Y_{(2)}<...<Y_{(m)}$ are the ordered $X_{i}s$ and $Y_{j}s$, we assume the following notations: $T_{(i)}=\ln X_{(i)}$, $Z_{(i)}=\frac{T_{(i)}-{\mu }_1}{-\sigma },i=1,...,n$ and $S_{(j)}=\ln Y_{(j)}$, $W_{(j)}=\frac{S_{(j)}-{\mu }_2}{-\sigma },j=1,...,m$, where ${\mu }_1=-\frac{1}{\alpha }\ln {\theta }_1$, ${\mu }_2=-\frac{1}{\alpha }\ln {\theta }_2$ and $\sigma =-\frac{1}{\alpha }$.

The log-likelihood function of the data $T_{(1)},...,T_{(n)}$ and $S_{(1)},...,S_{(m)}$ is

$$\begin{array}{c}\hfill L({\mu }_1,{\mu }_2,\sigma )\propto -(m+n)\ln (-\sigma )+\sum _{i=1}^n\ln g(z_{(i)})+\sum _{j=1}^m\ln g(w_{(j)}).\end{array}$$

(17)

Differentiating (17) in regard to ${\mu }_1,{\mu }_{2}and\sigma$, the score equations are obtained as

$$\begin{array}{c}\frac{\partial L}{\partial {\mu }_1}=\frac{1}{\sigma }\sum _{i=1}^n\frac{g^{\prime }(z_{(i)})}{g(z_{(i)})}=0,\hfill \end{array}$$

(18)

$$\begin{array}{c}\frac{\partial L}{\partial {\mu }_2}=\frac{1}{\sigma }\sum _{j=1}^m\frac{g^{\prime }(w_{(j)})}{g(w_{(j)})}=0,\hfill \end{array}$$

(19)

$$\begin{array}{c}\frac{\partial L}{\partial \sigma }=-(m+n)\frac{1}{\sigma }-\frac{1}{\sigma }\sum _{i=1}^n\frac{g^{\prime }(z_{(i)})}{g(z_{(i)})}{z}_{(i)}-\frac{1}{\sigma }\sum _{j=1}^m\frac{g^{\prime }(w_{(j)})}{g(w_{(j)})}{w}_{(j)}=0.\hfill \end{array}$$

(20)

We note that the function $h(z_{(i)})=\frac{g^{\prime }(z_{(i)})}{g(z_{(i)})}$ makes the score Equation (18) nonlinear and intricate. Thus, we approximate the function $h(z_{(i)})=\frac{g^{\prime }(z_{(i)})}{g(z_{(i)})}$ by expanding it in a Taylor series around $c_i=E(Z_{(i)})$. Furthermore, we also approximate the function $h(w_{(j)})=\frac{g^{\prime }(w_{(j)})}{g(w_{(j)})}$ by expanding it in a Taylor series around $d_j=E(W_{(j)})$. From [18], it is known that

$$\begin{array}{c}\hfill G(Z_{(i)})\stackrel{d}{=}{U}_{(i)},\end{array}$$

where $U_{(i)}$ is the i-th order statistic from the uniform $U(0,1)$ distribution. Therefore,

$$\begin{array}{c}\hfill Z_{(i)}\stackrel{d}{=}{G}^{-1}(U_{(i)}),\end{array}$$

and

$$\begin{array}{c}\hfill c_i=E{Z}_{(i)}\approx G^{-1}(E{U}_{(i)})=G^{-1}(i/(n+1)).\end{array}$$

We use the following notations, $p_i=\frac{i}{n+1}$, ${\overline{p}}_j=\frac{j}{m+1}$; therefore, $c_i=G^{-1}(p_i)=-\ln (-\ln p_i)$, $d_j=G^{-1}({\overline{p}}_j)=-\ln (-\ln {\overline{p}}_j).$

Expanding the function $h(z_{(i)})$ and $h(w_{(j)})$ and keeping the first two terms, we have

$$\begin{array}{c}h(z_{(i)})=\frac{g^{\prime }(z_{(i)})}{g(z_{(i)})}\approx h(c_i)+h^{\prime }(c_i)(z_{(i)}-c_i)\equiv a_i-b_i{z}_{(i)},i=1,...,n,\hfill \\ h(w_{(j)})=\frac{g^{\prime }(w_{(j)})}{g(w_{(j)})}\approx h(d_j)+h^{\prime }(d_j)(w_{(j)}-d_j)\equiv {\overline{a}}_j-{\overline{b}}_j{z}_{(j)},j=1,...,m,\hfill \end{array}$$

where

$$\begin{array}{c}{a}_i=h(c_i)-c_i{h}^{\prime }(c_i)=\ln p_i(\ln (-\ln p_i)-1)-1,\hfill \\ b_i=-h^{\prime }(c_i)=-\ln p_i,\hfill \\ {\overline{a}}_j=h(d_j)-d_j{h}^{\prime }(d_j)=\ln {\overline{p}}_j(\ln (-\ln {\overline{p}}_j)-1)-1,\hfill \\ {\overline{b}}_j=-h^{\prime }(d_j)=-\ln {\overline{p}}_j.\hfill \end{array}$$

Therefore, (18)–(20) can be represented as

$$\begin{array}{c}\frac{\partial L}{\partial {\mu }_1}\approx \frac{1}{\sigma }\sum _{i=1}^n(a_i-b_i{z}_{(i)})=0,\hfill \end{array}$$

(21)

$$\begin{array}{c}\frac{\partial L}{\partial {\mu }_2}\approx \frac{1}{\sigma }\sum _{j=1}^m({\overline{a}}_j-{\overline{b}}_j{w}_{(j)})=0,\hfill \end{array}$$

(22)

$$\begin{array}{c}\frac{\partial L}{\partial \sigma }\approx -\frac{1}{\sigma }\left(m+n+\sum _{i=1}^n(a_i-b_i{z}_{(i)})z_{(i)}+\sum _{j=1}^m({\overline{a}}_j-{\overline{b}}_j{w}_{(j)})w_{(j)}\right)=0.\hfill \end{array}$$

(23)

The estimator of ${\mu }_1$ and ${\mu }_2$ can be obtained from Equations (21) and (22), as

$$\begin{array}{c}{\tilde{\mu }}_1=A_1+B_1\tilde{\sigma },\hfill \end{array}$$

(24)

$$\begin{array}{c}{\tilde{\mu }}_2=A_2+B_2\tilde{\sigma },\hfill \end{array}$$

(25)

where

$$\begin{array}{c}\hfill A_1=\frac{{\sum }_{i=1}^n{b}_i{T}_{(i)}}{{\sum }_{i=1}^n{b}_i},A_2=\frac{{\sum }_{j=1}^m{\overline{b}}_j{S}_{(j)}}{{\sum }_{j=1}^m{\overline{b}}_j},B_1=\frac{{\sum }_{i=1}^n{a}_i}{{\sum }_{i=1}^n{b}_i},B_2=\frac{{\sum }_{j=1}^m{\overline{a}}_j}{{\sum }_{j=1}^m{\overline{b}}_j}.\end{array}$$

(26)

The estimator of $\sigma <0$ can be determined as the unique solution of the quadratic equation

$$\begin{array}{c}\hfill C{\sigma }^2+D\sigma -E=0,\end{array}$$

(27)

where

$$\begin{array}{c}C=(m+n)+B_1\sum _{i=1}^n{a}_i+B_2\sum _{j=1}^m{\overline{a}}_j-B_1^2\sum _{i=1}^n{b}_i-B_2^2\sum _{j=1}^m{\overline{b}}_j=m+n,\hfill \\ D=\sum _{i=1}^n{a}_i(A_1-T_{(i)})-2B_1\sum _{i=1}^n{b}_i(A_1-T_{(i)})+\sum _{j=1}^m{\overline{a}}_j(A_2-S_{(j)})-2B_2\sum _{j=1}^m{\overline{b}}_j(A_2-S_{(j)}),\hfill \\ E=\sum _{i=1}^n{b}_i{(T_{(i)}-A_1)}^2+\sum _{j=1}^m{\overline{b}}_j{(S_{(j)}-A_2)}^2>0.\hfill \end{array}$$

Therefore,

$$\begin{array}{c}\hfill \tilde{\sigma }=\frac{-D-\sqrt{D^2+4E(m+n)}}{2(m+n)},\end{array}$$

(28)

Once $\tilde{\sigma }$ is obtained, $\widetilde{{\mu }_1}$ and $\widetilde{{\mu }_2}$ can be derived immediately. Hence, the AMLE of R is given by

$$\begin{array}{c}\hfill \tilde{R}=\frac{{\tilde{\theta }}_2}{{\tilde{\theta }}_1+{\tilde{\theta }}_2},\end{array}$$

(29)

where

$$\begin{array}{c}\hfill \tilde{\alpha }=-\frac{1}{\tilde{\sigma }},\widetilde{{\theta }_1}=e^{\frac{1}{\tilde{\sigma }}(A_1+B_1\tilde{\sigma })},\widetilde{{\theta }_2}=e^{\frac{1}{\tilde{\sigma }}(A_2+B_2\tilde{\sigma })}.\end{array}$$

5. Confidence Intervals of R

In this section, we present an asymptotic confidence interval (C.I.) of R and two C.I.s based on the non-parametric bootstrap methods.

5.1. Asymptotic Confidence Interval of R

In this subsection, we derive the asymptotic distribution of the MLE $\widehat{\theta }=(\widehat{\alpha },\widehat{{\theta }_1},\widehat{{\theta }_2})$ and $\widehat{R}$. Based on the asymptotic distribution of $\widehat{R}$, the corresponding asymptotic confidence interval of R is obtained. We denote the exact Fisher information matrix of $\theta =(\alpha ,{\theta }_1,{\theta }_2)$ as $J(\theta )=-\mathbb{E}(I;\theta )$, where $I={\left[I_{ij}\right]}_{i,j=1,2,3}$, $I_{ij}={\partial }^{2}L/\partial {\theta }_i\partial {\theta }_j$ and L is given in (6):

$$\begin{array}{c}\hfill I(\theta )=\left(\begin{array}{ccc}\frac{{\partial }^{2}L}{\partial {\alpha }^2}& \frac{{\partial }^{2}L}{\partial \alpha \partial {\theta }_1}& \frac{{\partial }^{2}L}{\partial \alpha \partial {\theta }_2}\\ \frac{{\partial }^{2}L}{\partial {\theta }_1\partial \alpha }& \frac{{\partial }^{2}L}{\partial {\theta }_1^2}& \frac{{\partial }^{2}L}{\partial {\theta }_1\partial {\theta }_2}\\ \frac{{\partial }^{2}L}{\partial {\theta }_2\partial \alpha }& \frac{{\partial }^{2}L}{\partial {\theta }_2\partial {\theta }_1}& \frac{{\partial }^{2}L}{\partial {\theta }_2^2}\end{array}\right)=\left(\begin{array}{ccc}{I}_{11}& I_{12}& I_{13}\\ I_{21}& I_{22}& I_{23}\\ I_{31}& I_{32}& I_{33}\end{array}\right).\end{array}$$

It is easy to see that

$$\begin{array}{c}{I}_{11}=-\frac{m+n}{{\alpha }^2}-\frac{1}{{\theta }_1}\sum _{i=1}^n{(\ln x_i)}^2{x}_i^{-\alpha }-\frac{1}{{\theta }_2}\sum _{j=1}^m{(\ln y_j)}^2{y}_j^{-\alpha },\hfill \\ I_{12}=I_{21}=-\frac{1}{{\theta }_1^2}\sum _{i=1}^n{x}_i^{-\alpha }\ln x_i,\hfill \\ I_{13}=I_{31}=-\frac{1}{{\theta }_2^2}\sum _{j=1}^m{y}_j^{-\alpha }\ln y_j,\hfill \\ I_{22}=\frac{n}{{\theta }_1^2}-\frac{2}{{\theta }_1^3}\sum _{i=1}^n{x}_i^{-\alpha },\hfill \\ I_{33}=\frac{m}{{\theta }_2^2}-\frac{2}{{\theta }_2^3}\sum _{j=1}^m{y}_j^{-\alpha },\hfill \\ I_{23}=I_{32}=0.\hfill \end{array}$$

Moreover,

$$\begin{array}{c}{J}_{11}=-\mathbb{E}(\frac{{\partial }^{2}L}{\partial {\alpha }^2})=\frac{1}{{\alpha }^2}\left[(m+n)(1+{\gamma }^{″}(2))+n{(\ln {\theta }_1)}^2+m{(\ln {\theta }_2)}^2+2{\gamma }^{\prime }(2)(n\ln {\theta }_1+m\ln {\theta }_2)\right],\hfill \\ J_{12}=J_{21}=\frac{1}{{\theta }_1^2}\sum _{i=1}^n\mathbb{E}(x_i^{-\alpha }\ln x_i)=-\frac{n}{{\theta }_1\alpha }\left[\ln {\theta }_1+{\gamma }^{\prime }(2)\right],\hfill \\ J_{13}=J_{31}=\frac{1}{{\theta }_2^2}\sum _{j=1}^m\mathbb{E}(y_j^{-\alpha }\ln y_j)=-\frac{m}{{\theta }_2\alpha }\left[\ln {\theta }_2+{\gamma }^{\prime }(2)\right],\hfill \\ J_{22}=\frac{n}{{\theta }_1^2},\hfill \\ J_{33}=\frac{m}{{\theta }_2^2},\hfill \\ J_{23}=J_{32}=0,\hfill \end{array}$$

where $\gamma (\alpha )={\int }_0^{\infty }{x}^{\alpha -1}{e}^{-x}dx$ is the Gamma function.

Theorem 1.

As $n\to \infty$ and $m\to \infty$ and $\frac{n}{m}\to p$, then

$$\begin{array}{c}\left(\sqrt{m}(\widehat{\alpha }-\alpha ),\sqrt{n}(\widehat{{\theta }_1}-{\theta }_1),\sqrt{m}(\widehat{{\theta }_2}-{\theta }_2)\right)\stackrel{d}{\to }{N}_3\left(0,A^{-1}(\alpha ,{\theta }_1,{\theta }_2)\right),\hfill \end{array}$$

where

$$\begin{array}{c}\hfill A(\alpha ,{\theta }_1,{\theta }_2)=\left(\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& 0\\ a_{31}& 0& a_{33}\end{array}\right),\end{array}$$

and

$$\begin{array}{c}{a}_{11}=\frac{1}{{\alpha }^2}\left[(1+p)(1+{\gamma }^{\prime \prime }(2))+p{(\ln {\theta }_1)}^2+{(\ln {\theta }_2)}^2+2{\gamma }^{\prime }(2)(p\ln {\theta }_1+\ln {\theta }_2)\right]=\underset{n,m\to \infty }{\lim }\frac{J_{11}}{m},\hfill \\ a_{12}=a_{21}=-\frac{\sqrt{p}}{{\theta }_1\alpha }\left[ln{\theta }_1+{\gamma }^{\prime }(2)\right]=\underset{n,m\to \infty }{\lim }\frac{\sqrt{p}}{n}{J}_{12},\hfill \\ a_{13}=a_{31}=-\frac{1}{{\theta }_2\alpha }\left[ln{\theta }_2+{\gamma }^{\prime }(2)\right]=\underset{n,m\to \infty }{\lim }\frac{1}{m}{J}_{13},\hfill \\ a_{22}=\frac{1}{{\theta }_1^2}=\underset{n,m\to \infty }{\lim }\frac{1}{n}{J}_{22},\hfill \\ a_{33}=\frac{1}{{\theta }_2^2}=\underset{n,m\to \infty }{\lim }\frac{1}{m}{J}_{33},\hfill \\ a_{23}=a_{32}=0.\hfill \end{array}$$

Proof. 

We can use the asymptotic properties of MLEs to prove it. ☐

Theorem 2.

As $n\to \infty$ and $m\to \infty$ and $\frac{n}{m}\to p$, then

$$\begin{array}{c}\sqrt{n}(\widehat{R}-R)\to N(0,B),\hfill \end{array}$$

where

$$\begin{array}{c}\hfill B=\frac{1}{u_A{({\theta }_1+{\theta }_2)}^4}\left[{\theta }_1^2(a_{11}{a}_{22}-a_{12}{a}_{21})+{\theta }_2^2(a_{11}{a}_{33}-a_{13}{a}_{31})-2{\theta }_1{\theta }_2{a}_{12}{a}_{13}\right],\end{array}$$

(30)

and $u_A=a_{11}{a}_{22}{a}_{33}-a_{13}{a}_{22}{a}_{31}-a_{12}{a}_{21}{a}_{33}$.

Proof. 

By using Theorem 1 and the delta method, we immediately derive the asymptotic distribution of $\widehat{R}$ as follows:

$$\begin{array}{c}\sqrt{n}(\widehat{R}-R)\to N(0,B),\hfill \end{array}$$

where

$$\begin{array}{c}\hfill B=c_A^t{A}^{-1}{c}_A,\end{array}$$

with

$$\begin{array}{c}\hfill c_A=\left(\begin{array}{c}\frac{\partial R}{\partial \alpha }\\ \frac{\partial R}{\partial {\theta }_1}\\ \frac{\partial R}{\partial {\theta }_2}\end{array}\right)=\frac{1}{{({\theta }_1+{\theta }_2)}^2}\left(\begin{array}{c}0\\ -{\theta }_2\\ {\theta }_1\end{array}\right),\end{array}$$

$$\begin{array}{c}\hfill A^{-1}=\frac{1}{u_A}\left(\begin{array}{ccc}{a}_{22}{a}_{33}& -a_{12}{a}_{33}& -a_{22}{a}_{13}\\ -a_{21}{a}_{33}& a_{11}{a}_{33}-a_{13}{a}_{31}& a_{21}{a}_{13}\\ -a_{22}{a}_{31}& a_{12}{a}_{31}& a_{11}{a}_{22}-a_{12}{a}_{21}\end{array}\right),\end{array}$$

and

$$\begin{array}{c}\hfill u_A=a_{11}{a}_{22}{a}_{33}-a_{13}{a}_{22}{a}_{31}-a_{12}{a}_{21}{a}_{33}.\end{array}$$

Therefore,

$$\begin{array}{c}\hfill B=c_A^t{A}^{-1}{c}_A=\frac{1}{u_A{({\theta }_1+{\theta }_2)}^4}\left[(a_{11}{a}_{33}-a_{13}{a}_{31}){\theta }_2^2+(a_{11}{a}_{22}-a_{12}{a}_{21}){\theta }_1^2-2{\theta }_1{\theta }_2{a}_{12}{a}_{13}\right].\end{array}$$

☐

We can derive the $100(1-\gamma )\%$ confidence interval for R using Theorem 2 as

$$\begin{array}{c}\hfill \left(\widehat{R}-z_{1-\frac{\gamma }{2}}\frac{\widehat{B}}{\sqrt{n}},\widehat{R}+z_{1-\frac{\gamma }{2}}\frac{\widehat{B}}{\sqrt{n}}\right),\end{array}$$

where $z_r$ is 100$\gamma$th percentile of N (0, 1). The confidence interval of R can be derived by using the estimate of B in (30). To estimate B, we use the MLEs of $\alpha$, ${\theta }_1$ and ${\theta }_2$ and the following:

$$\begin{array}{c}{\widehat{a}}_{11}=\frac{1}{{\widehat{\alpha }}^2}\left[(1+p)(1+{\gamma }^{″}(2))+p{(\ln {\widehat{\theta }}_1)}^2+{(\ln {\widehat{\theta }}_2)}^2+2{\gamma }^{\prime }(2)(p\ln {\widehat{\theta }}_1+\ln {\widehat{\theta }}_2)\right],\hfill \\ {\widehat{a}}_{12}={\widehat{a}}_{21}=-\frac{\sqrt{p}}{{\widehat{\theta }}_1\widehat{\alpha }}\left[\ln {\widehat{\theta }}_1+{\gamma }^{\prime }(2)\right],\hfill \\ {\widehat{a}}_{13}={\widehat{a}}_{31}=-\frac{1}{{\widehat{\theta }}_2\widehat{\alpha }}\left[\ln {\widehat{\theta }}_2+{\gamma }^{\prime }(2)\right],\hfill \\ {\widehat{a}}_{22}=\frac{1}{{\widehat{\theta }}_1^2},\hfill \\ {\widehat{a}}_{33}=\frac{1}{{\widehat{\theta }}_2^2}.\hfill \end{array}$$

5.2. Bootstrap Confidence Intervals

In this subsection, two confidence intervals based on the non-parametric bootstrap methods are proposed: (i) the percentile bootstrap method (Boot-p) (see [19]) and (ii) the bootstrap-t method (Boot-t) (see [20]).

The algorithms for conducting the confidence intervals of R are presented as follows:

(i)

Boot-p method:

 

Step 1: From the sample $x_1,x_2,...,x_n$ and $y_1,y_2,...,y_m$ compute $\widehat{\alpha }$, $\widehat{{\theta }_1}$ and $\widehat{{\theta }_2}$.

Step 2: Generate a bootstrap sample $x_1^{\ast },...,x_n^{\ast }$ using $\widehat{\alpha }$ and $\widehat{{\theta }_1}$ and generate a bootstrap sample $y_1^{\ast },...,y_m^{\ast }$ using $\widehat{\alpha }$ and $\widehat{{\theta }_2}$. Based on $x_1^{\ast },...,x_n^{\ast }$ and $y_1^{\ast },...,y_m^{\ast }$, compute the estimate of R, say ${\widehat{R}}^{\ast }$.

Step 3: Repeat Step 2 NBOOT times.

Step 4: Let $H_1(x)=P({\widehat{R}}^{\ast }\le x)$ be the cumulative distribution function(cdf) of ${\widehat{R}}^{\ast }$. Define ${\widehat{R}}_{Boot-p}=H_1^{-1}(x)$ for a given x. Thus, the approximate $100(1-z)\%$ C.I. of R is given by:

$$\begin{array}{c}\hfill \left({\widehat{R}}_{Boot-p}(\frac{z}{2}),{\widehat{R}}_{Boot-p}(1-\frac{z}{2})\right).\end{array}$$

Note: In this paper, ${\widehat{R}}^{\ast }$ can be computed using (14) in Step 2.

(ii)

Boot-t method:

 

Step 1: From the sample $x_1,...,x_n$ and $y_1,...,y_m$, compute $\widehat{\alpha }$, $\widehat{{\theta }_1}$ and $\widehat{{\theta }_2}$.

Step 2: Generate a bootstrap sample $x_1^{\ast },...,x_n^{\ast }$ using $\widehat{\alpha }$ and $\widehat{{\theta }_1}$ and generate a bootstrap sample $y_1^{\ast },...,y_m^{\ast }$ using $\widehat{\alpha }$ and $\widehat{{\theta }_2}$. Based on $x_1^{\ast },...,x_n^{\ast }$ and $y_1^{\ast },...,y_m^{\ast }$, compute the estimate of R as ${\widehat{R}}^{\ast }$ and the statistic is

$$\begin{array}{c}\hfill T^{\ast }=\frac{\sqrt{n}({\widehat{R}}^{\ast }-\widehat{R})}{\sqrt{Var(\widehat{R^{\ast }})}}.\end{array}$$

Step 3: Repeat Step 2 NBOOT times.

Step 4: Let $H_2(x)=P(T^{\ast }\le x)$ be the cumulative distribution function(cdf) of $T^{\ast }$. Define ${\widehat{R}}_{Boot-t}=\widehat{R}+n^{-\frac{1}{2}}\sqrt{Var(\widehat{R})}{H}_2^{-1}(x)$ for a given x. Thus, the approximate $100(1-z)\%$ C.I. of R is given by:

$$\begin{array}{c}\hfill \left({\widehat{R}}_{Boot-t}(\frac{z}{2}),{\widehat{R}}_{Boot-t}(1-\frac{z}{2})\right).\end{array}$$

Note: In this paper, $Var({\widehat{R}}^{\ast })$ can be obtained using Theorem 2 in Step 2.

6. Bayesian Inference on $\mathit{R}$

In this section, the Bayes estimate of R can be obtained under assumption that the shape parameter $\alpha$ and scale parameters ${\theta }_1$ and ${\theta }_2$ are random variables. According to the likelihood function in Section 3, we note that $\alpha$ has a positive exponent and ${\theta }_1$ and ${\theta }_2$ have negative exponents, respectively, so we can assume that ${\theta }_1$ and ${\theta }_2$ have independent Inverse Gamma pdf and $\alpha$ follows Gamma distribution. We choose this family such that prior-to-posterior updating yields a posterior that is also in the family:

$$\begin{array}{c}\hfill \pi ({\theta }_1)=\frac{b_1^{a_1}}{\gamma (a_1)}{\theta }_1^{-(1+a_1)}{e}^{-\frac{b_1}{{\theta }_1}},{\theta }_1>0,\end{array}$$

(31)

$$\begin{array}{c}\hfill \pi ({\theta }_2)=\frac{b_2^{a_2}}{\gamma (a_2)}{\theta }_2^{-(1+a_2)}{e}^{-\frac{b_2}{{\theta }_2}},{\theta }_2>0,\end{array}$$

(32)

where all the hyper-parameters $a_i$ and $b_i$$(i=1,2)$ are assumed to be known and non-negative.

The prior density function of $\alpha$ is denoted as $\pi (\alpha )$, and we assume that it has a Gamma (0, 1) distribution.

We have the likelihood function based on the above assumptions as

$$\begin{array}{c}\hfill L(data|\alpha ,{\theta }_1,{\theta }_2)={\alpha }^{m+n}{\theta }_1^{-n}{\theta }_2^{-m}(\prod _{i=1}^n{x}_i^{-\alpha -1})(\prod _{j=1}^m{y}_j^{-\alpha -1})e^{-\frac{1}{{\theta }_1}{\sum }_{i=1}^n{x}_i^{-\alpha }}{e}^{-\frac{1}{{\theta }_2}{\sum }_{j=1}^m{y}_j^{-\alpha }}.\end{array}$$

(33)

The joint density of the data, $\alpha$, ${\theta }_1$ and ${\theta }_2$ becomes

$$\begin{array}{c}\hfill P(data,\alpha ,{\theta }_1,{\theta }_2)=L(data|\alpha ,{\theta }_1,{\theta }_2)\pi (\alpha )\pi ({\theta }_1)\pi ({\theta }_2).\end{array}$$

(34)

Therefore, the joint posterior density of $\alpha$, ${\theta }_1$ and ${\theta }_2$ given the data is

$$\begin{array}{c}\hfill P(\alpha ,{\theta }_1,{\theta }_2|data)=\frac{P(data,\alpha ,{\theta }_1,{\theta }_2)}{{\int }_0^{\infty }{\int }_0^{\infty }{\int }_0^{\infty }L(data,\alpha ,{\theta }_1,{\theta }_2)\mathrm{d}\alpha \mathrm{d}{\theta }_1\mathrm{d}{\theta }_2}.\end{array}$$

(35)

Since the expression (35) cannot be written in a closed form, the Bayes estimate of R and the corresponding credible interval of R are derived adopting the Gibbs sampling technique:

$$\begin{array}{cc}\hfill P(\alpha ,{\theta }_1,{\theta }_2|data)\propto & {\alpha }^{m+n-1}{\theta }_1^{-(n+1+a_1)}{\theta }_2^{-(m+1+a_2)}({\prod }_{i=1}^n{x}_i^{-\alpha -1})({\prod }_{j=1}^m{y}_j^{-\alpha -1})\\ & \exp \left\{-\frac{1}{{\theta }_1}({\sum }_{i=1}^n{x}_i^{-\alpha }+b_1)-\frac{1}{{\theta }_2}({\sum }_{j=1}^m{y}_j^{-\alpha }+b_2)-\alpha \right\}.\end{array}$$

The posterior pdfs of $\alpha$, ${\theta }_1$ and ${\theta }_2$ can be obtained based on the expression $P(\alpha ,{\theta }_1,{\theta }_2|data)$ as the following:

$$\begin{array}{c}{\theta }_1|\alpha ,{\theta }_2,data\sim IG(n+a_1,b_1+\sum _{i=1}^n{x}_i^{-\alpha }),\hfill \\ {\theta }_2|\alpha ,{\theta }_1,data\sim IG(m+a_2,b_2+\sum _{j=1}^m{y}_j^{-\alpha }),\hfill \\ f_{\alpha }(\alpha |{\theta }_1,{\theta }_2,data)\propto {\alpha }^{m+n-1}\prod _{i=1}^n{x}_i^{-\alpha -1}\prod _{j=1}^m{y}_j^{-\alpha -1}\exp \left\{-\frac{1}{{\theta }_1}(\sum _{i=1}^n{x}_i^{-\alpha }+b_1)-\frac{1}{{\theta }_2}(\sum _{j=1}^m{y}_j^{-\alpha }+b_2)-\alpha \right\}.\hfill \end{array}$$

The posterior pdf of $\alpha$ is not known. We use the Metropolis–Hastings method with normal proposal distribution to generate random numbers from this distribution.

The algorithm of Gibbs sampling is described as follows:

Step 1: Start with an initial guess $({\alpha }^{(0)},{\theta }_1^{(0)},{\theta }_2^{(0)})$.

Step 2: Set $t=1$.

Step 3: Generate ${\theta }_1^{(t)}$ form $IG(n+a_1,b_1+{\sum }_{i=1}^n{x}_i^{-{\alpha }^{(t-1)}})$.

Step 4: Generate ${\theta }_2^{(t)}$ form $IG(m+a_2,b_2+{\sum }_{j=1}^m{y}_j^{-{\alpha }^{(t-1)}})$.

Step 5: Using the Metropolis–Hastings method, generate ${\alpha }^{(t)}$ from $f_{\alpha }({\alpha }^{(t-1)}|{\theta }_1^{(t)},{\theta }_2^{(t)},data)$.

• Generate a new candidate value $\delta$ from $N(\ln {\alpha }^{(t-1)},1)$.
• Set ${\alpha }^{\prime }=\exp (\delta )$.
• calculate $p=\min (1,\frac{f_{\alpha }({\alpha }^{\prime }|{\theta }_1^{(t)},{\theta }_2^{(t)},data)}{f_{\alpha }({\alpha }^{(t-1)}|{\theta }_1^{(t)},{\theta }_2^{(t)},data)})$.
• Update ${\alpha }^{(t)}$ with probability p; otherwise set ${\alpha }^{(t)}={\alpha }^{(t-1)}$.
   
  Step 6: Compute $R^{(t)}$ from (14).
  Step 7: Set $t=t+1$.
  Step 8: Repeat step 3–7, M times.

The approximate posterior mean of R is

$$\begin{array}{c}\hfill \widehat{E}(R|data)=\frac{1}{M}\sum _{t=1}^M{R}^{(t)},\end{array}$$

(36)

and the approximate posterior variance of R is

$$\begin{array}{c}\hfill \widehat{Var}(R|data)=\frac{1}{M}\sum _{t=1}^M{(R^{(t)}-\widehat{E}(R|data))}^2.\end{array}$$

(37)

Using the method proposed by [21], we immediately construct the $100(1-\gamma )\%$ highest posterior density (HPD) credible interval as

$$\begin{array}{c}\hfill (R_{\left[\frac{\gamma }{2}M\right]},R_{\left[(1-\frac{\gamma }{2})M\right]}),\end{array}$$

(38)

where $R_{\left[\frac{\gamma }{2}M\right]}$ and $R_{\left[(1-\frac{\gamma }{2})M\right]}$ are the $\left[\frac{\gamma }{2}M\right]$-th smallest integer and the $\left[(1-\frac{\gamma }{2})M\right]$-th smallest integer of $\{R_t,t=1,2,...,M\}$.

7. Numerical Simulations and Data Analysis

In this section, we present a Monte Carlo simulation study and a real data set to illustrate different estimation methods proposed in the preceding sections.

7.1. Numerical Simulations Study

Since we cannot compare the performances of the different methods theoretically, some simulation results to compare the performances of the different methods are presented. We mainly compute the biases and mean square errors (MSEs) of the MLEs, AMLEs and Bayes estimates. The asymptotic C.I. of R and two C.I.s based on the non-parametric bootstrap methods are obtained. We also conduct the Bayes estimates and HPD credible intervals of R. Here, we assume that $a1=a2=b1=b2=0.0001$. We consider sample sizes $(n,m)=(10,10),(20,15),(25,25),(30,40),(50,50).$ For parameter values, we assume that ${\theta }_2=1$, ${\theta }_1=0.5,1,1.5,2,3$ and $\alpha =2$. All the results are based on the 1000 replications. For the bootstrap methods, we compute the confidence intervals based on 300 resampling. The Bayes estimates are based on 1000 sampling, namely, M = 1000. In each case, the nominal level for the C.I.s or the credible intervals is 0.95.

We also obtain the average biases and MSEs of the MLEs, AMLEs and Bayes estimates over 1000 replications in

Table 1. Biases and mean square errors(MSEs) of the maximum likelihood estimate (MLE), approximate maximum likelihood estimate (AMLE) and Bayes estimate of R, when $\alpha =2$, ${\theta }_2=1$ and for different values of ${\theta }_1$.

$(\mathit{n},\mathit{m})$	${\mathit{\theta }}_1=0.5$	${\mathit{\theta }}_1=1$	${\mathit{\theta }}_1=1.5$	${\mathit{\theta }}_1=2$	${\mathit{\theta }}_1=3$
$(10,10)$	0.0919(0.0165)	−0.0094(0.0079)	−0.0082(0.0108)	−0.0223(0.0137)	−0.0098(0.0090)
	0.0881(0.0158)	−0.0088(0.0078)	−0.0207(0.0137)	−0.0145(0.0143)	−0.0054(0.0091)
	0.0721(0.0120)	−0.0084(0.0065)	0.0018(0.0083)	−0.0036(0.0111)	0.0138(0.0078)
$(20,15)$	0.0281(0.0078)	0.0095(0.0074)	−0.0098(0.0079)	−0.0188(0.0091)	−0.0127(0.0065)
	0.0317(0.0078)	0.0061(0.0075)	−0.0068(0.0078)	−0.0161(0.0095)	−0.0079(0.0065)
	0.0207(0.0086)	0.0090(0.0064)	0.0034(0.0079)	−0.0056(0.0078)	0.0043(0.0060)
$(25,25)$	0.0132(0.0037)	0.0060(0.0025)	0.0025(0.0067)	0.0057(0.0045)	−0.0270(0.0049)
	0.0132(0.0038)	0.0062(0.0024)	0.0040(0.0068)	0.0047(0.0045)	−0.0242(0.0049)
	0.0055(0.0036)	0.0061(0.0022)	0.0076(0.0063)	0.0128(0.0041)	−0.0165(0.0045)
$(30,40)$	−0.0094(0.0014)	−0.0039(0.0010)	−0.0714(0.0051)	−0.0048(0.0024)	−0.0510(0.0038)
	−0.0092(0.0014)	−0.0064(0.0011)	−0.0760(0.0051)	−0.0060(0.0023)	−0.0504(0.0038)
	−0.0135(0.0015)	−0.0040(0.0009)	−0.065(0.0042)	0.0058(0.0023)	−0.0302(0.0035)
$(50,50)$	−0.0182(0.0005)	−0.0087(0.0007)	−0.0033(0.0038)	0.0007(0.0017)	−0.0035(0.0027)
	−0.0192(0.0005)	−0.0080(0.0006)	−0.0033(0.0039)	0.0012(0.0017)	−0.0022(0.0027)
	−0.0217(0.0006)	−0.0079(0.0006)	−0.0008(0.0036)	0.0051(0.0015)	0.0023(0.0025)

In each cell, the average biases are provided and corresponding MSEs are presented within brackets. The first to the third row corresponds the results for MLEs, AMLEs and Bayes estimates respectively.

. From Table 1, we can find that the Bayes estimates are almost as efficient as the MLEs and the AMLEs for all sample sizes. Interestingly, in most of the cases, the MSEs of the Bayes estimate are smaller than the MSEs of the MLEs or AMLEs. We can find that the biases and MSEs of the MLEs and AMLEs are very close. As the sample size $(n,m)$ increases, the MSEs of the estimates decrease as expected.

Table 2. Average confidence lengths and coverage probabilities.

$(\mathit{n},\mathit{m})$	${\mathit{\theta }}_1=0.5$	${\mathit{\theta }}_1=1$	${\mathit{\theta }}_1=1.5$	${\mathit{\theta }}_1=2$	${\mathit{\theta }}_1=3$
$(10,10)$	0.3879(0.90)	0.4205(0.91)	0.4074(0.90)	0.3880(0.93)	0.3377(0.89)
	0.3894(0.90)	0.4206(0.92)	0.4079(0.90)	0.3907(0.93)	0.3403(0.90)
	0.4044(0.89)	0.4526(0.92)	0.4303(0.91)	0.4042(0.92)	0.3469(0.91)
	0.5179(0.90)	0.4526(0.90)	0.5425(0.90)	0.5240(0.90)	0.3939(0.90)
	0.3997(0.95)	0.4113(0.95)	0.4028(0.93)	0.3931(0.95)	0.3546(0.94)
$(20,15)$	0.3201(0.92)	0.3480(0.92)	0.3427(0.92)	0.3176(0.92)	0.2862(0.91)
	0.3205(0.93)	0.3480(0.93)	0.3430(0.94)	0.3186(0.93)	0.2873(0.91)
	0.3305(0.92)	0.3657(0.93)	0.3577(0.94)	0.3275(0.93)	0.2915(0.92)
	0.3810(0.90)	0.4089(0.91)	0.4020(0.94)	0.3814(0.91)	0.3531(0.90)
	0.3281(0.95)	0.3427(0.93)	0.3390(0.95)	0.3207(0.96)	0.2953(0.96)
$(25,25)$	0.2501(0.93)	0.2729(0.94)	0.2634(0.93)	0.2509(0.96)	0.2244(0.95)
	0.2502(0.93)	0.2729(0.95)	0.2637(0.94)	0.2512(0.95)	0.2245(0.95)
	0.2510(0.93)	0.2774(0.94)	0.2696(0.92)	0.2555(0.96)	0.2250(0.96)
	0.2756(0.93)	0.2961(0.93)	0.2765(0.92)	0.2755(0.93)	0.2354(0.95)
	0.2537(0.96)	0.2700(0.96)	0.1884(0.95)	0.2512(0.95)	0.2286(0.96)
$(30,40)$	0.2278(0.96)	0.2496(0.94)	0.2437(0.94)	0.2261(0.96)	0.2059(0.95)
	0.2279(0.96)	0.2496(0.95)	0.2436(0.95)	0.2263(0.96)	0.2062(0.94)
	0.2296(0.95)	0.2522(0.95)	0.2482(0.95)	0.2275(0.97)	0.2050(0.94)
	0.2482(0.94)	0.2682(0.94)	0.2588(0.95)	0.2458(0.92)	0.2267(0.93)
	0.2284(0.97)	0.2470(0.95)	0.2423(0.97)	0?2272(0.97)	0.2090(0.95)
$(50,50)$	0.1788(0.96)	0.1948(0.96)	0.1887(0.95)	0.1780(0.97)	0.1553(0.96)
	0.1788(0.96)	0.1948(0.96)	0.1887(0.96)	0.1781(0.96)	0.1554(0.95)
	0.1789(0.96)	0.1952(0.95)	0.1890(0.94)	0.0812(0.96)	0.1598(0.96)
	0.1848(0.94)	0.2013(0.94)	0.1942(0.92)	0.1845(0.95)	0.1588(0.95)
	0.1784(0.97)	0.1935(0.96)	0.1880(0.96)	0.1778(0.96)	0.1620(0.97)

In each cell, the average confidence lengths are provided and the corresponding coverage probabilities are given within brackets. The first to the fifth row corresponds the results for MLEs, AMLEs, Boot-p method, Boot-t method and Bayes estimates respectively.

reports the results of 95% asymptotic C.I. of R, we also obtain C.I.s based on the bootstrap methods and the HPD credible interval. We represent the results of the average confidence/credible lengths and the coverage probabilities. From Table 2, the coverage probabilities reach the nominal level 95% with the increase of the sample sizes. We observe that the MLE method is the most valid procedure to obtain the confidence intervals. The AMLEs and the Bayes estimates are the second best confidence intervals. Interestingly, we find that the HPD credible intervals provide the most highest coverage probabilities. The Boot-p confidence intervals perform better than Boot-t confidence intervals, in terms of coverage probabilities. One point we should know is that the bootstrap method depends on the number of resampling. For small sample sizes $(n,m)$, the coverage probabilities for the MLEs and AMLEs are less than the nominal levels, with the increase of sample sizes $(n,m)$, they perform well.

7.2. Data Analysis

We consider a real data set to illustrate the methods of inference discussed in this article. These strength data sets in

Table 3. Strength measured in GPA for carbon fibers tested under tension at gauge lengths of 20 mm.

1.312	1.314	1.479	1.552	1.700	1.803	1.861	1.865	1.944	1.958
1.966	1.997	2.006	2.021	2.027	2.055	2.063	2.098	2.140	2.179
2.224	2.240	2.253	2.270	2.272	2.274	2.301	2.301	2.359	2.382
2.382	2.426	2.434	2.435	2.478	2.490	2.511	2.514	2.535	2.554
2.566	2.570	2.586	2.629	2.633	2.642	2.648	2.684	2.697	2.726
2.770	2.773	2.800	2.809	2.818	2.821	2.848	2.880	2.954	3.012
3.067	3.084	3.090	3.096	3.128	3.233	3.433	3.585	3.585

and

Table 4. Strength measured in GPA for carbon fibers tested under tension at gauge lengths of 10 mm.

1.901	2.132	2.203	2.228	2.257	2.350	2.361	2.396	2.397	2.445
2.454	2.474	2.518	2.522	2.525	2.532	2.575	2.614	2.616	2.618
2.624	2.659	2.675	2.738	2.740	2.856	2.917	2.928	2.937	2.937
2.977	2.996	3.030	3.125	3.139	3.145	3.220	3.223	3.235	3.243
3.264	3.272	3.294	3.332	3.346	3.377	3.408	3.435	3.493	3.501
3.537	3.554	3.562	3.628	3.852	3.871	3.886	3.971	4.024	4.027
4.225	4.395	5.020

were analyzed previously by [3,22]. We know that if the random variable X follows $W(\alpha ,\theta )$, the random variable $T=\frac{1}{X}$ has the $IW(\alpha ,\theta )$. Hence, we have the following data sets from the inverse Weibull distribution. These data are presented in

Table 5. Transformed Data Set 1.

0.762	0.761	0.676	0.644	0.588	0.555	0.537	0.536	0.514	0.511
0.509	0.501	0.499	0.495	0.493	0.487	0.485	0.477	0.467	0.459
0.450	0.446	0.444	0.441	0.440	0.440	0.435	0.435	0.424	0.420
0.420	0.412	0.411	0.411	0.404	0.402	0.398	0.398	0.394	0.392
0.390	0.389	0.387	0.380	0.380	0.379	0.378	0.373	0.371	0.367
0.361	0.361	0.357	0.356	0.355	0.354	0.351	0.347	0.339	0.332
0.326	0.324	0.324	0.323	0.320	0.309	0.291	0.279	0.279

and

Table 6. Transformed Data Set 2.

0.526	0.469	0.454	0.449	0.443	0.426	0.424	0.417	0.417	0.409
0.407	0.404	0.397	0.397	0.396	0.395	0.388	0.383	0.382	0.382
0.381	0.376	0.374	0.365	0.365	0.350	0.343	0.342	0.340	0.340
0.336	0.334	0.330	0.320	0.319	0.318	0.311	0.310	0.309	0.308
0.306	0.306	0.304	0.300	0.299	0.296	0.293	0.291	0.286	0.286
0.283	0.281	0.281	0.276	0.260	0.258	0.257	0.252	0.249	0.248
0.237	0.228	0.199

. We analyze the data by adding 0.5 from both data sets. We fit the inverse Weibull models to the two data sets separately. The estimated shape and scale parameters, log-likelihood values, Kolmogorov-Smirnov (K-S) distances and corresponding p-values are presented in

Table 7. Scale parameter, shape parameter, log-likelihood, K-S distances and p-values of the fitted inverse Weibull models to Data Sets 1 and 2.

Data Set	Scale Parameter	Shape Parameter	Log-Likelihood	K-S	p-Value
1	4.9497	12.6152	71.8967	0.0417	0.9997
2	19.0814	13.6228	79.4095	0.0846	0.7572

. The expected frequencies and the observed frequencies based on the fitted models are also presented in the

Table 8. Expected frequencies and observed frequencies for modified Data Set 1 when fitting the inverse Weibull model.

Intervals	Expected Frequencies	Observed Frequencies
0.00–0.25	0.03	0
0.25–0.40	32.15	33
0.40–0.50	24.21	24
0.50–0.70	11.23	10
0.70–∞	1.38	2

and

Table 9. Expected frequencies and observed frequencies for modified Data Set 2 when fitting the inverse Weibull model.

Intervals	Expected Frequencies	Observed Frequencies
0.00–0.19	0.01	0
0.19–0.30	21.06	19
0.30–0.40	29.48	32
0.40–0.50	9.23	11
0.50–∞	3.22	1;

. We also obtain the chi-square values of 5.9914 and 5.9915. Obviously, the inverse Weibull model fits very well to Data Set 1 and Data Set 2.

The K-S values and the corresponding p-values in

Table 10. Scale parameter, shape parameter, log-likelihood, K-S distances and p-values of the fitted inverse Weibull models to Data Sets 1 and 2. Here, we assume that the two shape parameters are identical.

Data Set	Scale Parameter	Shape Parameter	Log-Likelihood	K-S	p-Value
1	5.3471	13.0933	71.8159	0.0424	0.9996
2	16.7168	13.0933	79.3215	0.0732	0.8878

show that the inverse Weibull models with equal shape parameters fit reasonably well to the modified data sets. It is clear that we cannot reject the null hypothesis that the two shape parameters are equal.

Based on Equations (14) and (29), we provide that the MLE and AMLE of R are $0.7576$ and 0.7571. The 95% confidence intervals of MLE, AMLE, Boot-p method and Boot-t method are obtained as (0.6917, 0.8235), (0.6911, 0.8231), (0.6993, 0.8197), (0.7015, 0.8421), respectively.

The Bayesian estimate of R is also presented based on Equation (36). In the previous sections, we assume that ${\theta }_1$ and ${\theta }_2$ have independent IG priors, $\alpha$ has a Gamma (0, 1) prior and $a1=a2=b1=b2=0.0001$. Under certain assumptions, we can conduct the Bayesian estimate of R as 0.7437 and the 95% HPD credible interval of R can be obtained as (0.6690, 0.8102).

8. Inference on R with All Different Parameters

In the sections above, we assume that the shape parameters are taken to be equal. In order to expand the paper, in this section, we study the inference of R with all different parameters. We consider the problem of estimating $R=P(Y<X)$ under the assumption that $X\sim IW({\alpha }_1,{\theta }_1)$ and $Y\sim IW({\alpha }_2,{\theta }_2)$. Then, it can be easily calculated that

$$\begin{array}{c}\hfill R=P(Y<X)=1-{\int }_0^{\infty }\frac{{\alpha }_2}{{\theta }_2}{y}^{-{\alpha }_2-1}{e}^{-\frac{y^{-{\alpha }_2}}{{\theta }_2}}{e}^{-\frac{y^{-{\alpha }_1}}{{\theta }_1}}dy.\end{array}$$

(39)

To computer the MLE of R, we first obtain the MLEs of ${\theta }_1$ and ${\theta }_2$. Suppose $X_1,X_2,\cdots ,X_n$ is a random sample from $IW({\alpha }_1,{\theta }_1)$ and $Y_1,Y_2,...,Y_m$ is a random sample from $IW({\alpha }_2,{\theta }_2)$. The joint likelihood function is:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill l({\alpha }_1,{\theta }_1,{\alpha }_2,{\theta }_2)=\frac{{\alpha }_1^n}{{\theta }_1^n}\frac{{\alpha }_2^m}{{\theta }_2^m}(\prod _{i=1}^n{x}_i^{-{\alpha }_1-1})(\prod _{j=1}^m{y}_j^{-{\alpha }_2-1})e^{-{\sum }_{i=1}^n\frac{x_i^{-{\alpha }_1}}{{\theta }_1}}{e}^{-{\sum }_{j=1}^m\frac{y_j^{-{\alpha }_2}}{{\theta }_2}}.\end{array}\end{array}$$

Then, the log-likelihood function is

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill L({\alpha }_1,{\theta }_1,{\alpha }_2,{\theta }_2)=& n\ln {\alpha }_1+m\ln {\alpha }_2-n\ln {\theta }_1-m\ln {\theta }_2-({\alpha }_1+1)\sum _{i=1}^n\ln x_i-({\alpha }_2+1)\sum _{j=1}^m\ln y_j\hfill \\ & -\frac{1}{{\theta }_1}\sum _{i=1}^n{x}_i^{-{\alpha }_1}-\frac{1}{{\theta }_2}\sum _{j=1}^m{y}_j^{-{\alpha }_2}.\hfill \end{array}\end{array}$$

(40)

Then, similar to the previous approaches, we can obtain the point estimates and interval estimates of R by using the MLE, AMLE and Bayesian method. Furthermore, we also can get bootstrap confidence intervals of R.

9. Conclusions

In this paper, we have addressed the estimator of $R=P(Y<X)$ for the inverse Weibull distribution. We assume independent inverse Weibull random variables with equal shape parameters but different scale parameters.

We obtain the maximum likelihood estimator of R and its asymptotic distribution. Note that MLEs do not have explicit forms, and we propose the approximate maximum likelihood estimator of R. The confidence interval of R is obtained using the asymptotic distribution. Two bootstrap confidence intervals are also obtained. By using the Gibbs sampling technique, we present the Bayesian estimator of R and the corresponding credible interval. The Metropolis–Hastings algorithm with the normal proposal distribution is used to generate random numbers from a given density function. Monte Carlo simulations are conducted to compare the proposed methods. Analysis of a real dataset is performed. In the future, we will consider the MLE, AMLE, asymptotic C.I., bootstrap C.I.s and Bayesian inference of $R=P(Y<X)$ for inverse Weibull distribution based on incomplete data such as progressively type-II censored samples.