The SSReliabilityClaytonMWD package provides tools for
modeling and estimating stress–strength reliability under dependent
structures using the Modified Weibull Distribution (MWD) and a
copula-based dependence model.
This vignette presents an end-to-end workflow for estimating
reliability in dependent stress–strength models using the
SSReliabilityClaytonMWD package.
Stress–strength reliability (SSR) is a fundamental concept in reliability engineering and survival analysis, defined as: \[ R = P(X>Y) \] where \(X\) denotes the strength variable and \(Y\) denotes the stress variable.
The vignette focuses on:
Install the latest released version from CRAN:
Install the latest development version from GitHub:
Let the strength \(X\) and stress
\(Y\) variables follow the Modified
Weibull Distribution (MWD): \(X \sim
MWD(a_1,b_1,\lambda_1)\) and \(Y \sim
MWD(a_2,b_2,\lambda_2)\). Their dependence is modeled using the
two-dimensional Clayton copula, see Clayton.Copula,
\[C_{\theta}(u,v) = (u^{-\theta} +
v^{-\theta} - 1)^{-1/\theta}.\] Under this model, the
stress–strength reliability \(R\) is
given by: \[
R = P(X>Y) = \int_{0}^{1} t^{-(\theta +1)} \big( t^{-\theta} +
G_Y(F_{X}^{-1}(t))^{-\theta} -1 \big)^ {-\left(\frac{1}{\theta}+1
\right) } \tag{1},
\] where \(F_X(x) \equiv
F_X(x;a_1,b_1,\lambda_1) = 1- \exp(-a_1 x^{b_1} e^{\lambda_1
x})\) and \(G_Y(y) \equiv
G_Y(y;a_2,b_2,\lambda_2) = 1- \exp(-a_2 y^{b_2} e^{\lambda_2
y})\). It is calculated as in Reliability_Clayton_MWD
by numerical integral stats::integrate.
Further details can be found in the companion vignette Dependent Stress–Strength Reliability Model and in Kızılaslan (2026).
Based on these assumptions, suppose \(n\) systems are observed, yielding dependent samples \((x_i,y_i)\; i=1,\dots, n\).
The joint log-likelihood function can be decomposed as: \[ \ell(\boldsymbol{\Omega_1},\boldsymbol{\Omega_2},\theta; \mathbf{x},\mathbf{y}) = \ell_1(\boldsymbol{\Omega_1}; \mathbf{x}) + \ell_2(\boldsymbol{\Omega_2}; \mathbf{y}) + \ell_3(\boldsymbol{\Omega_1},\boldsymbol{\Omega_2},\theta; \mathbf{x},\mathbf{y}), \] where \[ \ell_1(\boldsymbol{\Omega_1}; \mathbf{x}) = n \log(a_1) + \sum_{i=1}^{n} \log(b_1+\lambda_1 x_i) + (b_1-1) \sum_{i=1}^{n} \log x_i + \lambda_1 \sum_{i=1}^{n} x_i - a_1 \sum_{i=1}^{n}x_i^{b_1} e^{\lambda_1 x_i}, \]
\[ \ell_2(\boldsymbol{\Omega_2}; \mathbf{y}) = n \log(a_2) + \sum_{i=1}^{n} \log(b_2+\lambda_2 y_i) + (b_2-1) \sum_{i=1}^{n} \log y_i + \lambda_2 \sum_{i=1}^{n} y_i - a_2 \sum_{i=1}^{n} y_i^{b_2} e^{\lambda_2 y_i}, \]
\[ \ell_3(\boldsymbol{\Omega_1},\boldsymbol{\Omega_2},\theta; \mathbf{x},\mathbf{y}) = \sum_{i=1}^{n} \log \big( c_{\theta} \big( F_{X}(x_i; \boldsymbol \Omega_1), G_Y(y_i; \boldsymbol \Omega_2) \big) \big), \] and \(c_{\theta}\) is the joint PDF of Clayton copula.
We employ the two-step maximum likelihood approach, also known as the
inference functions for margins (IFM) method. In the first step, the
marginal parameters are estimated separately by maximizing their
respective log-likelihood functions. In the second step, the dependence
parameter \(\theta\) is estimated by
maximizing \(\ell_3\) given the
estimated marginals. This separation simplifies the estimation procedure
and improves numerical stability. In both steps, optimization is carried
out using stats::optim.
Finally, the the MLE of \(R\), \(\hat{R}\), is obtained by substituting the parameter estimates into Equation (1).
The least squares (LS) estimator is obtained by minimizing the sum of squared differences between the theoretical cumulative distribution function (CDF) and the empirical CDF evaluated at the ordered sample points.
We apply the same two-step procedure as in the MLE case for both the LSE and WLSE methods.
Let \(x_{(1)}<x_{(2)}< \dots < x_{(n)}\) and \(y_{(1)}<y_{(2)}< \dots < y_{(n)}\) denote the ordered samples. The LS estimates of the marginal parameters \(\boldsymbol{\Omega_1}\) and \(\boldsymbol{\Omega_2}\) are obtained by minimizing:
\[ Q_1(\boldsymbol{\Omega_1}) = \sum_{i=1}^n \bigg( F(x_{(i)};a_1,b_1,\lambda_1) - \frac{i-0.3}{n+0.4} \bigg)^2 \tag{2}, \] and \[ Q_2(\boldsymbol{\Omega_2}) = \sum_{i=1}^n \bigg( G(y_{(i)};a_2,b_2,\lambda_2) - \frac{i-0.3}{n+0.4} \bigg)^2 \tag{3}. \] In the second step, the LSE of the Clayton copula parameter \(\theta\) is obtained by minimizing: \[ Q_3(\theta) = \sum_{i=1}^n \bigg( C_{\theta}(u_i,v_i;\theta) - \hat{H}_n(x_i,y_i) \bigg)^2 \tag{4}, \] where \(u_i= F(x_i; \hat{a}_{1,LSE}, \hat{b}_{1,LSE}, \hat{\lambda}_{1,LSE})\), \(v_i= G(y_i; \hat{a}_{2,LSE}, \hat{b}_{2,LSE}, \hat{\lambda}_{2,LSE})\) and \(\hat{H}_n(x_i,y_i) = \frac{1}{n}\sum_{j=1}^n \mathbf{I}(X_j \leq x_i, Y_j \leq y_i)\).
Finally, the LSE of \(R\), \(\hat{R}_{LSE}\), is obtained by substituting the parameter estimates into Equation (1).
In the weighted least squares approach, weights are introduced to improve estimation accuracy. The weights are defined as: \[ w_i = \frac{(n+1)^2(n+2)}{i(n-i+1)}, \; i=1,\dots,n. \] These weights are incorporated into Equations (2)–(4) as follows: \[ Q_1^W(\boldsymbol{\Omega_1}) = \sum_{i=1}^n w_i \bigg( F(x_{(i)};a_1,b_1,\lambda_1) - \frac{i-0.3}{n+0.4} \bigg)^2, \]
\[ Q_2^W(\boldsymbol{\Omega_2}) = \sum_{i=1}^n w_i \bigg( G(y_{(i)};a_2,b_2,\lambda_2) - \frac{i-0.3}{n+0.4} \bigg)^2, \]
\[ Q_3^W(\theta) = \sum_{i=1}^n w_i \bigg( C_{\theta}(u_i,v_i;\theta) - \hat{H}_n(x_i,y_i) \bigg)^2. \] In the first step, the weighted objective functions \(Q_1^W\) and \(Q_2^W\) are minimized to obtain the marginal parameter estimates. In the second step, the dependence parameter \(\theta\) is estimated by minimizing \(Q_3^W\).
The WLSE of \(R\), \(\hat{R}_{WLSE}\), is obtained by substituting the corresponding estimates into Equation (1).
The maximum product of spacings (MPS) method is based on maximizing the geometric mean of spacings between consecutive values of the fitted distribution functions.
Let \(x_{(1)}<x_{(2)}< \dots < x_{(n)}\) and \(y_{(1)}<y_{(2)}< \dots < y_{(n)}\) denote the ordered samples for strength and stress, respectively.
The spacings for the marginal distributions are defined as: \[ D_{1i}(a_1,b_1,\lambda_1) = F(x_{(i)};a_1,b_1,\lambda_1) - F(x_{(i-1)};a_1,b_1,\lambda_1), \] and \[ D_{2i}(a_2,b_2,\lambda_2) = G(y_{(i)};a_2,b_2,\lambda_2) - G(y_{(i-1)};a_2,b_2,\lambda_2), \] where \(i=1,\dots, n+1\), \(F(x_{(0)};a_1,b_1,\lambda_1) = G(y_{(0)};a_2,b_2,\lambda_2) =0\) and \(F(x_{(n+1)};a_1,b_1,\lambda_1) = G(y_{(n+1)};a_2,b_2,\lambda_2) =1\).
In the first step, the MPS estimates of the marginal parameters \(\boldsymbol{\Omega_1}\) and \(\boldsymbol{\Omega_2}\) are obtained by maximizing: \[ M_1(a_1,b_1,\lambda_1) = \sum_{i=1}^{n+1} \frac{1}{n+1} \log \big(D_{1i} (a_1,b_1,\lambda_1) \big), \] and \[ M_2(a_2,b_2,\lambda_2) = \sum_{i=1}^{n+1} \frac{1}{n+1} \log \big(D_{2i} (a_2,b_2,\lambda_2) \big). \] In the second step, the MPS estimate of the Clayton copula parameter \(\theta\) is obtained by maximizing: \[ M_3(\theta) = \sum_{i=1}^{n+1} \frac{1}{n+1} \log \big(D_{3i} (\theta ) \big) \] where \(D_{3i}(\theta) = C_{\theta}\big( F(x_{(i)}; \widehat{\boldsymbol \Omega}_{1,MPS}), G(y_{(i)};\widehat{\boldsymbol \Omega}_{2,MPS} ); \theta \big) - C_{\theta}\big( F(x_{(i-1)}; \widehat{\boldsymbol \Omega}_{1,MPS}), G(y_{(i-1)}; \widehat{\boldsymbol \Omega}_{2,MPS} ); \theta \big)\).
Finally, the MPS estimate of \(R\), \(\hat{R}_{MPS}\), is obtained by substituting the parameter estimates into Equation (1).
The parametric bootstrap percentile method is employed to be able to construct bootstrap confidence intervals (BCIs) for unknown parameters and \(R\) as well based on all used methods here. The following Algorithm is employed for obtaining BCIs of the interested parameters based on MLE, LSE, WLSE and MPS methods.
We construct bootstrap confidence intervals for \(R\) using the following algorithm:
Step 1. Compute the estimates \(\hat{a}_1, \hat{b}_1, \hat{\lambda}_1, \hat{a}_2, \hat{b}_2\) and \(\hat{\lambda}_2\) using one of the estimation methods: ML, LS, WLS, or MPS, based on the observed dependent sample.
Step 2. Generate a dependent bootstrap sample using the estimated parameters.
Step 3. Compute bootstrap estimates \(\hat{a}_1^{*}, \hat{b}_1^{*}, \hat{\lambda}_1^{*}, \hat{a}_2^{*}, \hat{b}_2^{*}\), and \(\hat{\lambda}_2^{*}\) based on the bootstrap sample.
Step 4. Compute the bootstrap estimate of reliability \(R\), denoted by \(\hat{R}^{*}\).
Step 5. Repeat Steps 2–4 independently \(B\) times to obtain \(\{\hat{R}_1^{*},\ldots,\hat{R}_B^{*}\}\).
Step 6. The two-sided \(100(1-\gamma)\%\) bootstrap confidence interval for \(R\) is: \(\left( \hat{R}^{*[\gamma/2]}, \hat{R}^{*[1-\gamma/2]} \right)\), where \(\hat{R}^{*[p]}\) denotes the \(p\)-th empirical quantile of the bootstrap replications. BCIs for the parameters \(a_1, b_1, \lambda_1, a_2, b_2,\) and \(\lambda_2\) are obtained in the same manner.
To illustrate the performance and applicability of the proposed methods, we consider both simulated and real data examples.
A simulation study is conducted to evaluate the performance of the proposed estimators under the assumed model structure.
# Load the package
library(SSReliabilityClaytonMWD)
# set seed
set.seed(123)
n <- 50
a1 <- 0.75; b1 <- 1.5; lambda1 <- 0.6
a2 <- 1.2; b2 <- 0.5; lambda2 <- 0.9
theta <- 3
# simulate data
dat <- rMweibull_Clayton(n, a1, b1, lambda1, a2, b2, lambda2, theta)# true stress-strength reliability value
R_true <- Reliability_Clayton_MWD(a1, b1, lambda1, a2, b2, lambda2, theta)
R_true$value
#> [1] 0.9428307We fit a dependent stress–strength reliability model assuming that both strength and stress follow the MWD, and their dependence is modeled using a Clayton copula.
fit <- fit.SSR.ClaytonMWD(
data = dat,
ACI = TRUE,
bootstrap = TRUE,
B = 10,
seed = 2026,
one.step = TRUE,
alpha = 0.05
)
#> Warning: executing %dopar% sequentially: no parallel backend registeredThe function returns point estimates of all model parameters together with \(95\%\) asymptotic confidence intervals and bootstrap confidence intervals for both the parameters and the reliability \(R\).
print(fit)
#>
#> Fitting Results of Stress-Strength Reliability Model with MWD Marginals via Clayton Copula
#>
#>
#> Table: Parameter Estimates
#>
#> | | a1| b1| lambda1| a2| b2| lambda2| theta| R|
#> |:----|-------:|-------:|-------:|-------:|-------:|-------:|-------:|-------:|
#> |mle | 0.91531| 1.77111| 0.33610| 1.15016| 0.51926| 1.07039| 2.98107| 0.96239|
#> |lse | 1.26967| 1.91057| 0.00001| 0.94711| 0.44689| 1.32354| 9.72676| 0.99555|
#> |wlse | 1.10497| 1.83511| 0.14438| 0.96263| 0.45805| 1.32415| 3.01138| 0.96637|
#> |mps | 0.94898| 1.68266| 0.27554| 1.14662| 0.48797| 0.97071| 3.28740| 0.95748|
#>
#>
#> Table: 95% Asymptotic Confidence Intervals (MLE)
#>
#> | | a1| b1| lambda1| a2| b2| lambda2| theta| R|
#> |:------|-------:|-------:|-------:|-------:|-------:|-------:|-------:|-------:|
#> |lower | 0.00000| 0.67220| 0.00000| 0.39511| 0.28553| 0.32254| 1.56639| 0.91894|
#> |upper | 2.19113| 2.87001| 1.59347| 1.90520| 0.75298| 1.81823| 4.39576| 1.00000|
#> |length | 2.19113| 2.19781| 1.59347| 1.51009| 0.46745| 1.49569| 2.82937| 0.08106|
#>
#>
#> Table: 95% Bootstrap CIs (MLE)
#>
#> | | a1| b1| lambda1| a2| b2| lambda2| theta| R|
#> |:------|-------:|-------:|-------:|-------:|-------:|-------:|-------:|-------:|
#> |lower | 0.26216| 0.77866| 0.00001| 0.68765| 0.33817| 0.78146| 1.90427| 0.95857|
#> |upper | 1.50917| 2.24596| 1.55287| 2.05227| 0.67876| 1.85251| 4.56875| 0.98110|
#> |length | 1.24701| 1.46731| 1.55286| 1.36462| 0.34059| 1.07104| 2.66448| 0.02253|
#>
#>
#> Table: 95% Bootstrap CIs (LSE)
#>
#> | | a1| b1| lambda1| a2| b2| lambda2| theta| R|
#> |:------|-------:|-------:|-------:|-------:|-------:|-------:|--------:|-------:|
#> |lower | 0.25455| 0.62619| 0.00001| 0.40943| 0.19347| 0.18018| 5.72943| 0.98075|
#> |upper | 1.63375| 2.22374| 1.59905| 2.96950| 0.89047| 2.54395| 21.17203| 0.99974|
#> |length | 1.37920| 1.59755| 1.59904| 2.56007| 0.69699| 2.36377| 15.44261| 0.01898|
#>
#>
#> Table: 95% Bootstrap CIs (WLSE)
#>
#> | | a1| b1| lambda1| a2| b2| lambda2| theta| R|
#> |:------|-------:|-------:|-------:|-------:|-------:|-------:|--------:|-------:|
#> |lower | 0.26556| 0.69049| 0.00001| 0.49786| 0.26531| 0.63087| 1.61708| 0.94712|
#> |upper | 1.59581| 2.30977| 1.50784| 2.09689| 0.72586| 2.22361| 11.47759| 0.99625|
#> |length | 1.33025| 1.61928| 1.50783| 1.59904| 0.46055| 1.59274| 9.86052| 0.04913|
#>
#>
#> Table: 95% Bootstrap CIs (MPS)
#>
#> | | a1| b1| lambda1| a2| b2| lambda2| theta| R|
#> |:------|-------:|-------:|-------:|-------:|-------:|-------:|-------:|-------:|
#> |lower | 0.30178| 0.67258| 0.00001| 0.69724| 0.29324| 0.64941| 1.58618| 0.92730|
#> |upper | 1.41056| 2.02045| 1.36527| 1.93937| 0.61315| 1.61662| 5.94632| 0.97949|
#> |length | 1.10878| 1.34787| 1.36526| 1.24213| 0.31991| 0.96721| 4.36014| 0.05219|
#>
#> Kendall's tau estimate for theta: 3.30303The proposed methodology is applied to real-world data to demonstrate its practical relevance. Specifically, monthly occupancy data from two of Istanbul’s largest dams are used. The data are obtained from the open data platform of the Istanbul Metropolitan Municipality: https://data.ibb.gov.tr/en/.
Within the stress–strength framework, the Terkos dam is treated as the strength variable \(X\), while the Omerli dam is considered as the stress variable \(Y\). Under this formulation, the stress–strength reliability \(R=P(X>Y)\) represents the probability that the occupancy level of the Terkos dam (European side) exceeds that of the Omerli dam (Anatolian side). This probability provides a measure of the relative water availability between the two regions. Further details can be found in Kızılaslan (2026).
# Load example data from the package
data(TerkosDam)
data(OmerliDam)
real_data <- list(X = TerkosDam, Y = OmerliDam)We fit the dependent stress–strength reliability model assuming that both variables follow the MWD, with dependence modeled via a Clayton copula.
fit_ssr <- fit.SSR.ClaytonMWD(
data = real_data,
ACI = TRUE,
bootstrap = TRUE,
B = 10,
seed = 2026,
one.step = TRUE,
alpha = 0.05
)The fitted results for the model parameters and the reliability \(R\) are presented below.
print(fit_ssr)
#>
#> Fitting Results of Stress-Strength Reliability Model with MWD Marginals via Clayton Copula
#>
#>
#> Table: Parameter Estimates
#>
#> | | a1| b1| lambda1| a2| b2| lambda2| theta| R|
#> |:----|-------:|-------:|-------:|-------:|-------:|-------:|-------:|-------:|
#> |mle | 0.98542| 2.66095| 2.11322| 0.02428| 0.45908| 6.33059| 0.50551| 0.50428|
#> |lse | 0.07196| 1.72091| 5.71667| 0.01166| 0.00001| 7.06623| 1.09884| 0.49349|
#> |wlse | 0.04261| 1.29079| 6.20629| 0.01154| 0.00001| 7.10656| 1.23631| 0.49824|
#> |mps | 0.68997| 2.31111| 2.39000| 0.02043| 0.26258| 6.43699| 0.92547| 0.51110|
#>
#>
#> Table: 95% Asymptotic Confidence Intervals (MLE)
#>
#> | | a1| b1| lambda1| a2| b2| lambda2| theta| R|
#> |:------|-------:|-------:|-------:|-------:|-------:|-------:|-------:|-------:|
#> |lower | 0.00000| 0.43992| 0.00000| 0.00000| 0.00000| 4.21722| 0.14734| 0.40310|
#> |upper | 4.33682| 4.88197| 5.79639| 0.06723| 1.39756| 8.44395| 0.86369| 0.60545|
#> |length | 4.33682| 4.44205| 5.79639| 0.06723| 1.39756| 4.22673| 0.71636| 0.20235|
#>
#>
#> Table: 95% Bootstrap CIs (MLE)
#>
#> | | a1| b1| lambda1| a2| b2| lambda2| theta| R|
#> |:------|-------:|-------:|-------:|------:|-------:|-------:|-------:|-------:|
#> |lower | 0.08013| 1.15190| 0.00001| 0.0127| 0.14660| 2.96074| 0.24032| 0.44208|
#> |upper | 7.84406| 4.17488| 4.99832| 0.6515| 2.74715| 7.36907| 0.76629| 0.57280|
#> |length | 7.76393| 3.02298| 4.99831| 0.6388| 2.60056| 4.40833| 0.52597| 0.13072|
#>
#>
#> Table: 95% Bootstrap CIs (LSE)
#>
#> | | a1| b1| lambda1| a2| b2| lambda2| theta| R|
#> |:------|--------:|-------:|-------:|-------:|-------:|-------:|-------:|-------:|
#> |lower | 0.00432| 0.00001| 0.08203| 0.00479| 0.00001| 3.38092| 0.77300| 0.40633|
#> |upper | 14.07044| 5.07941| 8.83643| 0.46871| 2.90343| 8.70306| 1.22167| 0.54812|
#> |length | 14.06612| 5.07940| 8.75439| 0.46392| 2.90342| 5.32214| 0.44867| 0.14178|
#>
#>
#> Table: 95% Bootstrap CIs (WLSE)
#>
#> | | a1| b1| lambda1| a2| b2| lambda2| theta| R|
#> |:------|-------:|-------:|-------:|-------:|-------:|-------:|-------:|-------:|
#> |lower | 0.00600| 0.08335| 2.07582| 0.00735| 0.00001| 4.57432| 1.02051| 0.40990|
#> |upper | 1.93694| 3.81682| 8.40554| 0.14362| 2.18957| 7.75921| 2.79450| 0.55191|
#> |length | 1.93094| 3.73347| 6.32972| 0.13627| 2.18956| 3.18489| 1.77400| 0.14201|
#>
#>
#> Table: 95% Bootstrap CIs (MPS)
#>
#> | | a1| b1| lambda1| a2| b2| lambda2| theta| R|
#> |:------|-------:|-------:|-------:|-------:|-------:|-------:|-------:|-------:|
#> |lower | 0.03681| 0.81306| 0.20976| 0.00853| 0.00001| 3.89184| 0.39588| 0.43467|
#> |upper | 4.90435| 3.55720| 5.68968| 0.23335| 1.80479| 7.53494| 5.25280| 0.62724|
#> |length | 4.86754| 2.74414| 5.47992| 0.22482| 1.80478| 3.64310| 4.85692| 0.19257|
#>
#> Kendall's tau estimate for theta: 1.25912Dependent Stress–Strength Reliability
Model
vignette("ssr-theory", package = "SSReliabilityClaytonMWD")
Modified Weibull Distribution (MWD)
vignette("mwd-distribution", package = "SSReliabilityClaytonMWD")
Fitting the Modified Weibull Distribution (MWD)
vignette("fitting-mwd", package = "SSReliabilityClaytonMWD")
Kızılaslan, Fatih. (2026). Reliability estimation in dependent stress–strength model with Clayton copula and modified Weibull margins.arXiv:2604.12130