1. Introduction
The stress-strength model is extensively utilized in the context of reliability analysis. This model characterizes the life of a component with a random strength denoted by X and subject to a random stress represented by Y. The component fails if the applied stress exceeds the strength (Y > X); conversely, it continues to function if the strength surpasses the stress (X > Y). Therefore, 
$R = P(X \gt Y)$ serves as a metric for the reliability of a component. Stress-strength models are applied in various research domains, particularly in engineering. Examples include modeling the degradation of concrete pressure vessels, the decay of rocket motors, the static fatigue of ceramic components, and the fatigue failure of aircraft structures [Reference Nadarajah35]. For additional applications spanning engineering, quality control, medicine, and psychology, please refer to [Reference Jeevanand22], [Reference Kotz, Lumelskii and Pensky29], and [Reference Nadarajah and Kotz36].
In the context of single-component stress-strength models, various lifetime distributions for the stress and strength random variables have been assumed to estimate reliability. Significant contributions to this topic include studies on the exponential distribution (ba14), Weibull distribution (kundu06), Burr XII distribution (lio12), Lindley distribution (al13), Kumaraswamy distribution (nada14), standard two-sided power distribution (kaya2019), and the truncated proportional hazard rate distribution (bai19). For more comprehensive information, please refer to the extensive coverage of this subject provided by [Reference Kotz, Lumelskii and Pensky29].
With the development of society and advancements in science and technology, the complexity of products and the diversity of services have led to an increasing number of complex systems, the vast majority of which are multicomponent systems. How to conduct reliability research on multicomponent systems has become a topic of widespread concern. One of the key issues is the stress-strength model for multicomponent systems. Hanagal [Reference Hanagal16] initially estimated the system reliability under the premise that the strengths of the k components 
$(X_1, X_2,\ldots, X_k)$ are influenced by a statistically independent common stress 
$X_{k+1}$. Moreover, Kotz et al. [Reference Kotz, Lumelskii and Pensky29] provided both practical and theoretical insights into the theory and application of stress-strength models within economic and industrial systems. In recent years, many researchers have considered the multicomponent stress-strength model reliability under complete samples, for example, Rao and Kantam [Reference Rao and Kantam41] addressed multicomponent stress-strength model reliability using the log-logistic distribution. Kızılaslan and Nadar [Reference Kizilaslan and Nadar26] focused on Weibull distribution for multicomponent stress-strength model reliability estimation. Nadar and Kızılaslan [Reference Nadar and Kızılaslan33] estimated the reliability of a multicomponent stress-intensity model with k independent and equally distributed strength components, each consisting of a pair of dependent elements following a Marshall–Olkin binary Weibull distribution, each exposed to a common random stress following a Weibull distribution. Rao et al. [Reference Rao, Muhammad and H42] investigated multicomponent stress-strength reliability using the exponentiated Weibull distribution. Hassan and Alohali [Reference Hassan and Alohali17] estimated multicomponent stress-strength reliability with a generalized linear failure rate distribution. Akgül [Reference Akgül1] explored the estimation of multicomponent stress-strength reliability through classical and Bayesian approaches, assuming both stress and strength followed the Topp–Leone distribution. Kayal et al. [Reference Kayal, Tripathi, Dey and Wu24] investigated the reliability of multicomponent stress-strength systems using the Chen distribution. Rasekhi et al. [Reference Rasekhi, Saber and Yousof43] employed Bayesian and classical inference techniques to estimate reliability in multicomponent stress-strength system, assuming the stress and strength followed a generalized logistic distribution. Bai et al. [Reference Bai, Li, Balakrishnan and He4] discussed inferential procedures for stress-strength reliability of multistate system using generalized survival signature under the assumption of Gumbel copula dependence between strength variables. Kotb and Raqab [Reference Kotb and Raqab28] investigated reliability estimation in multicomponent stress-strength system using a modified Weibull distribution. Jana and Bera [Reference Jana and Bera21] studied interval reliability estimation for a multicomponent stress-strength system with inverse Weibull-distributed stress and strength components. Bai et al. [Reference Bai, Zhang, He and Balakrishnan6] discussed inferential procedures for stress-strength reliability of the multistate system using the proposed improved generalized survival signature, under the stress and strength are independent, while stress is dependent between Gumbl and Clayton copula.
It is essential to highlight that most of the models discussed in the aforementioned literature are built on the assumption that stress and strength are independent variables. Nevertheless, this assumption may not always hold in real-world situations. For instance, mechanical structures may demonstrate different fatigue lives under various stress levels, and a structure’s fatigue life essentially reflects its capacity to resist structural failures. In addition, Domma and Giordano [Reference Domma and Giordano11] elaborated the possibility of interdependence between stress and strength from six aspects, including engineering, operations research, quality control, economics, education, and insurance; furthermore, they aimed to fill a gap by evaluating reliability of stress-strength model taking into account the association between stress and strength via Farlie–Gumbel–Morgenstern (FGM) copula. Consequently, the stress-strength model incorporating stress-dependent strength has emerged as a focal point of interest among researchers. Vaidyanathan and VA [Reference Vaidyanathan and VA45] derived the expression of R when stress and strength follow bivariate Lindley distribution. Papadrakakis et al. [Reference Papadrakakis, Papadopoulos and Stefanou39] investigated stress-strength reliability estimation using the Extended FGM and Ali–Mikhail–Haq (AMH) copula with exponential distribution. The stress-strength model is investigated by [Reference Patil and Naik-Nimbalkar40] with stress and strength margin belonging to the Pareto family and the dependency is represented using four different types of copulas. Pak and Gupta [Reference Pak and Gupta38] discussed on estimation of the parameter R where more realistically stress and strength are dependent random variables distributed as bivariate Rayleigh model. Zhu [Reference Zhu46] evaluated multicomponent system reliability, in which, the stress and strength are assumed to have dependent Kumaraswamy variable and unit Gompertz variable based on Clayton copula. James et al. [Reference James, Chandra and Sebastian20] obtained the estimates of reliability and dependence parameters under the assumption stress and strength are linked by FGM copula with Rayleigh marginals as the underlying distribution.
To the best of our knowledge, until now a similar task has never been attempted for evaluation of R where strengths are dependent, which are dependent of the common stress (see Table 1). However, this problem is common in daily life, for example, the various components of a system are often interdependent due to shared environmental factors, leading to correlated strengths among these components. In addition to the interdependence of strengths within the system, there also exists a relationship between the system’s strength and the applied stress it experiences (that is, stronger systems tend to withstand higher levels of stress). Generally speaking, the dependence between the strength of the components of the system is not the same as the dependence between the strength of the system and the stress to which it is subjected. In addressing this theoretical challenge, the main interest of this attempt is to formulate dependence stress-strength reliability model with hierarchical Archimedean copula (HAC) approach. HACs provide a natural way to model the dependence structure in multivariate distributions by incorporating multiple levels of dependence through a hierarchical framework, which introduced by [Reference Joe23] as an extension of Archimedean copulas (ACs). Our works will improve the corresponding results of multicomponent stress-strength model from the independent assumption on component’s strength to the statistically dependent setting under complete sample.
Table 1. The relation of stress and strength in stress-strength model.
The predominant focus of research pertaining to the estimation of reliability measures is predominantly situated within a parametric framework, encompassing both classical and Bayesian methodologies. It is widely acknowledged that statistical models rooted in traditional parametric distributions may lack the necessary adaptability to offer a dependable portrayal of survival data. As a result, nonparametric estimators have gained substantial traction within the field. Semi-parametric estimators, seen as an intermediary approach, endeavor to strike a harmonious balance between conventional parametric techniques and their nonparametric counterparts. In a concerted effort to address this gap, the purpose of this paper is to study the estimation of of multicomponent stress-strength reliability using four semi-parametric approach based on a complete sample when there is hierarchical dependence among stresses, as well as between stresses and strengths. The subsequent sections of this paper are structured as follows. Section 2 provides an exposition on HAC theory, model delineation, and fundamental assumptions. In Section 3, we introduce the semi-parametric estimation approach for the dependence parameter and the estimation method for R within the multicomponent stress-strength model and establish the properties of the estimator. Illustrative simulations and presentation of real data analysis in Section 4. Concluding remarks are outlined in Section 5.
2. Preliminaries
Before proceeding to the main results, let us first recall concepts of copula and HAC, which will be used in the sequel.
2.1. Archimedean copula
The copula function has demonstrated its versatility in characterizing the relationship between variables, irrespective of their individual marginal behaviors. For readers new to the concept of copula and its applications, the foundational sources are the monographs by [Reference Joe23] and [Reference Nelsen37]. Additionally, Durante and Sempi [Reference Durante and Sempi12] furnished a comprehensive compilation of references pertaining to copulas. Next, we provide the definition of a copula.
Definition 2.1. A copula is a function 
$C:I^n\to I$, where 
$I=[0,1]$, with the following properties
(i)
$C(v_1,v_2,\ldots,v_n)$ is increasing in
$v_i,\;i=1,2,\ldots,n,$(ii)
$C(1,\ldots,1, v_i,1,\ldots 1)=v_i,$ for all
$i=1,2,\ldots,n,$(iii)
$C(0,\ldots,0)=0$ and
$ C(1,\ldots,1)=1,$ and(iv) for any
$x_1,x_2,\ldots,x_n,\;y_1,y_2,\ldots,y_n,$ if
$x_j\leq y_j,\;j=1,2,\ldots,n,$ then
\begin{equation*}\sum_{i_1=1}^{2}\sum_{i_2=1}^{2}\cdots\sum_{i_n=1}^{2}(-1)^{i_1+i_2+\cdots+i_n}C(v_{1i_1},v_{2i_2},\ldots, v_{ni_n})\geq 0,\end{equation*}where
$v_{j1}=x_j,v_{j2}=y_j,\;j=1,2,\ldots,n.$
Another reason the copula modeling approach offers significant flexibility is due to the availability of various copula functions. Among these, the AC family is a frequently used group of copula functions, and its n-dimensional AC function is defined as follows:
\begin{equation}
C_{\varphi}\left(u_1, u_2,\ldots,u_n\right)=\varphi\left(\varphi^{-1}\left(u_1\right)+\varphi^{-1}\left(u_2\right)+\cdots+\varphi^{-1}\left(u_n\right)\right), \varphi \in \Phi,
\end{equation} where Φ represents a class of function families: 
$\varphi: I \rightarrow[0, \infty]$, is a completely monotonic function, so that it satisfies
\begin{equation*}
\varphi(0)=1, \varphi(\infty)=0, (-1)^{k}\frac{d^k}{dt^k}\varphi(t)\geq0, k\in N,\;0 \lt t \lt 1,
\end{equation*} the function φ is recognized as the copula generator, and its inverse is denoted as φ −1, commonly defined as 
$\varphi^{-1}(u)=\inf\{t:\varphi(t)=u\}$.
General degree of dependence is an important feature for the copula. This can be explained by the link between Kendall’s tau(τ) and the copula, which is defined as follows.
Definition 2.2. Let variable (X, Y) and 
$\left(X^{\prime}, Y^{\prime}\right)$ be independently identically distributed (i.i.d.), the Kendall’s tau is defined as
\begin{equation*}\tau=P\left[\left(X-X^{\prime}\right)\left(Y-Y^{\prime}\right) \geq 0\right]-P\left[\left(X-X^{\prime}\right)\left(Y-Y^{\prime}\right) \leq 0\right].\end{equation*}If X and Y are AC dependent, according to Definition 2.2 and [Reference Nelsen37], the Kendall’s tau of (X, Y) can be written as
\begin{equation}
\tau=1+4 \int_0^1 \frac{\varphi(t)}{\varphi^{\prime}(t)}\mathrm{d}t,
\end{equation} where φ is generator of AC and 
$\varphi^{\prime}$ is the derivative of φ.
It should be noted that when given the n-dimensional AC function 
$C(\cdot)$, for any 
$r,s(r \lt s)\in\{1,2,\ldots,n\},$ we have
\begin{align}
C^{(r,s)}(u_r, u_s)&= C(\underbrace{1,\ldots,1}_{r-1},u_r,\underbrace{1,\ldots,1}_{s-r-1}, u_s,\underbrace{1,\ldots,1}_{n-s})\nonumber \\
&= \varphi(\varphi^{-1}\left(u_r\right)+\varphi^{-1}\left(u_s\right)+(n-2)\varphi^{-1}(1))\nonumber\\
&= \varphi(\varphi^{-1}\left(u_r\right)+\varphi^{-1}\left(u_s\right)).
\end{align} This implies that the resulting two-dimensional (2D) marginal distribution function 
$C^{(r,s)}\left(u_r, u_s\right)$ is still an AC function and has the same generator function φ.
Next, for all pairs 
$(X_{r}, X_{s})$ of n-AC with generator φ, the Kendall distribution of 
$(X_{r}, X_{s})$, where 
$r, s \in \{1, 2, \cdots, n\}$ and r ≠ s will be defined as
Definition 2.3. For all 
$r, s \in \{1, 2, \ldots, n\}$ and r ≠ s, let 
$\mathbf{X} \mapsto V^{r,s}=H(\mathbf{X})=C\left(U_r,U_s\right)$, where 
$\left(U_r,U_s\right) \sim C$. Then the Kendall distribution of 
$C(U_r,U_S)$ is
\begin{equation}
K^{\{r,s\}}(v)=P(V^{\{r,s\}} \lt v), v \in I.
\end{equation} For a bivariate AC with generator φ, Kendall’s tau 
$(\tau)$ can be calculated using (2).
Genest and Rivest [Reference Genest and Rivest15] demonstrated that the function 
$\varphi(t)$ can be derived from the Kendall distribution function 
$K(t) = \operatorname{Pr}\{C(U,V) \leq t\}$. Intriguingly, a relationship exists between the function 
$\varphi(t)$ and K(t), which can be expressed as follows
\begin{equation*}
K(t) =t-\frac{\varphi(t)}{\varphi^{\prime}(t)},
\end{equation*} and the function 
$\varphi(t)$ can be determined as
\begin{equation*}
\varphi(t) = \varphi\left(t_0\right) \exp \left[\int_{t_0}^t \frac{1}{z-K(z)} \mathrm{d}z \right],
\end{equation*} where 
$0 \lt t_0 \lt 1$ is a constant. The function K(t) is a crucial factor in identifying the function 
$\varphi(t)$ and therefore determines the dependence structure of the AC family.
2.2. Hierarchical Archimedean copulas
As is well known, simple multivariate ACs exhibit exchangeable, meaning that 
$\left(U_1, \ldots, U_k\right)=_{st}\left(U_{j_1}, \ldots, U_{j_k}\right)$, where 
$(j_1,j_2,\ldots,j_k)$ is any of the 
$n!$ permutations 
$(1,2,\ldots,k)$, and 
$=_{st}$ represents the distributional equivalence of the two random variables before and after. This assumption is often challenging to uphold in practical scenarios. A considerably more adaptable approach is furnished by HACs. A copula 
$C(\cdot)$ is classified as HAC if it adheres to the AC structure, where its arguments may be substituted with other HACs. When 
$C(\cdot)$ is recursively defined by (1) for k = 2, and further, with the arguments possibly permuted,
\begin{equation}
C\left(u_1, u_2,\ldots,u_k;\varphi_0,\varphi_1,\ldots,\varphi_{k-2} \right)=\varphi_{k-2}\left(\varphi_{k-2}^{-1}\left(C\left(u_1,\ldots,u_{k-1},1;\varphi_0,\varphi_1,\ldots,\varphi_{k-3}\right)+\varphi_{k-2}^{-1}\left(u_k\right)\right)\right),
\end{equation} then 
$C(\cdot)$ is called a HAC.
Remark 2.4. For 
$k\geq 3$ in (5), 
$C(\cdot)$ is called fully HACs with k − 1 hierarchies. Otherwise, 
$C(\cdot)$ is called partially HACs. Fully and partially HACs are summarized as HACs.
In the special case of fully nested copulas, the copula function is given by
\begin{align}
C(u_1, u_2,\ldots,u_k;\varphi_0,\varphi_1,\ldots,\varphi_{k-2} )&= \varphi_{k-2}(\varphi_{k-2}^{-1}\circ(\varphi _ {k -3 } [\cdots(\varphi _ {1 } ^ {- 1 } \circ \varphi _ {0 }[\varphi_0^{-1}(u_1)+\varphi_0^{-1}(u_2)]\nonumber\\
&+\varphi_1^{-1}(u_3))+\cdots+\varphi_{k-3}^{-1}(u_{k-1})])+\varphi_{k-2}^{-1}(u_k)).
\end{align} Illustrated the structure of HACs, we first illustrate it in an example. Consider 
$C^{f}=C_{\varphi_1}(C_{\varphi_2}(C_{\varphi_3}(u_1,u_2),u_3),u_4)$ and 
$C^{p_1}=C_{\varphi_1}(C_{\varphi_2}(u_1,u_2),C_{\varphi_3}(u_3,u_4))$ and 
$C^{p_2}=C_{\varphi_1}(C_{\varphi_2}(u_1,u_2,u_3),u_4)$ for three generators 
$\varphi_{1},\varphi_{2},\varphi_{3},$ its tree structure is depicted in Figure 1, where C f is called fully ACs and C p 1 and C p 2 are called partially ACs.
Figure 1. The structure tree-based complete and partial HACs.
The greatest strength of HAC lies in its capability to capture asymmetric dependence. Unlike conventional ACs, HACs delineate the entire dependence structure in a recursive manner. Note that the generators φi can belong to the same family with varying parameters, in which case HACs are referred to as homogeneous HACs. Alternatively, for greater flexibility, the generators can be derived from different families, in which case HAC is known as heterogeneous HACs.
2.3. Model description
Suppose that a system contains n parallel components which suffers from a common external stress. The component’s strengths 
$X_{1},X_{2},\ldots,X_{n}$ are dependent non-identically distributed random variables with cdf 
$F_{i}(\cdot)$, respectively. The stress Y is a random variables with cdf 
$G(\cdot)$, which is dependent with each Xi. Then the reliability of stress-strength model with parallel structure is given by
\begin{eqnarray*}
R&=&P(\max(X_{1},X_{2},\ldots,X_{n}) \gt Y)=P(X_{n:n} \gt Y).
\end{eqnarray*} Next, we will use the hierarchical copula to describe the dependence structure among variables. Without loss of generality, suppose that the dependence structure of Y and 
$X_{n:n}$ is represented by selected 2D copula 
$ C(u,v,\theta),$ and the dependence structure of 
$X_{1},X_{2},\ldots,X_{n}$ is represented by selected n-dimensional copula 
$C_{1}(u_1,u_2,\ldots,u_n,\theta_1).$ Therefore, R can then be written as
\begin{equation}
R=P(X_{n:n} \gt Y)= \int_{0}^{\infty}\frac{\partial C(F_{X_{n:n}}(x),G(x))}{\partial F_{X_{n:n}}(x)}{\mathrm{d}}F_{X_{n:n}}(x),
\end{equation} where 
$F_{X_{n:n}}(x)$ represents the CDF of 
$X_{n:n}$, as given by
\begin{equation}
F_{X_{n:n}}(x)=C_1(F_{X_1}(x),F_{X}(x_2),\ldots,F_{X_{n}}(x),\theta_1).
\end{equation} Let 
$H(x,y)$ be the joint cdf of 
$(X_{n:n},Y),$ then
\begin{align}
H(x,y)&= P(X_{n:n}\leq x,Y\leq y)\nonumber\\
&= C(C_1(F_{X_1}(x),F_{X_2}(x),\ldots,F_{X_{n}}(x),\theta_1),G(y),\theta).
\end{align} Note that 
$H(x,y)$ is actually a hierarchical copula that satisfies a specific structure, the structure tree-based hierarchical of copulas in Figure 2. The hierarchical copula is a powerful tool for modeling the dependence structure of multidimensional data. It involves a series of nested conditional copulas, each linking a subset of variables in the data. This allows for a more flexible and realistic modeling of complex dependence structures, since it can capture both global and local dependencies among the variables.
Figure 2. The structure tree-based 2-HAC.
Suppose that M items are put into a life testing experiment and the observed data are 
$X_{i1},X_{i2},\ldots,X_{in}$ and 
$Y_{i},\;i=1,2,\ldots,M.$ The observed values of the strength and stress random variables are designed as follows:
\begin{equation*}
\left(\begin{array}{cccc}
X_{11} & X_{12} & \cdots & X_{1 n} \\
\vdots & \vdots & \vdots & \vdots \\
X_{M 1} & X_{M 2} & \cdots & X_{M n}
\end{array}\right) \text {and }\left(\begin{array}{c}
Y_1 \\
\vdots \\
Y_M
\end{array}\right).
\end{equation*}3. Semi-parametric estimate
According to (7), if we obtain 
$\hat{F}_{X_{n:n}}(x)$, 
$\hat{G}(x)$, 
$\hat{\theta}_1$, and 
$\hat{\theta}$, the estimate of R is given by
\begin{equation}
\hat{R} = \int_{0}^{\infty}\frac{\partial C(\hat{F}_{X_{n:n}}(x),\hat{G}(x);\hat{\theta})}{\partial \hat{F}_{X_{n:n}}(x)}d\hat{F}_{X_{n:n}}(x).
\end{equation} It is important to note that HAC arises in the aforementioned solution. However, HAC is not directly observable, necessitating the creation of a semi-parametric HAC copula estimator in two stages: initially estimating the marginals and subsequently constructing the copula estimator using these estimated marginals. Due to the intricacies involved, in the subsequent phases, we adopt a stepwise estimation approach to estimate the pertinent parameters. The following Step 1–Step 2 can obtain 
$\hat{G}(x),$ 
$\hat{\theta}_1$ and 
$\hat{F}_{X_{n:n}}(x).$
Step 1. We derive nonparametric estimates of the marginal distributions of stress and strength separately using the empirical distribution function (EDF). Based on observed data 
$X_{i1},X_{i2},\ldots,X_{in}$ and 
$Y_{i},\;i=1,2,\ldots,M,$ using the marginal empirical method, the estimate of 
${F}_{X_{j}}(x_{ij})$ and 
$G(y_i),\;i=1,2,\ldots,M,\;j=1,2,\ldots,n$ are given by
\begin{equation}
\hat{u}_{ij}= \frac{M \times E F_{X_{j}}\left(x_{ij}\right)}{(M+1)}=\frac{R^{u}_{ij}}{M+1},
\end{equation}and
\begin{equation}
\hat{v}_{i}=\frac{M \times E F_{Y}\left(y_{i}\right)}{(M+1)}=\frac{R^{v}_{i}}{M+1},
\end{equation} respectively, where EFX and EFY are the marginal empirical cdf of X and Y, respectively, 
$R^{u}_{ij}$ denotes the rank of Xij among 
$X_{1j},X_{2j},\ldots,X_{Mj}$, and 
$R^{v}_{i}$ denotes the rank of Yi among 
$Y_{1},Y_{2},\ldots,Y_{M}$.
Step 2. Estimation of the dependence parameter θ 1
Generally speaking, the estimation of Copula dependence parameters is based on likelihood methods. However, for an n-dimensional Copula 
$C_1(\cdot,\theta_1)$, the likelihood function is very complex. Therefore, in this part, we consider methods based on Kendall’s τ and Bernstein polynomials to separately estimate the dependence parameter θ 1.
Step 2.1. Based on Kendall’s tau of 
$C_1(\cdot,\theta_1)$:
We employ the inversion of Kendall’s tau to estimate the dependence parameter θ 1. It was introduced by [Reference Kendall and Smith25] as a coefficient of agreement among 
$n\geq 2$ rankings. The Kendall’s tau of 
$\left(X_{1}, X_{2},\ldots,X_{n}\right)$ may be defined as the average value of Kendall’s tau taken over all possible pairs 
$(X_r, X_s)$ with 
$r,s={1,2,\ldots,n}$ of 
$\left(X_{1}, X_{2},\ldots,X_{n}\right)$
\begin{align}
\tau_{n}&=\frac{1}{n(n-1)}\sum_{r\neq s}\tau_{2}(X_r,X_s)\nonumber\\
&=\frac{1}{n(n-1)}\sum_{r\neq s}\left(4\int_{[0,1]^2}C_{1}(u_r, u_s;\theta_{1})dC_{1}(u_r, u_s;\theta_{1})-1\right).
\end{align}Further, the consistent estimator of τn can be expressed as
\begin{align}
\hat{\tau}_{n,M}&=\frac{1}{n(n-1)}\sum_{r\neq s}\hat{\tau}_{2,M}(X_r,X_s)\nonumber\\
&=\frac{1}{n(n-1)}\sum_{r\neq s}\left(\frac{4}{M(M-1)}\sum_{i\neq j}\mathbb{I}(X_{ir}\leq X_{jr},X_{is}\leq X_{js})-1\right),
\end{align} where 
$i,j=1,2,\ldots,M,r,s=1,2,\ldots,n$ and 
$\mathbb{I}(\cdot)$ is the indicative function.
Therefore, 
$\hat{\theta}_1$ can be obtained from the following nonlinear equation,
\begin{equation}
\theta_{\tau,M}=\xi^{(-1)}_{n}(\hat{\tau}_{n,M}),
\end{equation} where 
$\xi_n:\theta_1\mapsto\tau_n.$
Based on [Reference Genest, Nešlehová and Ben Ghorbal14], it is known that when the mapping ξn is both invertible and differentiable, the consistency and asymptotic normality of the estimator θn can be established by observing that 
$\hat{\tau}_{n,M}$ converges to a U-statistic that is unbiased and Gaussian as 
$n \rightarrow \infty$.
Step 2.2. Based on Bernstein estimate of Kendall distribution function of 
$C_1(\cdot,\theta_1),$
For over all possible pairs 
$(X_r, X_s)$ with 
$r,s=\{1,2,\ldots,n\}$ of 
$\left(X_{1}, X_{2},\ldots,X_{n}\right)$, the Kendall distribution of 
$(X_r, X_s)$ is given by
\begin{equation*}K^{\{r,s\}}(v)=P(C(U_r,U_s) \lt v), v \in {I},\end{equation*} where Ur and Us are random variables defined on the interval [0,1] following a uniform distribution, and the joint distribution of the random vector 
$(U_r,U_s)$ is also denoted by 
$C_1(u_r,u_s)$.
Furthermore, we can obtain the consistent estimator of 
$K^{\{r,s\}}(v)$
\begin{equation}
\hat{K}^{r,s}_M(v)=\frac{1}{M} \sum_{i=1}^M \left\{V^{r,s}_i \leq v\right\},
\end{equation} where 
$V^{r,s}_i=\frac{1}{M-1}\sum_{j=1}^{M}\mathbb{I}(X_{ir} \lt X_{jr},X_{is} \lt X_{js}),i=1,2,\ldots, M,$ and 
$(X_{ir},X_{is})|_{r,s=\{1,2,\ldots,n\}}$ is over all possible pairs of 
$(X_{i1},X_{i2},\ldots,X_{in}).$
Although EDF is a reliable estimator of the distribution function, it may not be suitable for estimating a continuous distribution function due to its jump discontinuities. To address this issue, a smooth estimator of the distribution function was introduced by [Reference Babu and Chaubey3], employing Bernstein polynomials. This novel estimator offers a continuous approximation of the conventional EDF.
Now, based on the pseudo-observations data, we use Bernstein method to estimate the univariate distribution function 
$K^{\{r,s\}}(v)$. The Bernstein estimator of order m > 0 for the Kendall distribution function 
$K^{\{r,s\}}(v)$ is expressed as the following equation:
\begin{equation}
\hat{K}^{\{r,s\}}_{m, M}(v)=\sum_{k=0}^m \hat{K}^{\{r,s\}}_M(k / m) P_{k, m}(v),
\end{equation} where 
$\hat{K}^{\{r,s\}}_n$ is the EDF of 
$K^{\{r,s\}}(t)$ (as defined in (16)) and 
$P_{k,m}(t)$ denotes the binomial probability with parameters m and k.
It’s important to note that the choice of the order m for the Bernstein polynomial estimator is crucial. According to the recommendation in [Reference Babu and Chaubey3], the order is typically chosen as 
$m = \frac{M}{\log(M)}$.
Using the Bernstein estimator of 
$K^{\{r,s\}}(v)$, we can obtain an estimator for 
${\tau}_{2}(X_r,X_s)=3-4\int_{0}^{1}{K}^{\{r,s\}}(t){\mathrm{d}}t$ as follows:
\begin{equation}
\hat{\tau}^{B\{r,s\}}_{2,m,M}=3-4\int_{0}^{1}\hat{K}^{\{r,s\}}_{m,M}(t){\mathrm{d}}t=3-4\int_{0}^{1}\left(\sum_{k=0}^m \hat{K}^{\{r,s\}}_M(k / m) P_{k, m}(t)\right){\mathrm{d}}t.
\end{equation}Furthermore, we obtain the estimate of τn as
\begin{align}
\hat{\tau}^{B}_{n,m,M}&= \frac{1}{n(n-1)}\sum_{r\neq s}\hat{\tau}^{B\{r,s\}}_{2,m,M}\nonumber\\
&= \frac{1}{n(n-1)}\sum_{r\neq s}\left[3-4\int_{0}^{1}\left(\sum_{k=0}^m \hat{K}^{\{r,s\}}_M\left(\frac{k}{m}\right) P_{k, m}(t)\right){\mathrm{d}}t\right].
\end{align} Therefore, we can obtain 
$\hat{\theta}_1$ from the following nonlinear equation:
\begin{equation}
\theta_{1}= \xi^{(-1)}_{n}(\hat{\tau}^{B}_{n,m,M}),
\end{equation} where 
$\xi_n:\theta_1\mapsto\tau_n.$
Next, we prove that 
$\hat{\tau}^{B}_{n,m,M}$ is consistent and asymptotically unbiased. Before presenting the main results, we first introduce some symbols and useful lemmas.
For a bounded function F defined on 
$[0,1]$, we denote 
$\| F \| = \sup_{x \in [0,1]} |F(x)|$. Let 
$\mathbf{X}_{i}=\left(X_{i1}, X_{i2},\ldots,X_{in}\right), i=1,2, \ldots, M$ be a sample of capacity M for an n-dimensional random intensity vector 
$\mathbf{X}=\left(X_1, X_2,\ldots,X_{n}\right)$ defined on 
$[0,1]$.
Lemma 3.1. Let 
$ \hat{\tau}^{B\{r,s\}}_{2,m}=3-4\int_{0}^{1}\hat{K}^{\{r,s\}}_{m}(t){\mathrm{d}}t,$ where 
$\hat{K}^{\{r,s\}}_m(v)=\sum_{k=0}^m {K}^{\{r,s\}}(k / m) P_{k, m}(v)$. Then, as 
$m\rightarrow\infty$,
\begin{equation*}\hat{\tau}^{B\{r,s\}}_{2,m}\rightarrow{\tau}^{\{r,s\}}_{2}.\end{equation*}Proof. Note that
\begin{align}
\Vert\hat{\tau}^{B\{r,s\}}_{2,m}-{\tau}^{\{r,s\}}_{2}\Vert&= \left\Vert3-4\int_{0}^{1}\hat{K}^{\{r,s\}}_{m}(t){\mathrm{d}}t-3+4\int_{0}^{1}{K}^{\{r,s\}}(t){\mathrm{d}}t\right\Vert\nonumber\\
&\leq 4\int_{0}^{1}\Vert{K}^{\{r,s\}}(t)-\hat{K}^{\{r,s\}}_{m}(t)\Vert {\mathrm{d}}t.
\end{align} Since 
${K}^{\{r,s\}}(x)$ is a continuous bounded function on 
$[0,1]$, according to Theorem 1 ([Reference Feller13], Chapter VII.2), for all 
$t\in[0,1]$, as 
$m\rightarrow\infty$,
\begin{equation*}\Vert\hat{K}^{\{r,s\}}_{m}(t)dt-{K}^{\{r,s\}}(t)\Vert\rightarrow 0\end{equation*}almost everywhere. Thus, the lemma is proved.
Now, we will show that 
$\hat{\tau}^{B}_{n,m,M}$ is consistent.
Theorem 3.2. Let 
$\mathbf{X}_{1},\mathbf{X}_{2},\ldots,\mathbf{X}_{M}$ be a sequence of 
$M\geq2$ mutually independent observations from a continuous n-variate distribution with underlying copula 
$C_1(\cdot,\theta_1)$, then,
\begin{equation*}\Vert\hat{\tau}^{B}_{n,m,M}-{\tau}^{B}_{n}\Vert\rightarrow 0.\;\text{a.s.}\;m,M\rightarrow\infty.\end{equation*}Proof. We first note that
\begin{align*}
\Vert\hat{\tau}^{B}_{n,m,M}-{\tau}^{B}_{n}\Vert&=\left\Vert\frac{1}{n(n-1)}\sum_{r\neq s}\hat{\tau}^{B\{r,s\}}_{2,m,M}-\frac{1}{n(n-1)}\sum_{r\neq s}{\tau}^{B\{r,s\}}_{2}\right\Vert\nonumber\\
&\leq\frac{1}{n(n-1)}\sum_{r\neq s} \Vert\hat{\tau}^{B\{r,s\}}_{2,m,M}-{\tau}^{B\{r,s\}}_{2}\Vert.\nonumber
\end{align*} By introducing 
$\hat{\tau}^{B\{r,s\}}_{2,m}$, it follows from the triangle inequality that
\begin{equation} \Vert\hat{\tau}^{B\{r,s\}}_{2,m,M}-\hat{\tau}^{B\{r,s\}}_{2}\Vert\leq\Vert\hat{\tau}^{B\{r,s\}}_{2,m,M}-\hat{\tau}^{B\{r,s\}}_{2,m}\Vert+\Vert\hat{\tau}^{B\{r,s\}}_{2,m}-{\tau}^{B\{r,s\}}_{2}\Vert.
\end{equation}To show that the estimator is consistent, it needs to show that the sum of the two terms on the RHS of the inequality goes to zero as m and M approach infinity. Recall (18) and (19), we have
\begin{align}
\Vert\hat{\tau}^{B\{r,s\}}_{2,m,M}-\hat{\tau}^{B\{r,s\}}_{2,m}\Vert&\leq4\int_{0}^{1}\Vert{\hat{K}}^{\{r,s\}}_{m}(t)-\hat{K}^{\{r,s\}}_{m,M}(t)\Vert {\mathrm{d}}t\nonumber\\
&\leq4\int_{0}^{1}\left\Vert\sum_{k=0}^m \left({K}^{\{r,s\}}(k/m)-\hat{K}^{\{r,s\}}_M(k/m)\right)\right\Vert P_{k, m}(t) {\mathrm{d}}t\nonumber\\
&\leq4\left\Vert\sum_{k=0}^m \left({K}^{\{r,s\}}(k/m)-\hat{K}^{\{r,s\}}_M(k/m)\right)\right\Vert\nonumber\\
&\leq 4\max_{0\leq k\leq m}\left\vert{K}^{\{r,s\}}(k/m)-\hat{K}^{\{r,s\}}_M(k/m)\right\vert\nonumber\\
&\leq4\Vert{K}^{\{r,s\}}-\hat{K}^{\{r,s\}}_M\Vert.
\end{align} Since by Glivenko–Cantelli theorem, 
$\Vert{K}^{\{r,s\}}-\hat{K}^{\{r,s\}}_M\Vert\rightarrow 0$ almost surely as 
$M\rightarrow \infty$, the result follows from (22), (23), and Lemma 3.1.
Combine Theorem 3.2, (20), and Delta method, we have
\begin{equation}
\sup _{\theta_1\in \Theta_1}\left|\hat{\theta}_{1}-{\theta}_{1}\right| \rightarrow 0 .
\end{equation} The following theorem shows that 
$\hat{\tau}^{B}_{n,m,M}$ is asymptotically unbiased.
Theorem 3.3. Let 
$\mathbf{X}_{1},\mathbf{X}_{2},\ldots,\mathbf{X}_{M}$ be a sequence of 
$M\geq2$ mutually independent observations from a continuous n-variate distribution with underlying copula 
$C_1(\cdot,\theta_1)$, then, for 
$m\rightarrow\infty$
\begin{equation*}E(\hat{\tau}^{B}_{n,m,M})\rightarrow {\tau}^{B}_{n}.\end{equation*}Proof. According to (19), we have
\begin{eqnarray*}
E(\hat{\tau}^{B}_{n,m,M})&=&E\left(\frac{1}{n(n-1)}\sum_{r\neq s}\hat{\tau}^{B\{r,s\}}_{2,m,M}\right)\nonumber\\
&=&\frac{1}{n(n-1)}\sum_{r\neq s}\left[3-4\int_{0}^{1}\left(\sum_{k=0}^m E(\hat{K}^{\{r,s\}}_M(k / m)) P_{k, m}(t)\right){\mathrm{d}}t\right]\nonumber\\
&=&\frac{1}{n(n-1)}\sum_{r\neq s}\left[3-4\int_{0}^{1}\left(\sum_{k=0}^m {K}^{\{r,s\}}(k / m)) P_{k, m}(t)\right){\mathrm{d}}t\right]\nonumber\\
&=&\frac{1}{n(n-1)}\sum_{r\neq s}\hat{\tau}^{B\{r,s\}}_{2,m},\nonumber
\end{eqnarray*}the result follows from Lemma 3.1, this completes the proof.
Hence, based on 
$\hat{\theta}_1$ and 
$\hat{u}_{ij},i=1,2,\ldots,M,j=1,2,\ldots,n,$ the estimate of 
$F_{X_{n:n}}(x_1,x_2,\ldots,x_n)$ can be given by
\begin{equation}
\hat{F}_{X_{n:n}^{i}}={C}_1(\hat{u}_{i1},\hat{u}_{i2},\ldots,\hat{u}_{in},\hat{\theta}_1):=\hat{u}_{X_{n:n}^{i}},\;i=1,2,\ldots,M,
\end{equation} where 
$X_{n:n}^{i}=\mathrm{max}(X_{i1},X_{i2},\ldots,X_{in}).$
Step 3. Estimation of the dependence parameter θ:
Below, based on 
$\hat{F}_{X{n:n}}(x)$, 
$\hat{G}(x)$, 
$\hat{\theta}_1$, and copula 
$C(\cdot)$, we utilize the method of moments to obtain an estimate 
$\hat{\theta}$. Let 
$\mathbf{X}^{*}_{i}=(X^{i}_{n:n},Y_{i})$ for 
$i=1,2,\ldots,M$ be a sample of 2D strength-stress, and 
$\mathbf{X}^{*}=(X_{n:n},Y)$. According to Step 2.1, we can obtain an estimate of 
$\hat{\theta}$ as follows:
\begin{equation}
\theta = \xi^{(-1)}_{2}(\hat{\tau}_{2,M}),
\end{equation} where 
$\xi_2:\theta\rightarrow \tau_2$, and 
$\hat{\tau}_{2,M}$ satisfies:
\begin{equation}
\hat{\tau}_{2,M} = \frac{4}{M(M-1)}\left(\sum_{i\neq j}\mathbb{I}(X^{i}_{n:n}\leq X^{j}_{n:n},Y_{i}\leq Y_{j})-1\right),
\end{equation} where 
$\mathbb{I}(\cdot)$ denotes the indicator function.
Thus, we have obtained estimates for 
${F}_{X_{j}}(x_{ij})$ and 
$G(y_i)$ for 
$i=1,2,\ldots,M$, 
$j=1,2,\ldots,n$, as well as estimates for θ 1 and θ. Thus, for Eq. (7), the estimated value of R, denoted as 
$\hat{R}$, can be expressed as:
\begin{equation}
\hat{R}=\int_{0}^{\infty}\frac{\partial C(\hat{F}_{X_{n:n}}(x),\hat{G}(x),\hat{\theta})}{\partial \hat{F}_{X_{n:n}}(x)}{\mathrm{d}}\hat{F}_{X_{n:n}}(x).
\end{equation} Next, we are going to prove that 
$\hat{R}$ is a consistent estimator of R. Before the proof, let’s start with the following useful conditions.
C1: Suppose that strength random vector 
$\left(X_1,X_2, \ldots, X_n\right)$ is a nonnegative continuous random variable vector on the set 
$D=\left[0, \infty\right]^{n}$. The joint distribution function of 
$\max\left(X_1, X_2,\ldots, X_n\right)$, denoted by 
$F_{X_{n:n}}\left(x\right)$, is continuous.
C2: Suppose that the stress random variable Y is a nonnegative continuous random vector on the set 
$[0, \infty]$. The distribution function of Y, denoted by G(y), is continuous.
C3: Suppose that the dependent parametric 
$\theta_1(\theta)$ is defined in the set Θ.
C4: Suppose that 
$\frac{\partial C(u,v)}{\partial v}$ exist everywhere and continuous on 
$(0,1)^2$.
We are now ready for the following useful theorem.
Theorem 3.4. If conditions C1–C4 hold. Then 
$\hat{R}$ is a consistent estimator of R.
Proof. According to (28), we have
\begin{align}
\Vert\hat{R}-R\Vert&= \left\Vert\int_{0}^{\infty}\frac{\partial C(\hat{F}_{X_{n:n}}(x),\hat{G}(x);\hat{\theta})}{\partial \hat{F}_{X_{n:n}}(x)}{\mathrm{d}}\hat{F}_{X_{n:n}}(x)-\int_{0}^{\infty}\frac{\partial C({F}_{X_{n:n}}(x),{G}(x);{\theta})}{\partial {F}_{X_{n:n}}(x)}{\mathrm{d}}{F}_{X_{n:n}}(x)\right\Vert\nonumber\\
&\leq \left\Vert\int_{0}^{\infty}\left[\frac{\partial C(\hat{F}_{X_{n:n}}(x),\hat{G}(x);{\theta})}{\partial \hat{F}_{X_{n:n}}(x)}-\frac{\partial C({F}_{X_{n:n}}(x),{G}(x);{\theta})}{\partial {F}_{X_{n:n}}(x)}\right]{\mathrm{d}}{F}_{X_{n:n}}(x)\right\Vert\nonumber\\
&+\left\Vert\int_{0}^{\infty}\frac{\partial C(\hat{F}_{X_{n:n}}(x),\hat{G}(x),\hat{\theta})}{\partial \hat{F}_{X_{n:n}}(x)}{\mathrm{d}}[(\hat{F}_{X_{n:n}}(x)-{F}_{X_{n:n}}(x))]\right\Vert.
\end{align}From (8), we have
\begin{equation*}
\left\Vert\hat{F}_{X_{n:n}}-F_{X_{n:n}}\right\Vert=\left\Vert{C}(\hat{F_1},\hat{F_2},\ldots,\hat{F_n};\hat{\theta}_1)-{C}(F_1,F_2,\ldots,F_n;\theta_1)\right\Vert.
\end{equation*}Hence
\begin{equation*}
\left|\hat{F}_{X_{n:n}}(x;\hat{\theta}_1)-F_{X_{n:n}}(x;\theta_1)\right|=\left|{C}(\hat{F_1}(x),\hat{F_2}(x),\ldots,\hat{F_n}(x);\hat{\theta}_1)-{C}(F_1(x),F_2(x),\ldots,F_n(x);\theta_1)\right|.\nonumber
\end{equation*} Note that 
${C}_1(u_1,u_2,\ldots,u_n)$ is a continuous function defined on 
$[0,1]^n$. According to the properties of copula functions, the function 
$C_{u}(u,v,\theta)=\frac{\partial C(u,v,\theta)}{\partial u}$ is continuous on 
$[0,1]^2$. Thus, by the continuous mapping theorem and the convergence properties 
$\sup _{x_i\in D}\left|\hat{F}_{X_i}(x)-F_{X_{i}}(x)\right| \rightarrow 0$, 
$\sup _{x\in D}\left|\hat{G}(x)-G(x)\right| \rightarrow 0$, and (24), we can conclude that:
\begin{equation*}\sup _{x \in D}\left|\hat{F}_{X_{n:n}}(x;\hat{\theta}_1)-F_{X_{n:n}}(x;\theta_1)\right| \rightarrow 0,\end{equation*}and
\begin{equation*}\sup _{x \in D}\left|C_{u}(\hat{F}_{X_{n:n}}(x),\hat{G}(x),\hat{\theta})-C_{u}({F}_{X_{n:n}}(x),{G}(x),{\theta})\right|= 0.\end{equation*}This theorem is thus proved.
4. Numerical results
In this section, we investigate the finite-sample behavior of semi-parametric estimation for multicomponent stress-strength model under HAC. In the simulation studies, we consider four AC families with a wide range of parameter values are considered.
Case 1: the Clayton copula (C copula) family defined as
\begin{equation*}
C_C(u_1,u_2,\ldots,u_n)=\left(\sum_{i=1}^{n}u^{-\theta}_i-n+1\right)^{-1 / \theta},
\end{equation*} where 
$\theta \in (0, \infty).$
Case 2: the Gumbel copula (G copula) family defined as
\begin{equation*}
C_G(u_1,u_2,\ldots,u_n)=\exp \left(-\left(\sum_{i=1}^{n}(-\ln (u_i))^\theta\right)^{1 / \theta}\right),
\end{equation*} where 
$\theta \in[1, \infty)$.
Case 3: the Frank copula (F copula) defined as
\begin{equation*}C_{F}\left(u_1, \ldots, u_n\right)=-\frac{1}{\theta} \log \left(1+\frac{\prod_{i=1}^n\left(\exp \left(-\theta u_i\right)-1\right)}{(\exp (-\theta)-1)^{n-1}}\right),\end{equation*} where 
$\theta \in[-\infty,\infty).$
Case 4: the Joe copula (J copula) defined as
\begin{equation*}
C_J\left(u_1, \ldots, u_n\right)=1-\left(1-\prod_{i=1}^{n}\left[1-\left(1-u_i\right)^\theta\right]\right)^{1 / \theta},
\end{equation*} where 
$\theta\in[1,\infty)$.
The copula generators and Kendall’s τ corresponding to the well-known cases of Clayton copula, Gumbel copula, Frank copula, and Joe copula are provided in Table 2.
Table 2. The copula generators and Kendall’s τ of Clayton, Gumbel, Frank, and Joe copula.
According to the conclusion in [Reference McNeil32], as long as C 1 and C belong to the same family and 
$\theta_1 \lt \theta$, the constructed hierarchical copula remains a copula. Therefore, when 
$\theta_1 \lt \theta$, the C-C HAC, G-G HAC, F-F HAC, and J-J HAC are all appropriate copulas, where the C-C HAC, G-G HAC, F-F HAC, and J-J HAC satisfy the following condition:
\begin{equation*}
H(u_1,u_2,\ldots,u_n,v) = C(C_1(u_1,u_2,\ldots,u_n;\theta_1),v;\theta).
\end{equation*}4.1. Simulation studies
Now, we will use Monte Carlo (Markov chain) simulations to validate the described estimation methods.
Consider M identical multicomponent products (from the same manufacturer and manufacturing process). Assume the products are composed of n different components in parallel and are subjected to a single stress. The reliability of such a product is given by:
\begin{equation*}
R=P(\max(X_1,X_2,\ldots,X_{n}) \gt Y),
\end{equation*} where 
$X_i|_{i=1,2,\ldots,n}$ represents the strength of the ith component and Y represents the stress on the product. Without loss of generality, assume Xi follows an exponential distribution with parameter λi, and Y follows an exponential distribution with parameter λ. To analyze the reliability of these M multicomponent products, four types of test products with 
$n=2,3,4,8$ are selected, denoted as Model 1, Model 2, Model 3, and Model 4. Set the true model parameters as
$(\lambda_1,\lambda_2,\lambda)=(1,2,3),$ 
$(\lambda_1,\lambda_2,\lambda_3,\lambda)=(1,2,3,4),$ 
$(\lambda_1,\lambda_2,\lambda_3,\lambda_4,\lambda)=(1,2,3,4,5)$ 
$(\lambda_1,\lambda_2,\lambda_3,\lambda_4,\lambda_5,\lambda_6,\lambda_7,\lambda_8,\lambda)=(1,2,3,4,5,6,7,8,9)$, and the dependency relationships between components and between strength and stress are characterized by C-C HAC, G-G HAC, F-F HAC, and J-J HAC. Additionally, based on the reviewers’ suggestions, we have also considered the independent scenario. Three sample sizes 
$M=\{50,100,200\}$ are chosen for the experimental samples.
Although the focus of this subsection is on the dependency parameters θ and θ 1 in HAC, the meanings and ranges of the dependency parameters vary across different copula functions. Therefore, they cannot be directly compared, whereas Kendall’s rank correlation coefficient can be directly compared across different models. Additionally, functions such as copClayton@iTau, copGumbel@iTau, copFrank@iTau, and copJoe@iTau in the R software can conveniently convert Kendall’s τ to dependency parameters. Hence, for consistent comparison of dependency parameters, we will convert all dependency parameters to Kendall’s τ. In this section, we set 
$(\tau_1,\tau)=(0.2,0.5)$, which evidently satisfies 
$\theta_1 \lt \theta$ in C-C HAC, G-G HAC, F-F HAC, and J-J HAC.
Based on the above settings, we first present an algorithm (Algorithm 4.1) for generating simulated data for a multicomponent stress-strength model with a sample size of M.
Based on the pseudo-observation data 
$\mathscr{D}$ obtained from Algorithm 4.1, the reliability R and the dependent parameters θ 1 and θ of the four types of multicomponent stress-strength models 
$M_1, M_2, M_3, M_4$ can be estimated using the methods proposed in this chapter. In order to evaluate the effectiveness of the estimation methods, we conducted a total of N = 1000 repeated sampling experiments, and then assessed the estimation performance using bias and Mean Squared Error (MSE). The results are shown in Tables 3–7. In these tables, subscript 
${1K}$ indicates estimation of dependent parameter θ 1 based on Kendall’s τ, subscript 
${1B}$ indicates estimation of dependent parameter θ 1 based on the Bernstein method, subscript KK indicates estimation of dependent parameter θ using Kendall’s τ method on the basis of estimating θ 1 using Kendall’s τ, subscript KB indicates estimation of dependent parameter θ using the Bernstein method on the basis of estimating θ 1 using Kendall’s τ, subscript BK indicates estimation of dependent parameter θ using Kendall’s τ method on the basis of estimating θ 1 using the Bernstein method, and subscript BB indicates estimation of dependent parameter θ using the Bernstein method on the basis of estimating θ 1 using the Bernstein method.
Table 3. The bias and MSE of dependent parameters 
$(\theta,\theta_1)$ and MSS R under C-C copula.
Table 4. The Bias and MSE of dependent parameters 
$(\theta,\theta_1)$ and MSS R under G-G copula.
Table 5. The bias and MSE of dependent parameters 
$(\theta,\theta_1)$ and MSS R under F-F copula.
Table 6. The bias and MSE of dependent parameters 
$(\theta,\theta_1)$ and MSS R under J-J copula.
Table 7. The bias and MSE of R under independence case.
Based on the results presented in Tables 3, 4, 5, 6, and 7, all estimation methods exhibit satisfactory performance regardless of whether the C-C HAC, G-G HAC, F-F HAC, or J-J HAC method are used to describe hierarchical dependence or independence case. It is worth noting that these methods maintain excellent efficiency even with limited sample sizes. As the sample size M increases gradually, most cases demonstrate a gradual improvement in estimation accuracy. When the sample size remains constant, the increase in the number of strength variables generally leads to improved estimation accuracy. Specifically, in the estimation of reliability R, the accuracy improves significantly with the increase in the number of strength variables, as evidenced by the trends in Figures 3, 4, and 5. Additionally, it can be clearly observed that RKB consistently outperforms RKK when the sample size remains constant. Overall, these findings demonstrate that selecting the RKB method can yield optimal estimates of R. These analyses not only enrich our understanding of hierarchical dependence modeling but also provide theoretical support for decision-making in practical applications.
Figure 3. The bias and MSE of R under G-G HACs with M = 50.
Figure 4. The bias and MSE of R under G-G HACs with M = 100.
Figure 5. The bias and MSE of R under G-G HACs with M = 200.
Furthermore, based on Model 1, Model 2, Model 3, and Model 4, Table 8 presents the discrepancies when the model is dependent but the dependence structure is ignored, that is independent copula, with the discrepancies represented by the percentage error (PE). The calculation formula is as follows:
\begin{equation*}PE=\frac{R_{dependent}-R_{independent}}{R_{independent}}\times100\%.\end{equation*}Table 8. The discrepancies in the reliability estimates of the model when the dependence structure is ignored.
4.2. Goodness-of-fit test and empirical study
In this section, We carry out an empirical study on four real datasets and demonstrate the above mentioned methods can be applied in practice. These four dates are as follows:
Data 1: 
$693.73, 704.66, 323.83, 778.17, 123.06, 637.66, 383.43, 151.48, 108.94, 50.16, 671.49,$
$ 183.16, 727.23, 257.44, 291.27, 101.15, 376.42, 163.40, 141.38, 700.74, 262.90, 353.24, 422.11,$
$ 43.93, 590.48, 212.13, 303.90, 506.60, 530.55, 177.25$
Data 2: 
$71.46, 419.02, 284.64, 585.57, 456.60, 688.16, 662.66, 113.85, 187.85, 45.58, 578.62,$
$756.70, 594.29, 166.49, 707.36, 99.72, 765.14, 187.13, 145.96, 350.70, 547.44, 116.99, 375.81,$
$119.86, 581.60, 48.01, 200.16, 36.75, 244.53, 83.55$
Data 3: 
$6.53, 7, 10.42, 14.48, 16, 10, 22.70, 34, 41.55, 42, 45.28, 49.40, 53.62, 63, 83, 84, 91, 108,$
$ 112, 129, 133, 133, 139, 140, 140, 146, 149, 154, 157, 160, 160, 165, 173, 176, 218, 225, 241, 248,$
$ 273, 277, 297, 405, 417, 420, 440, 523, 583, 594, 1101, 1146, 1417$
Data 4: 
$12.20, 23.56, 23.74, 25.87, 31.98, 37, 41.35, 47.38, 55.46, 58.36, 63.47, 68.46, 74.48, 78.26,$
$ 81.43, 84, 92, 94, 110, 112, 119, 127, 130, 133, 140, 146, 155, 159, 173, 179, 194, 195, 209, 249, 281,$
$319, 339, 432, 469, 519, 633, 725, 817, 1776$
Remark 4.1. The Data 1 and Data 2 have been analyzed by [Reference Saraçoğlu, Kinaci and Kundu44] for stress-strength reliability with exponential distribution. Data 3 and Data 4 have been discussed by [Reference Iranmanesh, Fathi Vajargah and Hasanzadeh19] for estimations of stress-strength reliability with inverted gamma distribution. All of them have been analyzed by [Reference Bai, Li, Balakrishnan and He4] for dependent stress-strength reliability of multistate system.
To investigate dependencies more effectively, it is essential to have an equal sample size for each dataset. Currently, Data 1 and Data 2 have a sample size of 30, while Data 3 and Data 4 have sample sizes greater than 30. To address this, we utilized the sample() function in R to randomly extract two groups of 30 data points from Data 3 and Data 4, respectively, denoted as Data
$^*$3 and Data
$^*$4.
Data
$^*$ 3: 
$173.00, 16.00, 140.00, 49.40, 22.70, 133.00, 146.00, 273.00, 583.00, 112.00, 523.00,$
$277.00, 241.00, 10.42, 42.00, 176.00, 218.00, 1417.00, 417.00, 14.48, 140.00, 1146.00,594.00,$
$297.00, 53.62, 154.00, 10.00, 7.00, 165.00, 45.28$
Data
$^*$ 4: 
$133.00, 1776.00, 633.00, 469.00, 92.00, 58.36, 119.00, 146.00, 281.00, 194.00, 179.00,$
$ 74.48, 112.00, 432.00, 159.00, 195.00, 68.46, 339.00, 47.38, 94.00, 84.00, 140.00, 41.35, 319.00,$
$725.00, 63.47, 78.26, 173.00, 110.00, 31.98$
To exemplify, we examine the reliability 
$R=P(\max\{X_1,X_2\} \gt Y)$. Here, X 1 represents Data 1, X 2 corresponds to Data 2, and Y denotes Data
$^*$ 3. In our initial step, we employed the cor.test(,method=”Kendall”) function in R to execute a Kendall’s tau-based test of dependence between X 1, X 2, and Y. Based on Figure 6, it is evident that there exists a pronounced interdependence among variables 
$X_1,X_2$, and Y, with varying degrees of dependence between them. This observation is further corroborated by the findings presented in Table 9, highlighting the strongest dependency between 
$X_1,X_2$, and Y. As a result, we can employ a hierarchical copula to accurately model the interdependencies among 
$X_1,X_2$, and Y. The dependent structure can be found in Figure 7. Furthermore, to discern the dependence structure between X 1, X 2, and Y, we utilize a goodness-of-fit test for copula. This test is rooted in the multiplier central limit theorems and was introduced by [Reference Kojadinovic, Segers and Yan27]. We present the results of the goodness-of-fit test for various copula models applied to stress variables X 1 and X 2, which constitute the first hierarchical layer, in Table 10. Since 
$\tau(\max\{X_1,X_2\},Y)=-0.1197 \lt 0,$ in order to have comparability, we chose AMH copula in addition to Frank copula mentioned above to fit dependent structure of stress and strength variables 
$\max\{X_1,X_2\}$ and Y, the goodness-of-fit results can be fined in Table 11. The results of the goodness-of-fit test rely on the empirical copula process using the copula distribution function. These results are presented under the null hypothesis 
$H_0: C \in C_{\theta}$ against the alternative hypothesis 
$H_1: C \not\in C_{\theta}$, where the null hypothesis signifies that the data belong to the parametric family. The analysis indicates that the F-F HAC model exhibits a superior fit for X 1, X 2, and Y. In Figure 
$8,$ we give a comparison graph between the actual data (left) and the simulated data (right). Through comparison, it is found that the simulated data are close to the actual data, and the usability of the model is further verified.
Figure 6. Scatter plots of 
$X_1,X_2$, and Y.
Figure 7. The dependent structure of 
$X_1,X_2$, and Y.
Figure 8. Scatter plots of 
$X_1,X_2$, and Y (left) versus simulated data (right).
Table 9. The Kendall tau of Data 1, Data 2, and Data
$^*$ 3.
Table 10. Goodness-of-fit results for (
$X_1,X_2$)， where the significance of the bold value indicates the optimal result.
Table 11. Goodness-of-fit results for (
$\max\{X_1,X_2\}$,Y), where the significance of the bold value indicates the optimal result.
Table 12 showcases the estimated reliability for the scenario 
$R=P(\max\{X_1,X_2\} \gt Y)$, acquired through both empirical and semi-parametric methodologies. In this instance, it’s notable that the empirical and semi-parametric reliabilities do not exhibit close proximity to each other.
Table 12. Reliability results of stress-strength model for Data 1, Data 2, and Data
$^{*}$ 3.
Through further research, we have discovered that when considering Data 1 and Data 2 as strength variables and Data
$^*$ 4 as stress variable, the interdependence structures between Data 1, Data 2, and Data
$^*$ 4 still differ (see Table 13 and Figure 9). However, after conducting goodness-of-fit tests, we found that the interdependence structures between Data 1 and Data 2, as well as between the pair (Data 1, Data 2) and Data
$^*$ 4, are distinct and cannot be characterized using the same family of copulas, that is, combine with Tables 10 and 14, we can use F-J HAC to characterize the interdependence structures between Data 1, Data 2, and Data
$^*$ 4. Similar phenomena have also been observed in the case of Data 2 and Data 3 as strength variables and Data
$^*$ 4 as stress variable (see Table 15 and Figure 10), combine with Tables 16 and 17, we can use C-F HAC to characterize the interdependence structures between Data 1, Data 2, and Data
$^*$ 4. Therefore, in the future, we will study multicomponent stress-strength reliability based on heterogeneous HACs methods.
Figure 9. Scatter plots of 
$X_2,X_3$, and X 4.
Figure 10. Scatter plots of 
$X_1,X_2$, and X 4.
Table 13. The Kendall tau of Data 1, Data 2, and Data
$^*$ 4.
Table 14. Goodness-of-fit results for (Data 1, Data 2) Data
$^*$ 4, where the significance of the bold value indicates the optimal result.
Table 15. The Kendall tau of Data 2, Data
$^*$ 3, and Data
$^*$ 4.
Table 16. Goodness-of-fit results for (Data 2, Data
$^*$ 3), where the significance of the bold value indicates the optimal result.
Table 17. Goodness-of-fit results for (Data 2, Data
$^*$ 3) Data
$^*$ 4, where the significance of the bold value indicates the optimal result.
5. Conclusion
This study underscores the significance of accounting for dependence among strengths and the interdependent between strengths and stress in reliability analysis within multicomponent stress-strength models. The conventional assumption of independence might not hold in various practical scenarios, prompting the development of innovative methodologies.
By introducing an AC-based hierarchical dependence approach, this paper has presented a novel framework for effectively modeling these interdependence. The application of four distinct semi-parametric methods to estimate reliability, accompanied by the determination of dependence parameters, has advanced the understanding of stress-strength relationships. Moreover, the derived asymptotic properties of the estimator provide a foundation for its robustness and applicability.
The validation of the proposed methods through Monte Carlo simulations and real-life data sets underscores their practical utility. This research not only enhances the comprehension of complex interrelationships within multicomponent stress-strength models but also paves the way for more accurate and reliable reliability estimation methodologies in the face of realistic dependence scenarios. The integration of HACs and semi-parametric methods showcases a noteworthy advancement in the field, fostering more informed decision-making in diverse engineering and scientific applications.
Looking ahead, this paper primarily considers partially HACs. In future work, we will focus on the application of fully hierarchical HACs in multicomponent stress-strength models. Additionally, while this study emphasizes the reliability of stress-strength models with parallel structures, we will explore more general systems, specifically coherent systems, to further enrich our understanding and methodologies in this area.
Acknowledgments
This work is supported by National Natural Science Foundation of China (Nos. 12361060 and 11861058) and Science and Technology Program of Gansu Province (No. 24JRRA144) and College Teachers Innovation Foundation Project of Gansu Provincial Education Department (No. 2024A-002). The authors are very grateful to editor and two anonymous reviewers for their careful reading of the manuscript and their helpful comments.
Competing interest
The authors declare no competing interest.
Open access
