Interval Estimation for the Stress-Strength Reliability with Bivariate Normal Variables

We propose a procedure to obtain accurate confidence intervals for the stress-strength reliability R = P (X > Y) when (X, Y) is a bivariate normal distribution with unknown means and covariance matrix. Our method is more accurate than standard methods as it possesses a third-order distributional accuracy. Simulations studies are provided to show the performance of the proposed method relative to existing ones in terms of coverage probability and average length. An empirical example is given to illustrate its usefulness in practice.


Introduction
Let ( ) , 1, , i i x y i n =  , be a sample from a bivariate normal distribution with mean ( ) x y µ µ and covariance matrix where , 0 x y σ σ > , and 1 ρ < .The stress-strength reliability of a system where X is the strength and Y is the stress is defined by where ( ) Φ ⋅ is the cumulative distribution function of the standard normal distribution; ( ) , , , , x y x y θ µ µ σ σ ρ ′ = denotes the parameter vector of the model.
Nguimkeu et al. [1] recently proposed a third-order method for inference about R for the case where the normal variables X and Y are independents; that is, 0 ρ = .However, in many empirical applications, the variables of interest are correlated, either directly or through their dependence over a common auxiliary variable.For example, in financial risk-management one may want to compare the stock returns from two companies.If these companies operate in the same industry the prices of these stocks are likely to be correlated.In welfare economics, comparing households' income to households' expenditures can be useful to test households saving capacity or their financial vulnerability.Applying the Nguimkeu et al. [1] test in such context could be misleading.
In this paper, we modify the procedure proposed by Nguimkeu et al. [1] to account for possible correlation between the stress and the strength variables when they are sampled from normal populations.Simulations are used to compute coverage properties of the statistic and compare its performance with existing alternative methods.An empirical example is provided to illustrate the usefulness of the method in practice.

The Procedure
From a sample of observations ( ) x y i n =  the log-likelihood function of , , , , )

∑
The maximum likelihood estimate (MLE) of R can then be obtained by ( ) ˆˆ,

ˆˆˆ2
x y where the MLE of the parameter vector, ( ) From the above log-likelihood function and MLEs, two standard methods for confidence interval estimation of the parameter ψ (and hence R) can be derived: the standardized maximum likelihood estimate method (also known as the Wald method) and the signed log-likelihood ratio method.The Wald method is based on the statistic (q) defined by where the delta method can be applied to estimate the variance of ψ by is the observed information matrix evaluated at θ .With the regularity condition stated in Cox and Hinkley [2] (Chapter 9), q is asymptotically distributed as standard normal, and a ( ) α − confidence interval for ψ can be approximated by where z α is the ( ) percentile of the standard normal.Although the Wald method is simple, it is not invariant to parameterization.
The signed log-likelihood ratio method is based on the statistic (r) defined by The constrained MLE ˆψ θ and the estimated Lagrange multiplier λ  can then be obtained by numerically solving the first-order conditions: ( ) ( ) ) The tilted log-likelihood function is defined by ( ) ( ) ( ) and is the same as the log-likelihood function when evaluated at the constrained MLE of ψ , i.e. ( ) ( ) . The observed information matrix of the tilted log-likelihood function evaluated at ˆψ θ , denoted ( ) Again, with the regularity conditions stated in Cox and Hinkley [2] (Chapter 9), ( ) r ψ is asymptotically dis- tributed as standard normal.Hence, a ( ) It is well known that both the Wald method and the signed log-likelihood ratio method are first-order methods; that is, both ( ) q ψ and ( ) r ψ converge in distribution to the standard normal distribution with rate of conver- gence ( ) O n − .Note that, computationally, confidence intervals for ψ can easily be obtained from ( 9) but the methodology is not invariant to reparameterization.While confidence intervals for ψ obtained from (15) generally require the use of numerical methods, the method is parameterization invariant.Doganaksoy and Schmee [3] showed that confidence intervals obtained from (15) have better coverage properties than those obtained from (9).
To improve the accuracy of the first-order methods, Barndorff-Nielsen [4] [5] introduced the modified signed log-likelihood ratio statistic where ( ) r ψ is the signed log-likelihood ratio statistic defined in (10), and ( ) Q ψ is a statistic that is based on the log-likelihood function and an ancillary statistic.Barndorff-Nielsen [4] [5] showed that ( ) * r ψ is asymp- totically distributed as standard normal with third-order accuracy.Fraser and Reid [6] showed that for the exponential family model,

( )
Q ψ is the standardized maximum likelihood estimate calculated in the canonical pa- rameter scale.Reid [7] and Severeni [8] provide a detailed overview of this development.
The stress-strength reliability with dependent normal random variables correspond to an exponential family model with canonical parameter given by ( ) which denote the derivatives of ( ) ϕ θ and ( ) ψ θ with respect to θ .The ( ) For the parameter of interest ( ) By change-of-basis then, ( ) ( ) calculated in the scale of the canonical parameter, ( ) The standardized MLE of ψˆ calculated in the ( ) The modified signed log-likelihood ratio statistic can then be obtained from (16).It is asymptotically distributed as standard normal with a ( ) O n − distributional accuracy.The ( ) Practically, this confidence interval is obtained by numerically solving for ψ in the inequality ( ) using a sufficient number of grid points of ψ chosen in an appropriate range.It follows that the ( ) α − confidence interval for R is then given by ( ) , where ( ) Φ ⋅ is the cumulative density function of the standard normal distribution.

Numerical Studies
An empirical example and Monte Carlo simulation studies are considered in this section.The aim of the empirical example is to illustrate how the various methods considered in this paper can produce quite different confidence intervals for ψ .Simulation studies are then performed to examine the statistical properties of the pro- posed method in terms of central coverage and average confidence interval length at the nominal size of 95%.To illustrate the accuracy of our proposed method (Proposed), we compare it with the commonly used asymptotic methods, i.e., the signed log-likelihood ratio test (r) and the Wald test (Q).We also compare it with approximations that were recently discussed in Barbiero [9].Specifically, Barbiero [9] proposed results based on approximate confidence intervals from the asymptotic variance of ψ (denoted AN2), approximate confidence intervals based on a logit transformation of R (denoted LOGIT) and a bootstrap bias-corrected and accelerated percentile confidence interval for R (denoted BCAPB).Since AN2 and LOGIT yield very similar results, only the AN2 and BCAPB procedures of Barbiero [9] are reported here.

An Example
Table 1 is a data set from Azen and Reed [10] which shows absorbance values of substances analyzed in 19 runs of a laboratory test for serum concentration level of an enzyme, leucine amino peptidase.Controls 1 and 2 are two control samples from the same pool analyzed in the same run.We estimate the confidence interval of the probability that Control 1 (X) is more than Control 2 (Y), assuming a bivariate normal distribution for the data.
For these data, .Using the methods developed in the previous sections, the 90%, 95% and 99% CIs for R are computed and reported in Table 2.The confidence intervals obtained from the five methods are quite different.Hence, simulation studies are performed to examine the accuracy of the five methods.

Design of the Monte Carlo Simulation Studies
The simulation set up is similar to Barbiero [9] [11].An array of eight different scenarios has been considered, each corresponding to a different combination of distribution parameters (and thus different reliabilities).These scenarios have been coded with a progressive number which is reported in Table 3.Without any loss in generality, the mean and standard deviation for Y has been set to 0 y µ = and 2 1 y σ = while the parameters x µ , 2 x σ and ρ vary.The correlation coefficient ρ takes two values, 0.5 and 0.8.The five parameters have been jointly set in order to assure values higher than 0.5 for the reliability in effort to reflect practice where there is concern for high reliability for the study component.The analyzed scenarios however cover a large range of reliability, since R goes from 0.614 to 0.943.Different sample sizes (n = 10, 20, 30, 50) are used in order to examine the reliance of the procedure on small samples.The number of Monte Carlo replications has been fixed at N = 10,000.The simulation study for the comparison of interval estimators proceeds as follows: 1) Set the parameters , , , , for the bivariate random variable ( ) , X Y and compute the corresponding R from Equation (1) (see Table 3); 2) Draw a random sample of size n from ( ) , X Y ; 3) Estimate R and a CI for R, using each listed procedures; 4) Check if this CI contains R; compute its length; 5) Repeat the two precedent steps N = 10,000 times and compute the overall CI coverage (proportion of the CIs containing R) and average length (computed over the 10,000 replications) for each interval estimator.

Results of the Monte Carlo Simulation Studies
The results of the simulation studies are reported in Table 4 and Table 5.Table 4 records the coverage probability produced by each method and Table 5 records the average CI length produced by each method.The accuracy of the proposed method is striking.Inspecting Table 4 reveals how accurate the coverage probability of the proposed method is across all sample sizes and scenarios.Figure 1 presents a graphical illustration of the coverage probabilities of each method graphed against the various scenarios.This figure provides visual confirmation of the accuracy of the proposed method.These findings are consistent with the fact that the proposed method has theoretically a third-order distributional accuracy.Close to our method in terms of accuracy is the bias-corrected and accelerated percentile bootstrap (BCAPB).The approximate estimator (AN2) also performs of the difference of the two variables; and

Figure 1 .
Figure 1.Monte Carlo coverage for various methods.

Table 1 .
Sample data for the empirical example.

Table 2 .
Confidence intervals and lengths for the example.

Table 3 .
Parameter values for the Monte Carlo simulation.

Table 5 .
Monte Carlo simulation results: average CI length.