The following R code may be used for constructing bias-corrected percentile confidence intervals for stress-strength models. These bias-corrected percentile confidence intervals are obtained with a bootstrapping method.

The code focuses at first on a situation in which stress and strength both follow a normal distribution. However, any other combination of distributions may be used in the R code (e.g, having stress follow a lognormal distribution, while strength a Weibull distribution). In fact, it will be demonstrated how to change the code such that the stress and strength both follow a Weibull distribution.

A comprehensive account of bootstrap confidence intervals for stress-strength models is given in Barbiero, A. (2011), Confidence Intervals for Reliability of Stress-Strength Models in the Normal Case, *Communications in Statistics – Simulation and Computation*, 40, 907-925. Note though that the bootstrap confidence intervals described/reported by Alessandro Barbiero in his 2011 paper are standard (or simple) percentile confidence intervals and * not* bias-corrected percentile confidence intervals.

Alessandro Barbiero also maintains an R package called StressStrength. However, the package is limited to situations in which stress and strength are both normally distributed.

See this blog post for computing likelihood based confidence intervals for stress-strength models.

A nice introduction to stress-strength models is given by Kapur and Lamberson in their 1977 book *Reliability in engineering design*.

Suggestions and/or questions? Please contact Stefan Gelissen (email: info at datall-analyse.nl).

See this page for an overview of all of Stefan’s *R code* blog posts.

**R code**

```
library(survival)
library(SPREDA)
library(StressStrength)
##stress and strength both follow a normal distribution
#generate artificial data
#strength data
strength <- rnorm(20, 6000, 600)
#stress data
stress <- rnorm(20, 5000, 190)
#fit stress and strength distribution (note: both follow a normal distribution)
sts <- survreg(Surv(stress)~1, dist="gaussian")
str <- survreg(Surv(strength)~1, dist="gaussian")
#numerical integration
integrand <- function(x, mu1, mu2, s1, s2) {
#pdf for stress
f1 <- dnorm((x-mu1)/s1) / s1
#cdf for strength
f2 <- pnorm((x-mu2)/s2)
f1*f2
}
#determine upper limit of the integral
uplimit <- qnorm(.99999999)*sts$scale+coef(sts)
#estimate of unreliability
Uhat <- integrate(integrand, 0, uplimit, mu1=coef(sts), mu2=coef(str),
s1=sts$scale, s2=str$scale)$value
#estimate of reliability (=Rhat)
(Rhat <- 1-Uhat)
#compare the numerical solution of the integral with the exact solution
SSR(c(coef(str), str$scale), c(coef(sts), sts$scale))
#confidence interval for reliability estimate
#initialize bootstrap
nboot <- 4000 #number of bootstrap samples
nobsStress <- length(stress)
nobsStrength <- length(strength)
simVec <- vector(mode="numeric", length=nboot)
#perform bootstrap
for(j in 1:nboot){
#draw bootstrap samples that mimic the original samples
Lboot <- coef(sts) + sts$scale*rnorm(nobsStress)
Sboot <- coef(str) + str$scale*rnorm(nobsStrength)
#fit stress and strength distribution to bootstrap samples
stsboot <- survreg(Surv(Lboot)~1, dist="gaussian")
strboot <- survreg(Surv(Sboot)~1, dist="gaussian")
#numerical integration
uplimit <- qnorm(.99999999)*stsboot$scale+coef(stsboot)
simVec[j] <- integrate(integrand, 0, uplimit,
mu1=coef(stsboot), mu2=coef(strboot),
s1=stsboot$scale, s2=strboot$scale)$value
}
#reliability estimates for bootstrap samples
Rhatboot <- 1-simVec
#bias-corrected percentile confidence interval
#calculate proportion of Rhatboot values that are less than Rhat
(q0 <- mean(Rhatboot <= Rhat))
#two sided confidence interval
alpha <- .05
zl.bias2 <- round(nboot*pnorm(2*qnorm(q0)+qnorm(alpha/2)))
zu.bias2 <- round(nboot*pnorm(2*qnorm(q0)+qnorm(1-alpha/2)))
c(sort(Rhatboot)[zl.bias2], sort(Rhatboot)[zu.bias2])
#confidence interval by Barbiero's bootstrap method
#(this is the bootstrap method described/applied in Barbiero [2011])
estSSR(strength, stress, type="B", alpha=.05, B=4000)$CI
#confidence interval by the Reiser-Guttman method
estSSR(strength, stress, type="RG", alpha=.05)$CI
#Barbiero (2011) showed that the Reiser-Gutmann method yields on average good results
#that is, the coverage probability of these intervals is often close to the nominal level
#note that the bias-corrected percentile interval lies closer to the Reiser-Guttman
#interval than the interval obtained with Barbiero's bootstrap method
#one sided confidence interval
alpha <- .05
zl.bias1 <- round(nboot*pnorm(2*qnorm(q0)+qnorm(alpha)))
sort(Rhatboot)[zl.bias1]
#confidence interval by Barbiero's bootstrap method
#(this is the bootstrap method described/applied in Barbiero [2011])
estSSR(strength, stress, type="B", alpha=.05, B=4000, twoside=FALSE)$CI
#confidence interval by the Guo-Krishnamoorthy method
estSSR(strength, stress, type="GK", alpha=.05, twoside=FALSE)$CI
#Barbiero (2011) demonstrated that the Guo-Krishnamoorthy method yields
#on average excellent results
#that is, the coverage probability of these intervals is often very close to the nominal level
#note that the bias-corrected percentile interval lies closer to the Guo-Krishnamoorthy
#interval than the interval obtained with Barbiero's bootstrap method
##stress and strength both follow a Weibull distribution
#generate artificial data
#strength data
strength <- exp(8 + .6*rsev(30))
#stress data
stress <- exp(7 + .3*rsev(40))
#fit stress and strength distribution (note: both follow a Weibull distribution)
sts <- survreg(Surv(stress)~1, dist="weibull")
str <- survreg(Surv(strength)~1, dist="weibull")
#numerical integration
integrand <- function(x, mu1, mu2, s1, s2) {
#pdf for stress
f1 <- dsev((log(x)-mu1)/s1) / (s1*x)
#cdf for strength
f2 <- psev((log(x)-mu2)/s2)
f1*f2
}
#determine upper limit of the integral
uplimit <- exp(qsev(.99999999)*sts$scale+coef(sts))
#estimate of unreliability
Uhat <- integrate(integrand, 0, uplimit, mu1=coef(sts), mu2=coef(str),
s1=sts$scale, s2=str$scale)$value
#estimate of reliability (=Rhat)
(Rhat <- 1-Uhat)
#compare this reliability estimate with that obtained by a Monte Carlo simulation:
#draw samples (of size ns) from the fitted stress and strength
#distributions, and subsequently calculate the proportion where the
#stress (L) is larger than the strength (S)
ns <- 100000
L <- exp(coef(sts) + sts$scale*rsev(ns))
S <- exp(coef(str) + str$scale*rsev(ns))
n <- ifelse(S > L, TRUE, FALSE)
#reliability
mean(n)
#confidence interval for reliability estimate
#initialize bootstrap
nboot <- 4000 #number of bootstrap samples
nobsStress <- length(stress)
nobsStrength <- length(strength)
simVec <- vector(mode="numeric", length=nboot)
#perform bootstrap
for(j in 1:nboot){
#draw bootstrap samples that mimic the original samples
Lboot <- exp(coef(sts) + sts$scale*rsev(nobsStress))
Sboot <- exp(coef(str) + str$scale*rsev(nobsStrength))
#fit stress and strength distribution to bootstrap samples
stsboot <- survreg(Surv(Lboot)~1, dist="weibull")
strboot <- survreg(Surv(Sboot)~1, dist="weibull")
#numerical integration
uplimit <- exp(qsev(.99999999)*stsboot$scale+coef(stsboot))
simVec[j] <- integrate(integrand, 0, uplimit,
mu1=coef(stsboot), mu2=coef(strboot),
s1=stsboot$scale, s2=strboot$scale)$value
}
#reliability estimates for bootstrap samples
Rhatboot <- 1-simVec
#bias-corrected percentile confidence interval
#calculate proportion of Rhatboot values that are less than Rhat
(q0 <- mean(Rhatboot <= Rhat))
#two sided confidence interval for reliability
alpha <- .05
zl.bias2 <- round(nboot*pnorm(2*qnorm(q0)+qnorm(alpha/2)))
zu.bias2 <- round(nboot*pnorm(2*qnorm(q0)+qnorm(1-alpha/2)))
c(sort(Rhatboot)[zl.bias2], sort(Rhatboot)[zu.bias2])
#one sided confidence interval for reliability
alpha <- .05
zl.bias1 <- round(nboot*pnorm(2*qnorm(q0)+qnorm(alpha)))
sort(Rhatboot)[zl.bias1]
```