# R code for constructing likelihood based confidence intervals for stress-strength models

In my previous blog post I demonstrated how to compute bootstrap confidence intervals for stress-strength models.
The following R-code may be used for computing likelihood based confidence intervals for stress-strength models.

Definition of the unreliability in case of a stress-strength model
A component may fail when the stress (or load) exceeds the strength. Accordingly, the unreliability of the component is defined as: $U_{hat} =\int^{+\infty}_{-\infty} f_l(l)F_s(l) \ dl$

where Uhat is the unreliability of the component, fl(·) the probability density function of the load, and Fs(·) the cumulative density function of the strength.

Note that both the load and strength random variables never take negative values. Since both these random variables are non-negative, the lower limit of the above integral changes. That is, the lower limit changes from -∞ to 0, since fl(·) and Fs(·) will be 0 for negative values (i.e., for l<0). As a consequence, the above integral is given by: $U_{hat} =\int^{+\infty}_{0} f_l(l)F_s(l) \ dl$

The reliability of the component is Rhat=1-Uhat.

Computing likelihood based intervals for Uhat
For computing the likelihood, a location-scale parameterization for the load and strength distributions will be employed. The location of the distribution for the load random variable will be denoted by μl, and its scale by σl. In addition, μs and σs are, respectively, the location and scale of strength random variable.
A consequence of using a location-scale based distribution is that fl(·) in the above integral becomes a function of μl and σl. Similarly, Fs(·) will be some function of μs and σs.
Furthermore, the likelihood function for the load (denoted by Ll) can now be defined in terms of μl and σl, and the likelihood function for the strength (denoted by Ls) in terms of μs and σs.
Finally, under the assumption that the load and strength random variables are independent, the likelihood function for stress-strength model is given by Lss=Ll(μl, σl)*Ls(μs, σs).

The likelihood function Lss will be used in computing the likelihood based intervals for Uhat. However, for computing these Uhat confidence intervals we need to reparameterize Lss. For this reparameterization to work, remember that Uhat was defined as some function of the parameters μl, σl, μs, and σs , denoted by Uhat=f(μl, σl, μs, σs). As a consequence, it is also possible to express μl as some function of Uhat, σl, μs, and σs , denoted by μl=g(Uhat, σl, μs, σs).

The reparameterization of Lss requires that μl is substituted by g(Uhat, σl, μs, σs) in the original likelihood function for Lss. In that way, Lss becomes a function of Uhat, σl, μs, σs. This latter parameterization of Lss makes it possible to construct likelihood based confidence intervals for Uhat.

For our specific reparameterization of Lss, it is thus needed to express μl as a function of the parameters Uhat, σl, μs, and σs. In other words, we need to compute μl given any values for Uhat, σl, μs, and σs.
The definition of Uhat showed that Uhat=f(μl, σl, μs, σs), and therefore that Uhat f(μl, σl, μs, σs) = 0. This latter equation may yield our sought for value for μl, given any values for the remaining parameters in the equation.
One way of obtaining this value for μl is by relying on numerical optimization. More specifically, μl can be found for any given values Uhat, σl, μs, and σs by minimizing abs(Uhat – f(μl, σl, μsσs)) over all possible values of μl. The following R code focuses at first on a situation in which the stress (or load) and strength both follow a normal distribution. However, any other combination of distributions may be used in the R code (e.g, having the load 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.

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

##data

#generate artificial data

#stress data
stress <- rnorm(20, 5000, 190)

#strength data
strength <- rnorm(20, 6000, 600)

#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")

#obtain MLE for mu and sigma
#Lmu
stsMu <- coef(sts)
names(stsMu) <- "stsMu"
#Lsigma
stsSigma <- sts$scale names(stsSigma) <- "stsSigma" #Smu strMu <- coef(str) names(strMu) <- "strMu" #Ssigma strSigma <- str$scale
names(strSigma) <- "strSigma"

##compute likelihood based interval

#specify integrand for 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
}

#estimate of unreliability using numerical integration
uplimit <- qnorm(.99999999)*sts$scale+coef(sts) #upper limit of the integral 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))

#function for computing abs(Uhat-f(stsMu,stsSigma,strMu,stsSigma))
stsMuMin <- function (stsMuPar, stsSigma, strMu, strSigma, Uhat){
uplimit <- qnorm(.99999999)*stsSigma+stsMuPar
abs(Uhat-integrate(integrand, 0, uplimit, mu1=stsMuPar, mu2=strMu,
s1=stsSigma, s2=strSigma)$value)} #log-likelihood function stsMuSE <- sqrt(diag(vcov(sts))) #extract standard error for stsMu llik.ss <- function (modpar, stsObs, strObs, Uhat) { stssigma <- exp(modpar) strsigma <- exp(modpar) strmu <- modpar stsmu <- optim(par=stsMu, stsMuMin, method="Brent", lower=stsMu-20*stsMuSE, upper=stsMu+20*stsMuSE, Uhat=Uhat, stsSigma=stssigma, strMu=strmu, strSigma=strsigma)$par

z1 <- (stsObs-stsmu)/stssigma
z2 <- (strObs-strmu)/strsigma

#assuming that the stress and strength random variables are independent,
#the log-likelihood function is given by:
sum(log(dnorm(z1)/stssigma)) +
sum(log(dnorm(z2)/strsigma))
}

#profile function
profile.ss <- function(stsObs, strObs, Uhats) {
plkuhat <- rep(0, length(Uhats))

for (i in 1:length(Uhats)) {
temp <- optim(c(log(stsSigma), log(strSigma), strMu), llik.ss,
stsObs=stsObs, strObs=strObs, Uhat=Uhats[i] )
plkuhat[i] <- temp$value } plkuhat } #MLE for Rhat Rhat #generate fixed values for Rhat Rhats <- seq(.75, .999, length.out=25) Uhats <- 1-Rhats #likelihood profile (=computationally intensive) #note: in case of warning messages, or irregular profile plot, #decrease the width of the interval for possible values of Xmu #in the llik.ss function #that is, change the bounds to, for instance, stsMu+/-10*stsMuSE ll <- profile.ss(Uhats=Uhats, stsObs=stress, strObs=strength) #plot profile likelihood for Rhat plot(Rhats,ll, type='l', xlab="Reliability", ylab="log-likelihood") #include lower limit for log-likelihood values loglikw <- llik.ss(modpar=c(log(stsSigma), log(strSigma), strMu), stsObs=stress, strObs=strength, Uhat=Uhat) alpha <- .05 #95% confidence interval limit <- loglikw - .5*qchisq(1-alpha, df=1) abline(h=limit,col="blue",lty=2) #include MLE for Rhat abline(v=Rhat, col="red", lty=2) #compute limits of confidence interval by hand (=computationally intensive) limitlkh <- function(stsObs, strObs, Uhats) { profile.ss(stsObs, strObs, Uhats) - limit } lowerci <- uniroot(limitlkh, c(Uhat, 1-.75), stsObs=stress, strObs=strength) upperci <- uniroot(limitlkh, c(1-1, Uhat), stsObs=stress, strObs=strength) round(c(1-lowerci$root, 1-upperci$root), 8) #include confidence bounds in plot abline(v=1-lowerci$root, col="darkgreen", lty=2)
abline(v=1-upperci$root, col="darkgreen", lty=2) #compare this confidence interval with that computed by the Reiser-Guttman method estSSR(strength, stress, type="RG", alpha=.05)$CI
#the Reiser-Gutmann method yields on average good results
#that is, the coverage probability of these intervals is often very close to the nominal level

#lower bound of one sided interval
alpha <- .05 #95% confidence interval
limit <- loglikw - .5*qchisq(1-2*alpha, df=1)

lowerci <- uniroot(limitlkh, c(Uhat, 1-.75), stsObs=stress, strObs=strength)
1-lowerci$root #compare confidence interval with that computed by the Guo-Krishnamoorthy method estSSR(strength, stress, type="GK", alpha=.05, twoside=FALSE)$CI
#this Guo-Krishnamoorthy method yields on average excellent results
#that is, the coverage probability of these intervals is often very close to the nominal level

##stress and strength both follow a Weibull distribution

##data

#generate artificial data

#stress data
stress <- exp(7 + .3*rsev(40))

#strength data
strength <- exp(8 + .6*rsev(30))

#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")

#obtain MLE for mu and sigma
#Lmu
stsMu <- coef(sts)
names(stsMu) <- "stsMu"
#Lsigma
stsSigma <- sts$scale names(stsSigma) <- "stsSigma" #Smu strMu <- coef(str) names(strMu) <- "strMu" #Ssigma strSigma <- str$scale
names(strSigma) <- "strSigma"

##compute likelihood based interval

#specify integrand for 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
}

#estimate of unreliability using numerical integration
uplimit <- exp(qsev(.99999999)*sts$scale+coef(sts)) #upper limit of the integral 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)

#function for computing abs(Uhat-f(stsMu,stsSigma,strMu,stsSigma))
stsMuMin <- function (stsMuPar, stsSigma, strMu, strSigma, Uhat){
uplimit <- exp(qsev(.99999999)*stsSigma+stsMuPar)
abs(Uhat-integrate(integrand, 0, uplimit, mu1=stsMuPar, mu2=strMu,
s1=stsSigma, s2=strSigma)$value)} #log-likelihood function stsMuSE <- sqrt(diag(vcov(sts))) #extract standard error for stsMu llik.ss <- function (modpar, stsObs, strObs, Uhat) { stssigma <- exp(modpar) strsigma <- exp(modpar) strmu <- modpar stsmu <- optim(par=stsMu, stsMuMin, method="Brent", lower=stsMu-20*stsMuSE, upper=stsMu+20*stsMuSE, Uhat=Uhat, stsSigma=stssigma, strMu=strmu, strSigma=strsigma)$par

z1 <- (log(stsObs)-stsmu)/stssigma
z2 <- (log(strObs)-strmu)/strsigma
#assuming that the stress and strength random variables are independent,
#the log-likelihood function is given by:
sum(log(dsev(z1)/(stsObs*stssigma))) +
sum(log(dsev(z2)/(strObs*strsigma)))
}

#profile function
profile.ss <- function(stsObs, strObs, Uhats) {
plkuhat <- rep(0, length(Uhats))

for (i in 1:length(Uhats)) {
temp <- optim(c(log(stsSigma), log(strSigma), strMu), llik.ss,
plkuhat[i] <- temp$value } plkuhat } #MLE for Rhat Rhat #generate fixed values for Rhat Rhats <- seq(.60, .99, length.out=25) Uhats <- 1-Rhats #likelihood profile (=computationally intensive) #note: in case of warning messages, or irregular profile plot, #decrease the width of the interval for possible values of Xmu #in the llik.ss function #that is, change the bounds to, for instance, stsMu+/-10*stsMuSE ll <- profile.ss(Uhats=Uhats, stsObs=stress, strObs=strength) #plot profile likelihood for Rhat plot(Rhats,ll, type='l', xlab="Reliability", ylab="log-likelihood") #include lower limit for log-likelihood values loglikw <- llik.ss(modpar=c(log(stsSigma), log(strSigma), strMu), stsObs=stress, strObs=strength, Uhat=Uhat) alpha <- .05 #95% confidence interval limit <- loglikw - .5*qchisq(1-alpha, df=1) abline(h=limit,col="blue",lty=2) #include MLE for Rhat abline(v=Rhat, col="red", lty=2) #compute limits of confidence interval by hand (=computationally intensive) limitlkh <- function(stsObs, strObs, Uhats) { profile.ss(stsObs, strObs, Uhats) - limit } lowerci <- uniroot(limitlkh, c(Uhat, 1-.60), stsObs=stress, strObs=strength) upperci <- uniroot(limitlkh, c(1-.99, Uhat), stsObs=stress, strObs=strength) round(c(1-lowerci$root, 1-upperci$root), 8) #include confidence bounds in plot abline(v=1-lowerci$root, col="darkgreen", lty=2)
abline(v=1-upperci$root, col="darkgreen", lty=2) #lower bound of one sided interval alpha <- .05 #95% confidence interval limit <- loglikw - .5*qchisq(1-2*alpha, df=1) lowerci <- uniroot(limitlkh, c(Uhat, 1-.65), stsObs=stress, strObs=strength) 1-lowerci$root