# Inference onPr(X<Y) in the Two-Parameter Weibull Model Based on Records

Authors: Ayman Baklizi
Journal: ISRN Probability and Statistics
 (2012)

---

## 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.


---

## Body

## 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 W(θ,β) are given by
(1.1)f(x)=θβxβ-1e-θxβ,x>0,θ>0,β>0F(x)=1-e-θxβ.LetX1,X2,… be an infinite sequence of iid random variables. An observation Xj is called a record if its value is greater than all previous observations, that is Xj>Xi for every i<j. We want to estimate the stress-strength reliability Pr(X<Y) using the data on records. The stress-strength reliability Pr(X<Y) arises in life-testing experiments when X and Y represent the lifetimes of two devices, and it gives the probability that the device with life time  X  fails before the other. As an example, Hall [6] studied a situation where the breakdown voltage Y of a capacitor must exceed the voltage output  X 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. LetX~W(θ1,β) and Y~W(θ2,β) be independent random variables. Let R=Pr(X<Y) be the stress strength reliability. Let  r_=r0,…,rn be a set of records from W(θ1,β), and let  s_=s0,…,sm be an independent set of records from   W(θ2,β). The likelihood functions are given by [3] as follows:
(2.1)L1(β,θ1∣r_)=f(rn)∏i=0n-1f(ri)(1-F(ri)),L2(β,θ2∣s_)=g(sm)∏i=0m-1g(si)(1-G(si)),
where  f and  F  are the pdf and cdf of  X~W(θ1,β)  and  g and  G  are the pdf and cdf of  Y~W(θ2,β). The likelihood function of  (β,θ1,θ2)  based on  r_,s_  is given by
(2.2)L(β,θ1,θ2∣r_,s_)=(θ1β)n+1e-θ1rnβ∏i=0nriβ-1(θ2β)m+1e-θ2smβ∏i=0msiβ-1.Taking the natural logarithm we get the log-likelihood function(2.3)l(β,θ1,θ2)=(n+m+2)ln⁡β+(n+1)ln⁡θ1+(m+1)ln⁡θ2-θ1rnβ-θ2smβ+(β-1)[∑i=1nln⁡ri+∑i=1mln⁡si].The first partial derivatives are given by(2.4)∂l(β,θ1,θ2)∂β=(n+m+2)β-θ1rnβln⁡rn-θ2smβln⁡sm+[∑i=1nln⁡ri+∑i=1mln⁡si],(2.5)∂l(β,θ1,θ2)∂θ1=(n+1)θ1-rnβ,∂l(β,θ1,θ2)∂θ2=(m+1)θ2-smβ.Equating these partial derivatives to zero and solving simultaneously we obtain(2.6)β^=(n+m+2)(n+1)ln⁡rn+(m+1)ln⁡sm-[∑i=0nln⁡ri+∑i=0mln⁡si],θ^1=n+1rnβ^,θ^2=m+1smβ^.Hence the MLE ofR is given by  R^=θ^2/θ^1+θ^2. The study of the exact distribution of  R^  is apparently rather complicated, so we will consider the asymptotic distribution. We need the asymptotic joint distribution of  β^, θ^1  and  θ^2.  The second partial derivatives of the log-likelihood function are given by
(2.7)∂2l(β,θ1,θ2)∂β2=-(n+m+2)β2-θ1rnβ(ln⁡rn)2-θ2smβ(ln⁡⁡sm)2,∂2l(β,θ1,θ2)∂θ12=-(n+1)θ12,∂2l(β,θ1,θ2)∂θ22=-(m+1)θ22,∂2l(β,θ1,θ2)∂β∂θ1=-rnβln⁡⁡rn,∂2l(β,θ1,θ2)∂β∂θ2=-smβln⁡⁡sm,∂2l(β,θ1,θ2)∂θ1∂θ2=0.The Fisher information matrix is given by(2.8)I=(I11I12I13I21I22I23I31I32I33)=-(E(∂2l∂β2)E(∂2l∂β∂θ1)E(∂2l∂β∂θ2)E(∂2l∂β∂θ1)E(∂2l∂θ12)E(∂2l∂θ1∂θ2)E(∂2l∂β∂θ2)E(∂2l∂θ1∂θ2)E(∂2l∂θ22)).Now we will find the entries of the information matrix. We need to findE(rnβlnrn), E(rnβ(lnrn)2), E(smβlnsm) and E(smβ(lnsm)2). Using result (2.4.3) of Arnold et al. [3] we have U=θ1rnβ~Γ(n+1,1),  and it follows that(2.9)E(rnβln⁡rn)=1βθ1E(θ1rnβ(ln⁡θ1rnβ-ln⁡θ1))=1βθ1[E(Uln⁡⁡U)-ln⁡⁡θ1E(U)]=(n+1)βθ1[ψ(n+2)-ln⁡θ1],(2.10)E(rnβ(ln⁡rn)2)=1θ1β2E(θ1rnβ(ln⁡θ1rnβ-ln⁡θ1)2)=1θ1β2E(θ1rnβ((ln⁡θ1rnβ)2-2ln⁡θ1lnθ1rnβ+(ln⁡θ1)2))=1θ1β2[E(U(ln⁡U)2)-2ln⁡θ1E(Uln⁡U)+(ln⁡θ1)2E(U)]=(n+1)θ1β2[(ψ(n+2)2+ψ′(n+2))-2ln⁡θ1ψ(n+2)+(ln⁡θ1)2].Similarly, sinceV=θ2smβ~Γ(m+1,1),  we have
(2.11)E(smβln⁡sm)=(m+1)βθ2[ψ(m+2)-ln⁡θ2],E(smβ(ln⁡sm)2)=(m+1)θ2β2[(ψ(m+2)2+ψ′(m+2))-2ln⁡θ2ψ(m+2)+(ln⁡θ2)2].LetR=θ1/(θ1+θ2)=h(θ1,θ2),  the maximum likelihood estimator of R is given by  R^=θ^1/(θ^1+θ^2)=h(θ^1,θ^2), and it follows that
(2.12)n(R^-R)=n[h(θ^1,θ^2)-h(θ1,θ2)]→DN(0,η2),
as n→∞, m→∞ such that m/n→p, 0<p<1. The asymptotic variance is given by
(2.13)η2(∂h(θ1,θ2)∂θ1)2var(θ^1)+(∂h(θ1,θ2)∂θ2)2var(θ^2)p+2cov⁡(θ^1,θ^2)(∂h(θ1,θ2)∂θ1)(∂h(θ1,θ2)∂θ2).A(1-α)% confidence interval for R based on this asymptotic result is given by
(2.14){R^-z1-α/2η^,R^+z1-α/2η^},
where η^ is obtained by substituting  m/n  for p and the MLEs of θ1 and  θ2 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  R^. The new interval is given by
(2.15){R^-z1-α/2η^^,R^+z1-α/2η^^}.
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(R/1-R), and the maximum likelihood of  λ  is given by λ^=ln(R^/1-R^). The asymptotic variance of  n(λ^-λ)  is given by  η2(dλ/dR)2. This variance can be estimated by substituting the maximum likelihood estimator instead of the parameter. A (1-α)% confidence interval for λ is given by (λ^L,λ^U), where λ^L=λ^-z1-α/2η^dλ/dR∣R=R^  and  λ^U=λ^+z1-α/2η^dλ/dR∣R=R^,  z1-α/2 is the (1-α/2) quantile of the standard normal distribution. Therefore A (1-α)% confidence interval for R is given by
(2.16)(eλ^L1+eλ^L,eλ^U1+eλ^U).


## 3. Bayesian Inference

As in Kundu and Gupta [11], we will use the conjugate gamma priors for the scale parameters θ1 and θ2 and a squared error loss function. The prior density of  θj therefore is given by
(3.1)πj(θj∣aj,bj)=bjajΓ(aj)θjaj-1e-bjθj,j=1,2,
where aj>0, bj>0 are the parameters of the prior distributions. Assuming that θ1 and θ2 are independent, the joint prior distribution of  θ1 and θ2 is given by  π(θ1,θ2)=∏j=12πj(θj). Assume that the prior distribution of  β,  denoted by ζ(β),  has a support on (0,∞). Recall that the likelihood function of  (β,θ1,θ2) based on  r_=r0,…,rn  and  s_=s0,…,sm is given by
(3.2)L(β,θ1,θ2∣r_,s_)=(θ1β)n+1e-θ1rnβ∏i=0nriβ-1(θ2β)m+1e-θ2smβ∏i=0msiβ-1.The joint posterior density function of(β,θ1,θ2) therefore is given by;
(3.3)π*(β,θ1,θ2)=π(θ1,θ2)ζ(β)L(β,θ1,θ2∣r_,s_)∭0∞π(θ1,θ2)ζ(β)L(β,θ1,θ2∣r_,s_)dβdθ1dθ2,
where the numerator is given by
(3.4)βn+m+2ζ(β)b1a1Γ(a1)θ1n+a1e-θ1(rnβ+b1)∏i=0nriβ-1b2a2Γ(a2)θ2m+a2e-θ2(smβ+b2)∏i=0msiβ-1.This expression forπ*(β,θ1,θ2) is difficult or even impossible to find in closed form. A simulation technique is needed. It is clear that the conditional posterior distributions of  θ1  and θ2 given β are given by
(3.5)π1*(θ1∣β)=(rnβ+b1)n+a1+1Γ(n+a1+1)θ1n+a1e-θ1(rnβ+b1),(3.6)π2*(θ2∣β)=(smβ+b2)m+a2+1Γ(m+a2+1)θ2m+a2e-θ2(smβ+b2).Note that givenβ,  θ1, and θ2  are independent. From the joint posterior distribution of θ1, θ2, and β  we obtain
(3.7)ζ*(β)∝ζ(β)βn+m+2(∏i=0nri∏i=0msi)β(rnβ+b1)-(n+a1+1)(smβ+b2)-(m+a2+1).Takinga1=a2=0 and  b1=b2=0 in the prior distributions we obtain
(3.8)ζ*(β)∝ζ(β)βn+m+2(rn-(n+1)sm-(m+1)∏i=0nri∏i=0msi)β.If we can takeζ(β)=e-β, we obtain
(3.9)ζ*(β)∝βn+m+2(e-1rn-(n+1)sm-(m+1)∏i=0nri∏i=0msi)β.It follows that the posterior distribution ofβ is
(3.10)Γ(n+m+3,ln⁡(e-1rn-(n+1)sm-(m+1)∏i=0nri∏i=0msi)-1).We can use the following Monte Carlo method to find approximate point estimates.(1)
Generateβ1 from ζ*(β)  given by (3.10).(2)
Generateθ1,1 and  θ2,1  from  π1*(θ1∣β1)  and  π2*(θ2∣β1) given by (3.5) and (3.6), respectively, and calculate  R1=θ1,1/θ1,1+θ2,1.(3)
Repeat steps(1) and (2)   M times to get  R1,…,RM.(4)
Calculate the approximate Bayes estimator and the approximate posterior variance;(3.11)E^(R∣r_,s_)=∑i=1MRiM=R-,V^(R∣r_,s_)=∑i=1M(Ri-R-)2M.Approximate highest posterior density (HPD) intervals for the stress-strength reliabilityR 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 R. We used n=m=6, 8, 10, 15. We fixed θ1=1  and used θ2 = 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  X and Y. 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 Table1. 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=1.
BiasR^
BiasR~
MSER^
MSER~
BiasR^
BiasR~
MSER^
MSER~
θ
2
(
n
,
m
)
=
(
6,6
)
(
n
,
m
)
=
(
8,8
)
1
−0.004
−0.004
0.026
0.020
0.000
0.000
0.018
0.016
3
−0.020
0.011
0.015
0.012
−0.019
0.003
0.011
0.009
5
−0.021
0.027
0.010
0.008
−0.017
0.015
0.008
0.006
7
−0.018
0.029
0.006
0.007
−0.013
0.018
0.005
0.005
(
n
,
m
)
=
(
10,10
)
(
n
,
m
)
=
(
15,15
)
9
−0.017
0.031
0.004
0.006
−0.013
0.022
0.004
0.004
1
0.001
0.000
0.014
0.013
0.000
−0.000
0.009
0.009
3
−0.012
0.004
0.010
0.008
−0.009
−0.003
0.006
0.006
5
−0.015
0.007
0.006
0.005
−0.009
0.002
0.004
0.004
7
−0.013
0.012
0.004
0.004
−0.010
0.003
0.003
0.003
9
−0.012
0.014
0.003
0.003
−0.008
0.005
0.002
0.002Table 2
Simulated error probabilities and expected width of the intervals,θ1=1.
Interval
L
U
T
W
L
U
T
W
θ
2
(n, m) = (6,6)
(n, m) = (8,8)
1
ANE
0.079
0.080
0.159
0.492
0.065
0.062
0.127
0.443
1
ANO
0.079
0.080
0.159
0.492
0.065
0.062
0.127
0.443
1
ANT
0.043
0.046
0.088
0.460
0.040
0.037
0.076
0.418
1
HPD
0.027
0.065
0.092
0.545
0.023
0.068
0.091
0.507
1
Perc
0.063
0.063
0.126
0.562
0.053
0.048
0.101
0.498
1
BOOT
0.025
0.031
0.056
0.506
0.026
0.026
0.052
0.456
3
ANE
0.170
0.017
0.187
0.398
0.135
0.015
0.150
0.359
3
ANO
0.170
0.017
0.187
0.398
0.135
0.015
0.150
0.359
3
ANT
0.039
0.027
0.065
0.398
0.048
0.025
0.073
0.358
3
HPD
0.057
0.050
0.107
0.434
0.052
0.043
0.094
0.395
3
Perc
0.164
0.009
0.172
0.408
0.140
0.005
0.145
0.366
3
BOOT
0.029
0.015
0.043
0.435
0.036
0.015
0.051
0.387
5
ANE
0.211
0.009
0.220
0.313
0.186
0.008
0.194
0.287
5
ANO
0.211
0.009
0.220
0.313
0.186
0.008
0.194
0.287
5
ANT
0.051
0.023
0.074
0.339
0.049
0.025
0.073
0.304
5
HPD
0.037
0.042
0.079
0.364
0.051
0.037
0.088
0.322
5
Perc
0.200
0.004
0.204
0.296
0.178
0.004
0.181
0.273
5
BOOT
0.035
0.018
0.053
0.370
0.031
0.020
0.050
0.329
7
ANE
0.224
0.005
0.229
0.263
0.200
0.003
0.203
0.243
7
ANO
0.224
0.005
0.229
0.263
0.200
0.003
0.203
0.243
7
ANT
0.044
0.017
0.060
0.302
0.047
0.022
0.069
0.270
7
HPD
0.044
0.036
0.079
0.312
0.055
0.037
0.092
0.274
7
Perc
0.209
0.002
0.211
0.237
0.187
0.001
0.188
0.223
7
BOOT
0.033
0.014
0.047
0.331
0.030
0.020
0.050
0.293
9
ANE
0.261
0.003
0.264
0.222
0.212
0.004
0.215
0.208
9
ANO
0.261
0.003
0.264
0.222
0.212
0.004
0.215
0.208
9
ANT
0.049
0.013
0.062
0.271
0.035
0.017
0.051
0.242
9
HPD
0.038
0.036
0.074
0.278
0.048
0.033
0.080
0.244
9
Perc
0.237
0.001
0.237
0.195
0.188
0.003
0.191
0.186
9
BOOT
0.031
0.010
0.041
0.298
0.023
0.017
0.040
0.264
θ
2
(n, m) = (10,10)
(n, m) = (15,15)
1
ANE
0.053
0.059
0.111
0.404
0.048
0.049
0.096
0.339
1
ANO
0.050
0.058
0.108
0.404
0.048
0.049
0.096
0.340
1
ANT
0.032
0.042
0.074
0.385
0.038
0.036
0.074
0.328
1
HPD
0.015
0.068
0.083
0.475
0.009
0.066
0.074
0.415
1
Perc
0.043
0.053
0.095
0.449
0.045
0.041
0.086
0.368
1
BOOT
0.020
0.028
0.047
0.415
0.021
0.027
0.048
0.348
3
ANE
0.119
0.020
0.139
0.330
0.097
0.015
0.112
0.281
3
ANO
0.113
0.020
0.133
0.333
0.095
0.015
0.110
0.283
3
ANT
0.041
0.031
0.072
0.331
0.041
0.022
0.063
0.281
3
HPD
0.044
0.049
0.093
0.370
0.028
0.036
0.063
0.321
3
Perc
0.115
0.015
0.130
0.339
0.096
0.009
0.105
0.286
3
BOOT
0.028
0.024
0.051
0.355
0.023
0.020
0.043
0.296
5
ANE
0.152
0.007
0.159
0.263
0.125
0.009
0.134
0.228
5
ANO
0.147
0.006
0.153
0.267
0.120
0.009
0.129
0.231
5
ANT
0.045
0.019
0.063
0.279
0.038
0.022
0.060
0.237
5
HPD
0.055
0.034
0.089
0.294
0.040
0.034
0.074
0.255
5
Perc
0.143
0.003
0.145
0.255
0.116
0.007
0.123
0.223
5
BOOT
0.031
0.016
0.047
0.300
0.025
0.022
0.047
0.250
7
ANE
0.168
0.004
0.172
0.220
0.147
0.006
0.152
0.190
7
ANO
0.162
0.003
0.165
0.224
0.142
0.006
0.148
0.193
7
ANT
0.042
0.021
0.062
0.244
0.040
0.022
0.061
0.204
7
HPD
0.048
0.032
0.080
0.250
0.048
0.032
0.080
0.211
7
Perc
0.147
0.001
0.148
0.207
0.135
0.004
0.139
0.182
7
BOOT
0.027
0.022
0.049
0.263
0.026
0.029
0.055
0.216
9
ANE
0.201
0.004
0.204
0.188
0.143
0.003
0.146
0.166
9
ANO
0.198
0.003
0.201
0.192
0.139
0.003
0.141
0.168
9
ANT
0.039
0.016
0.055
0.218
0.033
0.020
0.053
0.183
9
HPD
0.051
0.031
0.081
0.219
0.053
0.023
0.075
0.183
9
Perc
0.176
0.002
0.178
0.174
0.124
0.002
0.126
0.156
9
BOOT
0.026
0.017
0.043
0.235
0.025
0.027
0.052
0.194
L: lower error probability, U: upper error probability, T: total error probability, W: expected width of the interval.


## 5. Conclusions

Results concerning the performance of the Bayes estimator and the maximum likelihood estimator are given in Table1. 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 R near 0.5 and small sample sizes. The situation is reversed for values of  R close to the extremes.The simulation results in Table2 show that the widths of all intervals are maximized when  R=0.5, and they become narrower as the true value of R 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.



---

*Source: 263612-2012-09-12.xml*