Inference on Pr(𝑋<𝑌) in the Two-Parameter Weibull Model Based on Records

Abstract

We consider the stress-strength reliability based on record values from the Weibull distribution. The Bayes estimator based on squared error loss and the maximum likelihood estimator are derived and their bias and mean squared error performance are studied. Likelihood-based confidence intervals as well as some bootstrap intervals are developed. We derived also the highest posterior density interval. Simulation studies are conducted to investigate and compare the performance of the estimators and intervals.

1. Introduction

Chandler [1] introduced and studied some properties of record values. Since then a considerable amount of the literature is devoted to the study of records. Ahsanullah [2] and Arnold et al. [3] provided a detailed account of theory of records and the inference issues associated with records. In this paper, we will consider confidence interval estimation of the stress-strength reliability based on record data when the underlying distribution is Weibull. The Weibull model is quite popular in analyzing reliability and life-testing data. Its flexibility and capability of modeling various forms of failure mechanisms and hazard function types gave him an important role in reliability literature; see Johnson et al. [4] and Murthy et al. [5]. The probability density function (pdf) and the cumulative distribution function (cdf) of the Weibull distribution 𝑊(𝜃,𝛽) are given by 𝑓(𝑥)=𝜃𝛽𝑥𝐹^𝛽−1𝑒(𝑥)=^{−𝜃𝑥𝛽},1−𝑥>𝑒0,_{−𝜃𝑥𝛽}𝜃>0,𝛽>0.(1.1)

Let 𝑋₁,𝑋₂,… be an infinite sequence of 𝑖𝑖𝑑 random variables. An observation 𝑋_𝑗 is called a record if its value is greater than all previous observations, that is 𝑋_𝑗>𝑋_𝑖 for every 𝑖<𝑗. We want to estimate the stress-strength reliability Pr(𝑋<𝑌) using the data on records. The stress-strength reliability Pr(𝑋<𝑌) arises in life-testing experiments when 𝑋 and 𝑌 represent the lifetimes of two devices, and it gives the probability that the device with life time 𝑋 fails before the other. As an example, Hall [6] studied a situation where the breakdown voltage 𝑌 of a capacitor must exceed the voltage output 𝑋 of a power supply in order for the component to work properly. Weerahandi and Johnson [7] considered another example on rocket motors. This probability has other interpretations in other disciplines, for example, in medical sciences; it is used as a measure of treatment effectiveness in studies involving comparison between control and treatment groups. Various other examples may be found in Kotz et al. [8].

Estimation of the stress-strength reliability based on record data was considered by Baklizi [9] for the exponential distribution with record values and by Baklizi [10] for the generalized exponential distribution. Kundu and Gupta [11] investigated this problem for simple random samples from the Weibull distribution. In Section 2 we derive the maximum likelihood estimator and the associated large sample intervals. Bayesian procedures based on records are derived in Section 3. Simulations to investigate the performance of the asymptotic inference procedures and to compare them with the bootstrap intervals are described in Section 4. The results and conclusions are given in Section 5.

2. Likelihood Inference

We shall consider the case when the shape parameters are equal. Let 𝑋∼𝑊(𝜃₁,𝛽) and 𝑌∼𝑊(𝜃₂,𝛽) be independent random variables. Let 𝑅=Pr(𝑋<𝑌) be the stress strength reliability. Let 𝑟=𝑟₀,…,𝑟_𝑛 be a set of records from 𝑊(𝜃₁,𝛽), and let 𝑠=𝑠₀,…,𝑠_𝑚 be an independent set of records from 𝑊(𝜃₂,𝛽). The likelihood functions are given by [3] as follows: 𝐿₁𝛽,𝜃₁∣𝑟=𝑓𝑟_𝑛^𝑛−1_𝑖=0𝑓𝑟_𝑖𝐿1−𝐹𝑟_𝑖₂,𝛽,𝜃₂∣𝑠=𝑔𝑠_𝑚^𝑚−1_𝑖=0𝑔𝑠_𝑖1−𝐺𝑠_𝑖,(2.1) where 𝑓 and 𝐹 are the pdf and cdf of 𝑋∼𝑊(𝜃₁,𝛽) and 𝑔 and 𝐺 are the pdf and cdf of 𝑌∼𝑊(𝜃₂,𝛽). The likelihood function of (𝛽,𝜃₁,𝜃₂) based on 𝑟,𝑠 is given by 𝐿𝛽,𝜃₁,𝜃₂∣𝑟,𝑠=𝜃₁𝛽^𝑛+1𝑒^{−𝜃1𝑟𝛽𝑛𝑛}_𝑖=0𝑟_𝑖𝛽−1𝜃₂𝛽^𝑚+1𝑒^{−𝜃2𝑠𝛽𝑚𝑚}_𝑖=0𝑠_𝑖𝛽−1.(2.2)

Taking the natural logarithm we get the log-likelihood function 𝑙𝛽,𝜃₁,𝜃₂=(𝑛+𝑚+2)ln𝛽+(𝑛+1)ln𝜃₁+(𝑚+1)ln𝜃₂−𝜃₁𝑟^𝛽𝑛−𝜃₂𝑠^𝛽𝑚+(𝛽−1)^𝑛_𝑖=1ln𝑟_𝑖+^𝑚_𝑖=1ln𝑠_𝑖.(2.3)

The first partial derivatives are given by 𝜕𝑙𝛽,𝜃₁,𝜃₂𝜕𝛽=(𝑛+𝑚𝛽+2)−𝜃¹𝑟^𝛽𝑛ln𝑟_𝑛−𝜃₂𝑠^𝛽𝑚ln𝑠_𝑚+^𝑛_𝑖=1ln𝑟_𝑖+^𝑚_𝑖=1ln𝑠_𝑖,(2.4)𝜕𝑙𝛽,𝜃₁,𝜃₂𝜕𝜃₁=(𝑛+1)𝜃₁𝜕𝑙𝛽,𝜃₁−𝑟^𝛽𝑛,,𝜃₂𝜕𝜃₂=(𝑚+1)𝜃₂−𝑠^𝛽𝑚.(2.5)

Equating these partial derivatives to zero and solving simultaneously we obtain 𝛽̂=(𝑛+𝑚+2)(𝑛+1)ln𝑟_𝑛+(𝑚+𝜃1̂₁)ln𝑠_𝑚−∑^𝑛𝑖=0ln𝑟_𝑖+∑^𝑚𝑖=0ln𝑠_𝑖,=𝑛+1𝑟𝛽_𝑛̂𝜃̂₂,=𝑚+1𝑠𝛽_𝑚̂.(2.6)

Hence the MLE of 𝑅 is given by  𝑅=𝜃̂₂/𝜃̂₁+𝜃̂₂. The study of the exact distribution of 𝑅 is apparently rather complicated, so we will consider the asymptotic distribution. We need the asymptotic joint distribution of 𝛽̂, 𝜃̂₁ and 𝜃̂₂. The second partial derivatives of the log-likelihood function are given by 𝜕²𝑙𝛽,𝜃₁,𝜃₂𝜕𝛽²=−(𝑛+𝑚+2)𝛽²𝜕₂−𝜃₁𝑟^𝛽𝑛ln𝑟_𝑛²−𝜃₂𝑠^𝛽𝑚ln𝑠_𝑚²,𝑙𝛽,𝜃₁,𝜃₂𝜕𝜃²¹=−(𝑛+1)𝜃²¹𝜕₂,𝑙𝛽,𝜃₁,𝜃₂𝜕𝜃²²=−(𝑚+1)𝜃²²𝜕₂,𝑙𝛽,𝜃₁,𝜃₂𝜕𝛽𝜕𝜃₁𝜕₂=−𝑟^𝛽𝑛ln𝑟_𝑛,𝑙𝛽,𝜃₁,𝜃₂𝜕𝛽𝜕𝜃₂𝜕₂=−𝑠^𝛽𝑚ln𝑠_𝑚,𝑙𝛽,𝜃₁,𝜃₂𝜕𝜃₁𝜕𝜃₂=0.(2.7)

The Fisher information matrix is given by 𝐈=⎛⎜⎜⎜⎜⎝𝐼¹¹𝐼₁₂𝐼₁₃𝐼₂₁𝐼₂₂𝐼₂₃𝐼₃₁𝐼₃₂𝐼₃₃⎞⎟⎟⎟⎟⎠=−⎛⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝𝐸𝜕²𝑙𝜕𝛽²𝐸𝜕²𝑙𝜕𝛽𝜕𝜃₁𝐸𝜕²𝑙𝐸𝜕₂𝜕𝛽𝜕𝜃₂𝑙𝜕𝛽𝜕𝜃₁𝐸𝜕²𝑙𝜕𝜃²¹𝐸𝜕²𝑙𝐸𝜕₂𝜕𝜃₁𝜕𝜃₂𝑙𝜕𝛽𝜕𝜃₂𝐸𝜕²𝑙𝜕𝜃₁𝜕𝜃₂𝐸𝜕²𝑙𝜕𝜃²²⎞⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠.(2.8)

Now we will find the entries of the information matrix. We need to find 𝐸(𝑟^𝛽𝑛ln𝑟_𝑛), 𝐸(𝑟^𝛽𝑛(ln𝑟_𝑛)²), 𝐸(𝑠^𝛽𝑚ln𝑠_𝑚) and 𝐸(𝑠^𝛽𝑚(ln𝑠_𝑚)²). Using result (2.4.3) of Arnold et al. [3] we have 𝑈=𝜃₁𝑟^𝛽𝑛∼Γ(𝑛+1,1), and it follows that𝐸𝑟^𝛽𝑛ln𝑟_𝑛=1𝛽𝜃₁𝐸𝜃₁𝑟^𝛽𝑛ln𝜃¹𝑟^𝛽𝑛−ln𝜃₁=1𝛽𝜃₁𝐸(𝑈ln𝑈)−ln𝜃₁=(𝑛+1)𝐸(𝑈)𝛽𝜃₁𝜓(𝑛+2)−ln𝜃₁,(2.9)𝐸𝑟^𝛽𝑛ln𝑟_𝑛²=1𝜃₁𝛽²𝐸𝜃₁𝑟^𝛽𝑛ln𝜃¹𝑟^𝛽𝑛−ln𝜃₁²=1𝜃₁𝛽²𝐸𝜃₁𝑟^𝛽𝑛ln𝜃¹𝑟^𝛽𝑛²−2ln𝜃₁𝑙𝑛𝜃₁𝑟^𝛽𝑛+ln𝜃₁²=1𝜃₁𝛽²𝐸𝑈(ln𝑈)²=(𝑛+1)−2ln𝜃₁𝐸(𝑈ln𝑈)+ln𝜃₁²𝐸(𝑈)𝜃₁𝛽²𝜓(𝑛+2)²+𝜓^(𝑛+2)−2ln𝜃₁𝜓(𝑛+2)+ln𝜃₁².(2.10)

Similarly, since 𝑉=𝜃₂𝑠^𝛽𝑚∼Γ(𝑚+1,1), we have 𝐸𝑠^𝛽𝑚ln𝑠_𝑚=(𝑚+1)𝐸𝑠_𝛽𝑚𝛽𝜃₂𝜓(𝑚+2)−ln𝜃₂,ln𝑠_𝑚²=(𝑚+1)𝜃₂𝛽²𝜓(𝑚+2)²+𝜓′(𝑚+2)−2ln𝜃₂𝜓(𝑚+2)+ln𝜃₂².(2.11)

Let 𝑅=𝜃₁/(𝜃₁+𝜃₂)=ℎ(𝜃₁,𝜃₂), the maximum likelihood estimator of 𝑅 is given by 𝑅=𝜃̂₁/(𝜃̂₁+𝜃̂₂)=ℎ(𝜃̂₁,𝜃̂₂), and it follows that √𝑛𝑅−𝑅=√𝑛ℎ𝜃̂₁,𝜃̂₂−ℎ𝜃₁,𝜃₂−→^𝐷𝑁0,𝜂²,(2.12) as 𝑛→∞, 𝑚→∞ such that 𝑚/𝑛→𝑝, 0<𝑝<1. The asymptotic variance is given by 𝜂²𝜕ℎ𝜃¹,𝜃₂𝜕𝜃₁²var𝜃̂₁+𝜕ℎ𝜃¹,𝜃₂𝜕𝜃₂²var𝜃̂₂+2cov𝜃̂₁𝑝,𝜃̂₂𝜕ℎ𝜃¹,𝜃₂𝜕𝜃₁𝜕ℎ𝜃¹,𝜃₂𝜕𝜃₂.(2.13)

A (1−𝛼)% confidence interval for 𝑅 based on this asymptotic result is given by 𝑅−𝑧^1−𝛼/2𝜂̂,𝑅+𝑧_1−𝛼/2𝜂̂,(2.14) where 𝜂̂ is obtained by substituting 𝑚/𝑛 for 𝑝 and the MLEs of 𝜃₁ and 𝜃₂ in the asymptotic standard deviation 𝜂. Another interval can be obtained by using the matrix of minus the second partial derivatives of the log-likelihood function. This matrix can be used in place of the Fisher information matrix to obtain 𝜂̂̂ as an estimator of the variance of 𝑅. The new interval is given by 𝑅−𝑧^1−𝛼/2𝜂̂̂,𝑅+𝑧_1−𝛼/2𝜂̂̂.(2.15) Many authors suggested that parameter transformation may improve the performance of intervals based on the asymptotic normality of the maximum likelihood estimator. For parameters representing probabilities, the logit transformation seems appropriate. Let 𝜆=ln(𝑅/1−𝑅), and the maximum likelihood of 𝜆 is given by 𝜆̂=ln(𝑅/1−𝑅). The asymptotic variance of √𝑛(𝜆̂−𝜆) is given by 𝜂²(𝑑𝜆/𝑑𝑅)². This variance can be estimated by substituting the maximum likelihood estimator instead of the parameter. A (1−𝛼)% confidence interval for 𝜆 is given by (𝜆̂_𝐿,𝜆̂_𝑈), where 𝜆̂_𝐿=𝜆̂−𝑧_1−𝛼/2𝜂̂𝑑𝜆/𝑑𝑅∣_𝑅=𝑅 and 𝜆̂_𝑈=𝜆̂+𝑧_1−𝛼/2𝜂̂𝑑𝜆/𝑑𝑅∣_𝑅=𝑅, 𝑧_1−𝛼/2 is the (1−𝛼/2) quantile of the standard normal distribution. Therefore A (1−𝛼)% confidence interval for 𝑅 is given by 1+𝑒𝜆𝑒̂^𝐿𝜆̂^𝐿,𝑒𝜆̂^𝑈1+𝑒𝜆̂^𝑈.(2.16)

3. Bayesian Inference

As in Kundu and Gupta [11], we will use the conjugate gamma priors for the scale parameters 𝜃₁ and 𝜃₂ and a squared error loss function. The prior density of 𝜃_𝑗 therefore is given by 𝜋_𝑗𝜃_𝑗∣𝑎_𝑗,𝑏_𝑗=𝑏^𝑎𝑗𝑗Γ𝑎_𝑗𝜃^{𝑎𝑗𝑗−1}𝑒^{−𝑏𝑗𝜃𝑗},𝑗=1,2,(3.1) where 𝑎_𝑗>0, 𝑏_𝑗>0 are the parameters of the prior distributions. Assuming that 𝜃₁ and 𝜃₂ are independent, the joint prior distribution of 𝜃₁ and 𝜃₂ is given by 𝜋(𝜃₁,𝜃₂)=∏^2𝑗=1𝜋_𝑗(𝜃_𝑗). Assume that the prior distribution of 𝛽, denoted by 𝜁(𝛽), has a support on (0,∞). Recall that the likelihood function of (𝛽,𝜃₁,𝜃₂) based on 𝑟=𝑟₀,…,𝑟_𝑛 and 𝑠=𝑠₀,…,𝑠_𝑚 is given by 𝐿𝛽,𝜃₁,𝜃₂∣𝑟,𝑠=𝜃₁𝛽^𝑛+1𝑒^{−𝜃1𝑟𝛽𝑛𝑛}_𝑖=0𝑟_𝑖𝛽−1𝜃₂𝛽^𝑚+1𝑒^{−𝜃2𝑠𝛽𝑚𝑚}_𝑖=0𝑠_𝑖𝛽−1.(3.2)

The joint posterior density function of (𝛽,𝜃₁,𝜃₂) therefore is given by; 𝜋^∗𝛽,𝜃₁,𝜃₂=𝜋𝜃¹,𝜃₂𝜁(𝛽)𝐿𝛽,𝜃₁,𝜃₂∣𝑟,𝑠∭^∞0𝜋𝜃₁,𝜃₂𝜁(𝛽)𝐿𝛽,𝜃₁,𝜃₂∣𝑟,𝑠𝑑𝛽𝑑𝜃₁𝑑𝜃₂,(3.3) where the numerator is given by 𝛽^𝑛+𝑚+2𝜁(𝛽)𝑏^𝑎11Γ𝑎₁𝜃^1𝑛+𝑎1𝑒^{−𝜃1(𝑟𝛽𝑛+𝑏1)𝑛}_𝑖=0𝑟_𝑖𝛽−1𝑏^𝑎2₂Γ𝑎₂𝜃^2𝑚+𝑎2𝑒^{−𝜃2(𝑠𝛽𝑚+𝑏2)𝑚}_𝑖=0𝑠_𝑖𝛽−1.(3.4)

This expression for 𝜋^∗(𝛽,𝜃₁,𝜃₂) is difficult or even impossible to find in closed form. A simulation technique is needed. It is clear that the conditional posterior distributions of 𝜃₁ and 𝜃₂ given 𝛽 are given by 𝜋^∗1𝜃₁∣𝛽=𝑟^𝛽𝑛+𝑏₁^𝑛+𝑎1+1Γ𝑛+𝑎₁+1𝜃^1𝑛+𝑎1𝑒^{−𝜃1(𝑟𝛽𝑛+𝑏1)},(3.5)𝜋^∗2𝜃₂∣𝛽=𝑠^𝛽𝑚+𝑏₂^𝑚+𝑎2+1Γ𝑚+𝑎₂+1𝜃^2𝑚+𝑎2𝑒^{−𝜃2(𝑠𝛽𝑚+𝑏2)}.(3.6)

Note that given 𝛽, 𝜃₁, and 𝜃₂ are independent. From the joint posterior distribution of 𝜃₁, 𝜃₂, and 𝛽 we obtain 𝜁^∗(𝛽)∝𝜁(𝛽)𝛽^𝑛+𝑚+2^𝑛_𝑖=0𝑟_𝑖𝑚_𝑖=0𝑠_𝑖^𝛽𝑟^𝛽𝑛+𝑏₁^{−(𝑛+𝑎1+1)}𝑠^𝛽𝑚+𝑏₂^{−(𝑚+𝑎2+1)}.(3.7)

Taking 𝑎₁=𝑎₂=0 and 𝑏₁=𝑏₂=0 in the prior distributions we obtain 𝜁^∗(𝛽)∝𝜁(𝛽)𝛽^𝑛+𝑚+2𝑟^{𝑛−(𝑛+1)}𝑠^{𝑚𝑛−(𝑚+1)}_𝑖=0𝑟_𝑖𝑚_𝑖=0𝑠_𝑖^𝛽.(3.8)

If we can take 𝜁(𝛽)=𝑒^−𝛽, we obtain 𝜁^∗(𝛽)∝𝛽^𝑛+𝑚+2𝑒⁻¹𝑟^{𝑛−(𝑛+1)}𝑠^{𝑚𝑛−(𝑚+1)}_𝑖=0𝑟_𝑖𝑚_𝑖=0𝑠_𝑖^𝛽.(3.9)

It follows that the posterior distribution of 𝛽 is Γ⎛⎜⎜⎝𝑛+𝑚+3,ln𝑒⁻¹𝑟^{𝑛−(𝑛+1)}𝑠^{𝑚𝑛−(𝑚+1)}_𝑖=0𝑟_𝑖𝑚_𝑖=0𝑠_𝑖⁻¹⎞⎟⎟⎠.(3.10)

We can use the following Monte Carlo method to find approximate point estimates. (1)Generate 𝛽₁ from 𝜁^∗(𝛽) given by (3.10).(2)Generate 𝜃_1,1 and 𝜃_2,1 from 𝜋^∗1(𝜃₁∣𝛽₁) and 𝜋^∗2(𝜃₂∣𝛽₁) given by (3.5) and (3.6), respectively, and calculate 𝑅₁=𝜃_1,1/𝜃_1,1+𝜃_2,1.(3)Repeat steps (1) and (2)   𝑀 times to get 𝑅₁,…,𝑅_𝑀.(4)Calculate the approximate Bayes estimator and the approximate posterior variance; 𝐸𝑅∣𝑟,𝑠=^𝑀_𝑖=1𝑅_𝑖𝑀=𝑉𝑅∣𝑟𝑅,,𝑠=^𝑀_𝑖=1𝑅^𝑖−𝑅²𝑀.(3.11)

Approximate highest posterior density (HPD) intervals for the stress-strength reliability 𝑅 may be found using the algorithm of Chen and Shao [12].

4. A Simulation Study

We conducted simulations to compare the performance of the point and interval estimators developed in this paper. In our simulations we also included some bootstrap intervals [13], namely, the percentile interval and the bootstrap-t interval based on the logit transformation of 𝑅. We used 𝑛=𝑚=6, 8, 10, 15. We fixed 𝜃₁=1 and used 𝜃₂ = 1, 3, 5, 7, 9. The confidence level taken is (1−𝛼)=0.95. In each simulation run we generated 2000 samples of records from the distributions of 𝑋 and 𝑌. For each pair of samples we calculated the following intervals:(1)ANE: the interval based on the asymptotic normality of the MLE with variance estimate based on the Fisher information matrix, given by (2.14),(2)ANO: the interval based on the asymptotic normality of the MLE with variance estimate based on the observed information matrix, given by (2.15),(3)ANT: the interval based on the asymptotic normality of the MLE of the transformed parameter with variance estimate based on the observed information matrix, given by (2.16),(4)HPD: the approximate highest posterior density interval of Chen and Shao [12],(5)Perc: the percentile interval, [13], and(6)Boot: the bootstrap-t interval based on the logit transformation of the stress-strength reliability parameter, [13].

The biases and mean squared error of the Bayes and the maximum likelihood estimator are simulated and the results are given in Table 1. The lower (L), upper (U), and total (T) error rates and the expected widths (W) of the intervals are approximated using the results of the 2000 simulation replications. For the percentile and bootstrap-t intervals we used 1000 bootstrap resamples. The results of our simulations are given in Table 2.

Table 1 

Simulated biases and mean squared errors of the estimators, 𝜃₁=1.

Table 2 

Simulated error probabilities and expected width of the intervals, 𝜃₁=1.

5. Conclusions

Results concerning the performance of the Bayes estimator and the maximum likelihood estimator are given in Table 1. It appears that the MLE is negatively biased while the Bayes estimator is positively biased. However, the biases are very small, especially for large sample sizes. The mean squared errors of the estimators are very close, especially for larger samples. However, it appears that the MLE has higher MSE for values of 𝑅 near 0.5 and small sample sizes. The situation is reversed for values of 𝑅 close to the extremes.

The simulation results in Table 2 show that the widths of all intervals are maximized when 𝑅=0.5, and they become narrower as the true value of 𝑅 approaches the extremes. As expected, increasing the sample sizes also results in shorter intervals. The performance of the intervals (ANO) and (ANE) based on the asymptotic normality of the MLE are very similar. They tend to be anticonservative and are short. However, for the transformed intervals, the situation completely changed, and the results are very encouraging. This is true for (ANT) interval and the bootstrap-t interval based on the transformed parameter. The performance of bootstrap-t intervals for the original parameter is also investigated, but the results are very poor and not included in Table 2. The Bayes interval (HPD) appears to have better performance than (ANO), and (ANE). However it is dominated by the intervals based on the transformed parameter (ANT) and (Boot). The same is true for the percentile interval which tends to be anticonservative. Concerning the symmetry of lower (L) and upper (U) error rates, it appears that the (ANE), (ANO) and (Perc) intervals are highly asymmetric, especially for extreme values of the stress-strength reliability. On the other hand, the (ANT), (Boot), and (HPD) tend to be symmetric even for small samples. In conclusion, we would recommend the use of the bootstrap-t interval (Boot) based on the transformed parameter followed by (ANT) for all sample sizes.