Multicomponent Stress–Strength Reliability with Extreme Value Distribution Margins: Its Theory and Application to Hydrological Data

Abstract

This paper focuses on the estimation of multicomponent stress–strength models, an important concept in reliability analyses used to determine the probability that a system will function successfully under varying stress conditions. Understanding and accurately estimating these probabilities is essential in fields such as engineering and risk management, where the reliability of components under extreme conditions can have significant consequences. This is the case in applications where one seeks to model extreme hydrological events. Specifically, this study examines cases where the random variables X (representing strength) and Y (representing stress) follow extreme value distributions. New analytical expressions are derived for multicomponent stress–strength reliability (MSSR) when different classes of extreme distributions are considered, using the extreme value $\mathbb{H}$-function. These results are applied to three l-max stable laws and six p-max stable laws, providing a robust theoretical framework for multicomponent stress–strength analyses under extreme conditions. To demonstrate the practical relevance of the proposed models, a real dataset is analyzed, focusing on the monthly water capacity of the Shasta Reservoir in California (USA) during August and December from 1980 to 2015. This application showcases the effectiveness of the derived expressions in modeling real-world data.

1. Introduction

Stress–strength reliability (SSR) consists of the study of the probability of the failure of a system or component based on a comparison of the applied stress to the strength of the system. That is, SSR is defined as $P(X<Y)$, where Y is an observed random stress and X is an observed random strength. There are several applications of this theory, such as in engineering [1] and asset selection [2,3], among others. See Kotz et al. [4] for more details. The SSR theory can be extended to a multicomponent system that consists of k strength components under a common stress, where the system functions when at least s $(1\le s\le k)$ components simultaneously survive. The estimation of reliability in a multicomponent stress–strength model has been studied by several authors: exponential distributions were explored by Bhattacharyya and Johnson [5]; Wang et al. [6] studied the estimation of stress–strength reliability for multicomponent systems with the Rayleigh distribution; a Bayesian estimation of reliability for a multicomponent stress–strength model based on Topp–Leone distribution was studied by Akram and Yousif [7]; and Lio et al. [8] developed an inferential methodology for multicomponent stress–strength reliability when the strength and stress variables follow Burr XII distributions under Type-II censoring on strength data.

Multicomponent stress–strength analysis holds significant importance in hydrology because it allows for the evaluation of the reliability and resilience of water resource systems under extreme environmental conditions. This type of analysis is particularly useful in assessing how well a system, such as a reservoir or flood control infrastructure, can withstand the stresses caused by natural events like heavy rainfall, drought, or flooding. By comparing the system’s strength (e.g., reservoir capacity) against the stress imposed by hydrological events, researchers and engineers can predict the likelihood of system failure. This predictive capability is essential for designing robust water management strategies, ensuring that systems can continue to function effectively under varying and extreme conditions.

Although spatial information is absent, multicomponent stress–strength assessments still offer hydrologists a simple yet effective way to understand water system stability across broad scenarios. This approach allows hydrologists to simulate different conditions and quantify failure risks, supporting conservation strategies, drought response planning, and water management decisions. By focusing on the interactions among critical factors instead of spatial variations, multicomponent stress–strength studies provide a high-level perspective on system robustness, helping hydrologists identify potential system weaknesses and develop plans for sustainable resource management.

In the context of the present paper, let $X_1,X_2,\cdots ,X_k$ denote the strength variables of k components in a system and suppose they follow a common cumulative distribution function (CDF) $F_X(·)$. Each component is subject to an independent random stress Y, which follows the probability density function (PDF) $f_Y(·)$. The multicomponent stress–strength reliability (MSSR) model is given by

$$\begin{array}{ccc}\hfill R_{s,k}& =& P\left(\mathrm{at} \mathrm{least} \mathrm{s} \mathrm{of}(X_1,\cdots ,X_k)\mathrm{exceed}Y\right)\hfill \\ & =& \sum _{j=s}^k\left(\begin{array}{c}k\\ j\end{array}\right){\int }_{-\infty }^{\infty }{\left(1-F_X\left(u\right)\right)}^j{\left(F_X\left(u\right)\right)}^{k-j}{f}_Y\left(u\right)du.\hfill \end{array}$$

(1)

In the context of extreme events, it is of interest to model the random variables in Equation (1) using extreme value theory (EVT). In particular, let one consider the Fisher–Tippett–Gnedenko theorem, which describes the limiting behavior of the maximum of a sequence of independent and identically distributed (i.i.d.) random variables. Let ${\left\{X_i\right\}}_{i=1}^M$ be a M-sized sample of i.i.d. random variables with CDF $F(·)$. If there exist normalizing constants ${\{a_n>0\}}_{n=1}^M$ and ${\{b_n\in \mathbb{R}\}}_{n=1}^M$ such that the scaled maximum converges in its distribution to a non-trivial limit, then this limit must be one of three types of extreme value distributions, namely: Gumbel, Fréchet or Weibull. Mathematically,

$$\underset{n\to \infty }{lim}P\left({\displaystyle \frac{max\{X_1,\dots ,X_n\}-b_n}{a_n}}\le x\right)=G\left(x\right)$$

(2)

or, equivalently,

$$\underset{n\to \infty }{lim}{\left(F(a_{n}x+b_n)\right)}^n=G\left(x\right).$$

(3)

In other words, if the limit above converges, then, up to a linear change of coordinates, $G\left(x\right)$ will assume either the form of the CDF of the Gumbel, the Fréchet, or the Weibull distribution and can be referred to as an extreme value distribution or stable l-max law (see, for example, Resnick [9] and Mohan and Ravi [10]). The conditions under which the distribution function $G(·)$ for (2) holds true are well documented and referenced in Resnick [9], Galambos [11], and Haan and Ferreira [12]. Applications of EVT are found in finance, natural catastrophes, and equipment failures, among others (cf. Embrechts et al. [13]).

A more general normalization of the partial maximum $Ma{x}_n=\max {\left\{X_i\right\}}_{i=1}^M$ can also be considered. It was shown in Pantcheva [14] that under power normalization there are six classes of CDFs $H(·)$ such that

$$\underset{n\to \infty }{lim}{F}^n\left({\alpha }_n{\left|x\right|}^{{\beta }_n}sign\left(x\right)\right)=H\left(x\right),$$

(4)

where ${\alpha }_n,{\beta }_n>0$ are sequences of positive real numbers. These distributions are called extreme value distributions under power normalization or p-max stable laws. The expansion of definitions and diverse properties from l-max stable distributions to p-max stable distributions happens naturally. For a complete characterization of these distributions see Mohan and Ravi [10].

The goal of this paper is to study estimations of (1) when X and Y follow extreme value distributions. New analytical expressions are obtained for the MSSR of several classes of extreme distributions. In particular, for l-max and p-max cases, all the expressions for MSSR are novel and were developed in the present paper. Some of the partial results were obtained in other papers (e.g., the GEV case [3]) or have their proof related to previously published results (e.g., p-max [15] and l-max cases [16]). In general, except for the GEV case, where MSSR was covered, in [3], previous works have covered single-component SSR situations, while the present paper presents novel analytical results for multicomponent SSR scenarios. Therefore, the present paper uncovers unexplored results.

Its application to a real dataset is carried out to show the performance of such expressions in data modeling, illustrating how these new results are useful to hydrologists. The data represent the monthly water capacity of the Shasta Reservoir in California (USA) during August and December from 1980 to 2015 and were previously modeled in [6,8,17].

This paper is organized as follows: Section 2 presents general expressions for (1) in terms of a special function called the extreme-value $\mathbb{H}$-function. The results are presented for three l-max stable laws and six p-max stable laws. A summarizing result for the l-max stable laws using the GEV distribution is also presented in Section 2. Section 3, on the other hand, deals with real dataset applications of the analytical expressions found. The last section deals with conclusions.

2. Multicomponent System Reliability

Let $X_1,\cdots ,X_n$ be i.i.d. RVs with CDF $F_X(·)$ and let Y be an independent RV with PDF $f_Y(·)$. Set $M_n=max\{X_1,\cdots ,X_n\}$. Then, $P(M_n\le u)=F_X^n\left(u\right)$ and we find that

$$P(X_1<Y,\cdots ,X_n<Y)=P(M_n\le Y)={\int }_{-\infty }^{\infty }{F}_X^n\left(u\right)f_Y\left(u\right)du=:I_n.$$

(5)

Closed expressions for (5) are presented in the Results subsections for each family of extreme distributions.

In a broader context, the multicomponent stress–strength reliability (MSSR) is given by

$$\begin{array}{ccc}\hfill R_{s,k}& =& P\left(\mathrm{at} \mathrm{least} \mathrm{s} \mathrm{of}(X_1,\cdots ,X_k)\mathrm{exceed}Y\right)\hfill \\ & =& \sum _{j=s}^k\left(\begin{array}{c}k\\ j\end{array}\right){\int }_{-\infty }^{\infty }{\left(1-F_X\left(u\right)\right)}^j{\left(F_X\left(u\right)\right)}^{k-j}{f}_Y\left(u\right)du.\hfill \end{array}$$

Using a binomial expansion, we obtain

$$R_{s,k}=\sum _{j=s}^k\sum _{r=0}^j\left(\begin{array}{c}k\\ j\end{array}\right)\left(\begin{array}{c}j\\ r\end{array}\right){(-1)}^{j-r}{\int }_{-\infty }^{\infty }{\left(F_X\left(u\right)\right)}^{k-r}{f}_Y\left(u\right)du.$$

(6)

Note that the integral terms in (6) are particular cases of (5) when $n=k-r$. Therefore,

$$R_{s,k}=\sum _{j=s}^k\sum _{r=0}^j\left(\begin{array}{c}k\\ j\end{array}\right)\left(\begin{array}{c}j\\ r\end{array}\right){(-1)}^{j-r}{I}_{k-r}.$$

(7)

We are interested in obtaining $R_{s,k}$ for three l-max stable laws and six p-max stable laws. From (7), we know that it is enough to have a closed expression to $I_n$ defined in (5).

Recently, Rathie et al. [18] introduced the extreme-value $\mathbb{H}$-function as

$$\mathbb{H}(a_1,a_2,a_3,a_4,a_5,a_6):={\int }_0^{\infty }{y}^{a_6}exp\{-a_{1}y-{(a_2{y}^{a_3}+a_4)}^{a_5}\}dy,$$

(8)

where $\Re \left(a_1\right),\Re \left(a_2\right),\Re \left(a_4\right)\in {\mathbb{R}}_{\ge 0},a_3,a_5\in \mathbb{C}$, both $\Re \left(a_1\right)$ and $\Re \left(a_2\right)$ cannot be equal to zero at the same time, $\Re \left(a_6\right)>-1$ when $a_1\ne 0$ or $a_1=0$ and $\mathrm{sign}(\Re \left(a_3\right))=\mathrm{sign}(\Re \left(a_5\right))$, and $\Re \left(a_6\right)<-1$ when $a_1=0$ and $\mathrm{sign}(\Re \left(a_3\right))\ne \mathrm{sign}(\Re \left(a_5\right))$. In this paper, $\mathbb{R}$, $\mathbb{C}$, and ℜ denote real numbers, complex numbers, and the real part of a complex number, respectively.

The extreme-value $\mathbb{H}$-function is an important tool for obtaining closed expressions of $I_n$ for each extreme value distribution.

2.1. An Alternative Expression for MSSR

The MSSR (7) can be rewritten as

$$R_{s,k}=1-\sum _{\omega =k-s+1}^k\left(\begin{array}{c}k\hfill \\ w\hfill \end{array}\right){\int }_{-\infty }^{+\infty }{\left(F_X\left(u\right)\right)}^w{\left(1-F_X\left(u\right)\right)}^{k-w}{f}_Y\left(u\right)du.$$

(9)

Indeed,

$$\begin{array}{cc}\hfill R_{s,k}& =P\left(\mathrm{at}\mathrm{least}\mathrm{s}\mathrm{of}\left(X_1,\dots ,X_k\right)\mathrm{exceed}Y\right)\hfill \\ & =1-P\left(\mathrm{at}\mathrm{least}\mathrm{s}\mathrm{of}\left(X_1,\dots ,X_k\right)\mathrm{do}\mathrm{not}\mathrm{exceed}Y\right)\hfill \\ & ={\int }_{-\infty }^{+\infty }\left[1-\sum _{\omega =k-s+1}^k\left(\begin{array}{c}k\hfill \\ w\hfill \end{array}\right){\left(F_X\left(u\right)\right)}^w{\left(1-F_X\left(u\right)\right)}^{k-\omega }\right]f_Y\left(u\right)du,\hfill \end{array}$$

since ${\int }_{-\infty }^{+\infty }{f}_Y\left(u\right)du=1$.

Note that

$$1={\left(F_X\left(u\right)+\left(1-F_X\left(u\right)\right)\right)}^k=\sum _{n=0}^k\left(\begin{array}{c}k\hfill \\ n\hfill \end{array}\right){\left(F_X\left(u\right)\right)}^n{\left(1-F_X\left(u\right)\right)}^{k-n}.$$

Thus,

$$1-\sum _{w=k-s+1}^k\left(\begin{array}{c}k\hfill \\ w\hfill \end{array}\right){\left(F_X\left(u\right)\right)}^w{\left(1-F_X\left(u\right)\right)}^{k-w}=\sum _{w=0}^k\left(\begin{array}{c}k\hfill \\ w\hfill \end{array}\right){\left(F_X\left(u\right)\right)}^w{\left(1-F_X\left(u\right)\right)}^{k-w}.$$

Set $k-w=j$. Noting that $\left(\begin{array}{c}k\\ k-j\end{array}\right)=\left(\begin{array}{c}k\hfill \\ w\hfill \end{array}\right)$, we obtain

$$1-\sum _{w=k-s+1}^k\left(\begin{array}{c}k\hfill \\ w\hfill \end{array}\right){\left(F_X\left(u\right)\right)}^w{\left(1-F_X\left(u\right)\right)}^{k-w}=\sum _{j=s}^k\left(\begin{array}{c}k\hfill \\ s\hfill \end{array}\right){\left(F_X\left(u\right)\right)}^{k-j}{\left(1-F_X\left(u\right)\right)}^j,$$

In the next subsections, we describe the cumulative distribution function (CDF) and the probability density function (PDF) of l- and p-max stable laws and we obtain the corresponding expression of $I_n$ in each case.

2.2. l-Max Stable Laws

2.2.1. Fréchet

We say that a random variable (RV) X has a Fréchet distribution, and denote this as $X\sim$ Fréchet$(\mu ,\sigma ,\alpha )$, with a location parameter $\mu \in \mathbb{R}$, scale parameter $\sigma \in {\mathbb{R}}_{+}$, and shape parameter $\alpha \in {\mathbb{R}}_{+}$, if X has a CDF and PDF, respectively, given by

$$F(x;\mu ,\sigma ,\alpha )=\left\{\begin{array}{cc}0,& x<\mu ,\\ exp\left\{-{\left({\displaystyle \frac{x-\mu }{\sigma }}\right)}^{-\alpha }\right\},& x\ge \mu ,\end{array}\right.$$

and

$$f(x;\mu ,\sigma ,\alpha )={\displaystyle \frac{\alpha }{\sigma }}{\left({\displaystyle \frac{x-\mu }{\sigma }}\right)}^{-\alpha -1}F(x;\mu ,\sigma ,\alpha ){\mathbb{1}}_{(\mu ,\infty )},$$

(10)

where ${\mathbb{1}}_A$ denotes the indicator function on set A. Figure 1 shows the behavior of Fréchet PDFs when parameters vary.

Let $X_1,\cdots ,X_n$ be i.i.d. RVs with the distribution Fréchet$({\mu }_x,{\sigma }_x,{\gamma }_x)$ and Y be an independent RV with the distribution Fréchet$({\mu }_y,{\sigma }_y,{\gamma }_y)$. If ${\mu }_y\ge {\mu }_x$, then

$$I_n=\mathbb{H}\left(1,n^{-1/{\alpha }_x}{\displaystyle \frac{{\sigma }_y}{{\sigma }_x}},-{\displaystyle \frac{1}{{\alpha }_y}},n^{-1/{\alpha }_x}{\displaystyle \frac{{\mu }_y-{\mu }_x}{{\sigma }_x}},-{\alpha }_x,0\right),$$

(11)

Indeed,

$$\begin{array}{cc}\hfill I_n& =P\left(X_1\le Y,\dots ,X_n\le Y\right)\hfill \\ & =P\left(M_n\le Y\right)\hfill \\ & ={\int }_{-\infty }^{+\infty }{F}_X^n(u;{\mu }_x,{\sigma }_x,{\gamma }_x)f_Y(u;{\mu }_y,{\sigma }_y,{\gamma }_y)du.\hfill \end{array}$$

Set $M=max\left\{{\mu }_x,{\mu }_y\right\}$. Then,

$$I_n={\int }_M^{+\infty }exp\left\{-n{\left({\displaystyle \frac{u-{\mu }_x}{{\sigma }_x}}\right)}^{-{\alpha }_x}\right\}{\displaystyle \frac{{\alpha }_y}{{\sigma }_y}}{\left({\displaystyle \frac{u-{\mu }_y}{{\sigma }_y}}\right)}^{-{\alpha }_y-1}exp\left\{-{\left({\displaystyle \frac{u-{\mu }_y}{{\sigma }_y}}\right)}^{-{\alpha }_y}\right\}du.$$

Substituting $v={\left({\displaystyle \frac{u-{\mu }_y}{{\sigma }_y}}\right)}^{-{\alpha }_y}$, if $M={\mu }_y$, we obtain

$$I_n={\int }_0^{\infty }exp\left\{-v-n{\left({\displaystyle \frac{{\sigma }_y}{{\sigma }_x}}{v}^{-1/{\alpha }_y}+{\displaystyle \frac{{\mu }_y-{\mu }_x}{{\sigma }_x}}\right)}^{-{\alpha }_x}\right\}dv.$$

(12)

Hence, (11) follows from (8) and (12).

2.2.2. Reversed Weibull

We say that an RV X has a Reversed Weibull distribution, and denote this by $X\sim$ RWeibull$(\mu ,\sigma ,\alpha )$, with a location parameter $\mu \in \mathbb{R}$, scale parameter $\sigma \in {\mathbb{R}}_{+}$, and shape parameter $\alpha \in {\mathbb{R}}_{+}$, if X has a CDF and PDF, respectively, given by

$$F(x;\mu ,\sigma ,\alpha )=\left\{\begin{array}{cc}exp\left\{-{\left({\displaystyle \frac{\mu -x}{\sigma }}\right)}^{\alpha }\right\},& x<\mu ,\\ 1,& x\ge \mu ,\end{array}\right.$$

and

$$f(x;\mu ,\sigma ,\alpha )={\displaystyle \frac{\alpha }{\sigma }}{\left({\displaystyle \frac{\mu -x}{\sigma }}\right)}^{\alpha -1}F(x;\mu ,\sigma ,\alpha )1_{(-\infty ,\mu )}\left(x\right).$$

(13)

Figure 2 shows the behavior of Reversed Weibull PDFs when parameters vary.

Let $X_1,\cdots ,X_n$ be i.i.d. RVs with the distribution RWeibull$({\mu }_x,{\sigma }_x,{\gamma }_x)$ and Y be an independent RV with the distribution RWeibull$({\mu }_y,{\sigma }_y,{\gamma }_y)$. If ${\mu }_x\ge {\mu }_y$, then

$$I_n=\mathbb{H}\left(1,n^{1/{\alpha }_x}{\displaystyle \frac{{\sigma }_y}{{\sigma }_x}},{\displaystyle \frac{1}{{\alpha }_y}},n^{1/{\alpha }_x}{\displaystyle \frac{{\mu }_x-{\mu }_y}{{\sigma }_x}},{\alpha }_x,0\right).$$

(14)

The proof of Equation (14) follows a similar rationale as that in the case of Equation (11) and will be skipped for simplicity.

2.2.3. Gumbel

Finally, we say that a random variable X has a Gumbel distribution, and denote this by $X\sim$ Gumbel$(\mu ,\sigma )$, with location parameter $\mu \in \mathbb{R}$ and scale parameter $\sigma \in {\mathbb{R}}_{+}$, if X has a CDF and PDF, respectively, given by

$$F(x;\mu ,\sigma ,\alpha )=exp\left\{-exp\left\{-{\displaystyle \frac{x-\mu }{\sigma }}\right\}\right\},$$

and

$$f(x;\mu ,\sigma ,\alpha )={\displaystyle \frac{1}{\sigma }}exp\left\{-{\displaystyle \frac{x-\mu }{\sigma }}\right\}exp\left\{-exp\left\{-{\displaystyle \frac{x-\mu }{\sigma }}\right\}\right\}.$$

Figure 3 presents Gumbel PDFs generated using different parameter sets.

Let $X_1,\cdots ,X_n$ be i.i.d. RVs with the distribution Gumbel$({\mu }_x,{\sigma }_x)$ and Y be an independent RV with the distribution Gumbel$({\mu }_y,{\sigma }_y)$. Then,

$$I_n=exp\left\{{\displaystyle \frac{{\mu }_y}{{\sigma }_y}}\right\}\mathbb{H}\left(exp\left\{{\displaystyle \frac{{\mu }_y}{{\sigma }_y}}\right\},nexp\left\{{\displaystyle \frac{{\mu }_x}{{\sigma }_x}}\right\},{\displaystyle \frac{{\sigma }_y}{{\sigma }_x}},0,1,0\right).$$

(15)

The proof of Equation (15) follows a similar rationale as that in the case of Equation (11) and will be skipped for simplicity.

2.2.4. GEV Distribution

For statistical applications, the l-max stable distributions can be summarized in a single distribution function called the generalized extreme value distribution (GEV). With a reading of Embrechts et al. [13], it is possible to conduct a detailed study of the GEV distribution.

We denote, using $X\sim GEV(\mu ,\sigma ,\gamma )$, $\mu ,\gamma \in \mathbb{R}$, and $\sigma >0$, an RV with a CDF given by

$$F(x;\mu ,\sigma ,\gamma )=exp\left\{-{\left[1+{\displaystyle \frac{\gamma }{\sigma }}\left(x-\mu \right)\right]}^{-{\displaystyle \frac{1}{\gamma }}}\right\},1+{\displaystyle \frac{\gamma }{\sigma }}(x-\mu )>0,$$

(16)

where $\mu$ is the location parameter, $\sigma$ is the scale parameter, and $\gamma$ is the shape parameter. The corresponding PDF is given by

$$f(x;\mu ,\sigma ,\gamma )=F(x,\gamma ,\mu ,\sigma ){\displaystyle \frac{1}{\sigma }}{\left[1+{\displaystyle \frac{\gamma }{\sigma }}\left(x-\mu \right)\right]}^{-{\displaystyle \frac{1}{\gamma }}-1},1+{\displaystyle \frac{\gamma }{\sigma }}(x-\mu )>0.$$

(17)

Figure 4 shows the behavior of (17) with different parameters. Note that the location parameter shifts the curve, the scale controls dispersion, and the density changes according to the sign of the shape.

Let $X_1,\cdots ,X_n$ be i.i.d. RVs with the distribution $GEV({\mu }_x,{\sigma }_x,{\gamma }_x)$ and Y be an independent RV with the distribution $GEV({\mu }_y,{\sigma }_y,{\gamma }_y)$. Then,

• When ${\gamma }_j>0$ and $j=x,y$,

  $$I_n=\mathbb{H}\left(1,{\displaystyle \frac{{\gamma }_x{\sigma }_y}{{\sigma }_x{\gamma }_y}}{n}^{-{\gamma }_x},-{\gamma }_y,\left[1+{\displaystyle \frac{{\gamma }_x}{{\sigma }_x}}\left({\mu }_y-{\mu }_x-{\displaystyle \frac{{\sigma }_y}{{\gamma }_y}}\right)\right]n^{-{\gamma }_x},-{\displaystyle \frac{1}{{\gamma }_x}},0\right),$$

  (18)

  provided that ${\mu }_y-{\displaystyle \frac{{\sigma }_y}{{\gamma }_y}}\ge {\mu }_x-{\displaystyle \frac{{\sigma }_x}{{\gamma }_x}}$. When ${\mu }_y-{\displaystyle \frac{{\sigma }_y}{{\gamma }_y}}\le {\mu }_x-{\displaystyle \frac{{\sigma }_x}{{\gamma }_x}}$,

  $$I_n=1-\mathbb{H}\left(1,{\displaystyle \frac{{\gamma }_y{\sigma }_x}{{\sigma }_y{\gamma }_x}}{n}^{-{\gamma }_y},-{\gamma }_x,\left[1+{\displaystyle \frac{{\gamma }_y}{{\sigma }_y}}\left({\mu }_x-{\mu }_y-{\displaystyle \frac{{\sigma }_x}{{\gamma }_x}}\right)\right]n^{-{\gamma }_y},-{\displaystyle \frac{1}{{\gamma }_y}},0\right).$$

  (19)
• When ${\gamma }_j<0$ and $j=x,y$,

  $$I_n=\mathbb{H}\left(1,{\displaystyle \frac{{\gamma }_x{\sigma }_y}{{\sigma }_x{\gamma }_y}}{n}^{-{\gamma }_x},-{\gamma }_y,\left[1+{\displaystyle \frac{{\gamma }_x}{{\sigma }_x}}\left({\mu }_y-{\mu }_x-{\displaystyle \frac{{\sigma }_y}{{\gamma }_y}}\right)\right]n^{-{\gamma }_x},-{\displaystyle \frac{1}{{\gamma }_x}},0\right),$$

  (20)

  provided that ${\mu }_y-{\displaystyle \frac{{\sigma }_y}{{\gamma }_y}}\le {\mu }_x-{\displaystyle \frac{{\sigma }_x}{{\gamma }_x}}$. When ${\mu }_y-{\displaystyle \frac{{\sigma }_y}{{\gamma }_y}}\ge {\mu }_x-{\displaystyle \frac{{\sigma }_x}{{\gamma }_x}}$

  $$I_n=1-\mathbb{H}\left(1,{\displaystyle \frac{{\gamma }_y{\sigma }_x}{{\sigma }_y{\gamma }_x}}{n}^{-{\gamma }_y},-{\gamma }_x,\left[1+{\displaystyle \frac{{\gamma }_y}{{\sigma }_y}}\left({\mu }_x-{\mu }_y-{\displaystyle \frac{{\sigma }_x}{{\gamma }_x}}\right)\right]n^{-{\gamma }_y},-{\displaystyle \frac{1}{{\gamma }_y}},0\right).$$

  (21)
• When ${\gamma }_y={\gamma }_x=0$,

  $$I_n=exp\left\{{\displaystyle \frac{{\mu }_y}{{\sigma }_y}}\right\}\mathbb{H}\left(exp\left\{{\displaystyle \frac{{\mu }_y}{{\sigma }_y}}\right\},exp\left\{{\displaystyle \frac{{\mu }_x}{{\sigma }_x}}\right\}n,{\displaystyle \frac{{\sigma }_y}{{\sigma }_x}},0,1,0\right).$$

  (22)

  Proofs for the expressions (18) to (22) can be found in [3].

2.3. p-Max Stable Laws

For all the p-max stable laws, proofs of the expression presented follow from previously published results in [15] by considering standard parametric changes.

2.3.1. Case 1: $H_1$

Let $X_1,\cdots ,X_n$ be i.i.d. RVs with the distribution $H_1({\alpha }_x,{\beta }_x,{\gamma }_x)$ and Y be an independent RV with the distribution $H_1({\alpha }_y,{\beta }_y,{\gamma }_y)$, ${\alpha }_i,{\beta }_i,{\gamma }_i\in (0,\infty )$, $i\in \{x,y\}$, where $H_1(·)$ and $h_1(·)$ are the CDF and PDF, respectively, and can be defined as follows:

$$H_1(x;\alpha ,\beta ,\gamma )=\left\{\begin{array}{ccc}0,\hfill & \mathrm{if}& \gamma x^{\beta }<1,\hfill \\ exp\{-{(log\left(\gamma x^{\beta }\right))}^{-\alpha }\},\hfill & \mathrm{if}& \gamma x^{\beta }\ge 1,\hfill \end{array}\right.$$

$$h_1(x;\alpha ,\beta ,\gamma )=exp\left\{-{[log\gamma x^{\beta }]}^{-\alpha }\right\}{\displaystyle \frac{\alpha \beta }{x}}{[log\left(\gamma x^{\beta }\right)]}^{-\alpha -1}{\mathbb{1}}_{({\gamma }^{-1/\beta },\infty )}\left(x\right),$$

then

$$\begin{array}{cc}& I_n=\mathbb{H}\left(1,n^{-1/{\alpha }_x}{\displaystyle \frac{{\beta }_x}{{\beta }_y}},-{\displaystyle \frac{1}{{\alpha }_y}},n^{-1/{\alpha }_x}\left(log{\gamma }_x-{\displaystyle \frac{{\beta }_x}{{\beta }_y}}log{\gamma }_y\right),-{\alpha }_x,0\right),\hfill \end{array}$$

provided that ${\gamma }_y^{-1/{\beta }_y}\ge {\gamma }_x^{-1/{\beta }_x}$. Figure 5 presents the $H_1$ PDFs generated using different parameter sets.

2.3.2. Case 2: $H_2$

Let $X_1,\cdots ,X_n$ be i.i.d. RVs with the distribution $H_2({\alpha }_x,{\beta }_x,{\gamma }_x)$ and Y be an independent RV with the distribution $H_2({\alpha }_y,{\beta }_y,{\gamma }_y)$, ${\alpha }_i,{\beta }_i,{\gamma }_i\in (0,\infty )$, $i\in \{x,y\}$, where $H_2(·)$ and $h_2(·)$ are the CDF and PDF, respectively, and can be defined as follows:

$$H_2(x;\alpha )=\left\{\begin{array}{ccc}0,\hfill & \mathrm{if}& \gamma x^{\beta }<0,\hfill \\ exp\{-{(-log\left(\gamma x^{\beta }\right))}^{\alpha }\},\hfill & \mathrm{if}& 0\le \gamma x^{\beta }<1,\hfill \\ 1,\hfill & \mathrm{if}& \gamma x^{\beta }\ge 1,\hfill \end{array}\right.$$

$$h_2(x;\alpha ,\beta ,\gamma )=exp\left\{-{\left[-log\left(\gamma x^{\beta }\right)\right]}^{\alpha }\right\}{\left[-log\left(\gamma x^{\beta }\right)\right]}^{\alpha -1}{\displaystyle \frac{\beta \alpha }{x}}{\mathbb{1}}_{(0,{\gamma }^{-1/\beta })}\left(x\right),$$

then

$$I_n=\mathrm{H}\left(1,n^{1/{\alpha }_x}{\displaystyle \frac{{\beta }_x}{{\beta }_y}},{\displaystyle \frac{1}{{\alpha }_y}},n^{1/{\alpha }_x}\left(-log{\gamma }_x+{\displaystyle \frac{{\beta }_x}{{\beta }_y}}log{\gamma }_y\right),{\alpha }_x,0\right),$$

provided that ${\gamma }_y^{-1/{\beta }_y}\le {\gamma }_x^{-1/{\beta }_x}$. Figure 6 presents the $H_2$ PDFs generated using different parameter sets.

2.3.3. Case 3: $H_3$

Let $X_1,\cdots ,X_n$ be i.i.d. RVs with the distribution $H_3({\alpha }_x,{\beta }_x,{\gamma }_x)$ and Y be an independent RV with the distribution $H_3({\alpha }_y,{\beta }_y,{\gamma }_y)$, ${\alpha }_i,{\beta }_i,{\gamma }_i\in (0,\infty )$, $i\in \{x,y\}$, where $H_3(·)$ and $h_3(·)$ are the CDF and PDF, respectively, and can be defined as follows:

$$H_3(x;\alpha )=\left\{\begin{array}{ccc}0,\hfill & \mathrm{if}& -\gamma {(-x)}^{\beta }<-1,\hfill \\ exp\{-{(-log\left(\gamma {(-x)}^{\beta }\right))}^{-\alpha }\},\hfill & \mathrm{if}& -1\le -\gamma {(-x)}^{\beta }<0,\hfill \\ 1,\hfill & \mathrm{if}& -\gamma {(-x)}^{\beta }\ge 0,\hfill \end{array}\right.$$

$$h_3(x;\alpha ,\beta ,\gamma )=exp\left\{-{\left[-log\left(\gamma {(-x)}^{\beta }\right)\right]}^{-\alpha }\right\}{\displaystyle \frac{\alpha \beta }{-x}}{\left[-log\left(\gamma {(-x)}^{\beta }\right)\right]}^{-\alpha -1}{\mathbb{1}}_{(-1,0)}\left(x\right),$$

then

$$I_n=H\left(1,n^{-{\displaystyle \frac{1}{{\alpha }_x}}}{\displaystyle \frac{{\beta }_x}{{\beta }_y}},-{\displaystyle \frac{1}{{\alpha }_y}},n^{-{\displaystyle \frac{1}{{\alpha }_x}}}\left(-log{\gamma }_x+{\displaystyle \frac{{\beta }_x}{{\beta }_y}}log{\gamma }_y\right),-{\alpha }_x,0\right),$$

provided that ${\gamma }_y^{-1/{\beta }_y}\le {\gamma }_x^{-1/{\beta }_x}$. Figure 7 presents the $H_3$ PDFs generated using different parameter sets.

2.3.4. Case 4: $H_4$

Let $X_1,\cdots ,X_n$ be i.i.d. RVs with the distribution $H_4({\alpha }_x,{\beta }_x,{\gamma }_x)$ and Y be an independent RV with the distribution $H_4({\alpha }_y,{\beta }_y,{\gamma }_y)$, ${\alpha }_i,{\beta }_i,{\gamma }_i\in (0,\infty )$, $i\in \{x,y\}$, where $H_4(·)$ and $h_4(·)$ are the CDF and PDF, respectively, and can be defined as follows:

$$H_4(x;\alpha )=\left\{\begin{array}{ccc}exp\{-{(log\left(\gamma {(-x)}^{\beta }\right))}^{\alpha }\},\hfill & \mathrm{if}& -\gamma {(-x)}^{\beta }<-1,\hfill \\ 1,\hfill & \mathrm{if}& -\gamma {(-x)}^{\beta }\ge -1,\hfill \end{array}\right.$$

$$h_4(x;\alpha ,\beta ,\gamma )=exp\{-{\left[log\left(\gamma {(-x)}^{\beta }\right)\right]}^{\alpha }\}{\left[log\left(\gamma {(-x)}^{\beta }\right)\right]}^{\alpha -1}{\displaystyle \frac{\alpha \beta }{-x}}{\mathbb{1}}_{(-\infty ,-1)}\left(x\right),$$

then

$$I_n=\mathbb{H}\left(1,n^{{\displaystyle \frac{1}{{\alpha }_x}}}{\displaystyle \frac{{\beta }_x}{{\beta }_y}},{\displaystyle \frac{1}{{\alpha }_y}},n^{{\displaystyle \frac{1}{{\alpha }_x}}}\left(log{\gamma }_x-{\displaystyle \frac{{\beta }_x}{{\beta }_y}}log{\gamma }_y\right),{\alpha }_x,0\right),$$

provided that ${\gamma }_y^{1/{\beta }_y}\le {\gamma }_x^{1/{\beta }_x}$. Figure 8 presents the $H_4$ PDFs generated using different parameter sets.

2.3.5. Case 5: $H_5$

Let $X_1,\cdots ,X_n$ be i.i.d. RVs with the distribution $H_5({\beta }_x,{\gamma }_x)$ and Y be an independent RV with the distribution $H_5({\beta }_y,{\gamma }_y)$, ${\beta }_i,{\gamma }_i\in (0,\infty )$, $i\in \{x,y\}$, where $H_5(·)$ and $h_5(·)$ are the CDF and PDF, respectively, and can be defined as follows:

$$H_5\left(x\right)=\left\{\begin{array}{ccc}0,\hfill & \mathrm{if}& \gamma x^{\beta }<0,\hfill \\ exp\{-{\left(\gamma x^{\beta }\right)}^{-1}\}\hfill & \mathrm{if}& \gamma x^{\beta }\ge 0,\hfill \end{array}\right.$$

$$h_5(x;\beta ,\gamma )=exp\left\{-{\left(\gamma x^{\beta }\right)}^{-1}\right\}{\displaystyle \frac{\beta }{\gamma }}{x}^{-\beta -1}{\mathbb{1}}_{[0,\infty )}\left(x\right),$$

then

$$I_n={\displaystyle \frac{1}{{\gamma }_y}}\mathbb{H}\left({\gamma }_y^{-1},n{\gamma }_x^{-1},{\displaystyle \frac{{\beta }_x}{{\beta }_y}},0,1,0\right).$$

Figure 9 presents the $H_5$ PDFs generated using different parameter sets.

2.3.6. Case 6: $H_6$

Let $X_1,\cdots ,X_n$ be i.i.d. RVs with the distribution $H_6({\beta }_x,{\gamma }_x)$ and Y be an independent RV with the distribution $H_6({\beta }_y,{\gamma }_y)$, ${\alpha }_i,{\beta }_i,{\gamma }_i\in (0,\infty )$, $i\in \{x,y\}$, where $H_6(·)$ and $h_6(·)$ are defined as follows:

$$H_6\left(x\right)=\left\{\begin{array}{ccc}exp\{-\gamma {(-x)}^{\beta }\},\hfill & \mathrm{if}& -\gamma {(-x)}^{\beta }<0,\hfill \\ 1\hfill & \mathrm{if}& -\gamma {(-x)}^{\beta }\ge 0,\hfill \end{array}\right.$$

$$h_6(x;\beta ,\gamma )=exp\{-\gamma {(-x)}^{\beta }\}\gamma \beta {(-x)}^{\beta -1}{\mathbb{1}}_{(-\infty ,0)}\left(x\right),$$

then

$$I_n={\gamma }_y\mathbb{H}\left({\gamma }_y,n{\gamma }_x,{\displaystyle \frac{{\beta }_x}{{\beta }_y}},0,1,0\right).$$

Figure 10 presents the $H_6$ PDFs generated using different parameter sets.

3. Application to Hydrological Data

In order to evaluate the proposed framework, we modeled hydrological datasets as extreme value distributions and compared them in an MSSR sense. We denoted $X_1,X_2,X_3,X_4$ and $X_5$ as the observed strength data and Y the observed stress data, respectively. The data are frequently used in the literature and are also presented below (e.g., [6,8,17]):

$$\mathbf{X}=\left(\begin{array}{ccccc}{X}_1& X_2& X_3& X_4& X_5\\ 0.42& 0.56& 0.73& 0.81& 0.83\\ 0.29& 0.36& 0.37& 0.46& 0.48\\ 0.54& 0.56& 0.72& 0.74& 0.75\\ 0.52& 0.61& 0.67& 0.72& 0.76\\ 0.34& 0.43& 0.47& 0.72& 0.74\\ 0.29& 0.39& 0.46& 0.61& 0.80\end{array}\right)\mathrm{and}Y=\left(\begin{array}{c}0.70\\ 0.65\\ 0.46\\ 0.72\\ 0.53\\ 0.77\end{array}\right).$$

The data represent the monthly water capacity of the Shasta Reservoir in California, USA, during the months of August and December from 1980 to 2015. Following [6,8,17], one can assume that drought conditions will be avoided if the water capacity in December is lower than the water capacities in at least two Augusts within the following five years. This means that the reliability condition is met if, in at least three of the next five years, the water capacities in August are not lower than the water capacity recorded in the preceding December. In this real-world setting, we consider a $s=3$ and $k=5$ MSSR scenario where stress–strength components are obtained as follows: $Y_1=0.70$ is the capacity of December 1980 and $X_1=(0.42,0.29,0.54,0.72,0.5310,0.77)$ are the capacities of August from 1981 to 1985; $Y_2=0.65$ is the capacity of December 1986 and $X_2=(0.56,0.36,0.56,0.61,0.43,0.39)$ are the capacities of August from 1987 to 1991, and so on. It is worth highlighting that the values presented have been divided by 3,014,878 (the mean of the water capacity); for more details about the dataset, the reader may refer to [6] and the references therein.

Aiming at estimating probabilities of the type $R_{s,k}$ given in (7) via the $I_n$ expressions given in the previous section, we need find the extreme distribution which best fits the data. It is important to compare different distributions as candidate models for the modeling of the dataset. We compared the performance of the l-max and p-max stable laws to the Rayleigh distributions as models for monthly water capacity. Maximum likelihood (ML) estimates and Kolmogorov–Smirnov (KS) p-values are given in

Table 1. ML-estimated parameters of distributions and the p-value of the Kolmogorov–Smirnov test.

RV	Models	$\widehat{\mathit{\mu }}$	$\widehat{\mathit{\sigma }}$	$\widehat{\mathit{\gamma }}$	KS p-Value
X	Fréchet	−28.80	29.29	194.27	0.34
X	Gumbel	0.49	0.15	-	0.34
X	RWeibull	0.85	0.31	1.70	0.47
X	GEV	0.54	0.18	−0.59	0.47
X	Rayleigh	0.24	7.44	-	0.28
RV	Models	$\widehat{\mathit{\alpha }}$	$\widehat{\mathit{\beta }}$	$\widehat{\mathit{\gamma }}$	KS $\mathit{p}$-Value
X	$H_1$	109.73	0.03	2.78	0.48
X	$H_2$	1.22	2.23	1.50	0.63
X	$H_5$	-	3.20	11.54	0.49
RV	Models	$\widehat{\mathit{\mu }}$	$\widehat{\mathit{\sigma }}$	$\widehat{\mathit{\gamma }}$	KS $\mathit{p}$-Value
Y	Fréchet	−20.84	21.42	197.44	0.71
Y	Gumbel	0.58	0.11	-	0.71
Y	RWeibull	0.77	0.13	0.71	0.82
Y	GEV	0.64	0.18	−1.40	0.87
Y	Rayleigh	0.41	15.85	-	0.68
RV	Models	$\widehat{\mathit{\alpha }}$	$\widehat{\mathit{\beta }}$	$\widehat{\mathit{\gamma }}$	KS $\mathit{p}$-Value
Y	$H_1$	39.43	0.14	2.93	0.63
Y	$H_2$	0.63	6.22	5.22	0.71
Y	$H_5$	-	5.36	20.17	0.63

. It is worth highlighting that estimating the parameters of extreme value distributions is not a straightforward task, as discussed in [19,20]. To do this, we considered the differential evolution method, a stochastic optimization strategy implemented in the python package Scipy [21] following the algorithm proposed by [22].

For each dataset, the empirical CDF (ECDF) is compared to the theoretical models in Figure 11 and Figure 12.

The Kolmogorov–Smirnov (KS) test was used to assess the adequacy of the theoretical cumulative distribution compared to the empirical distribution according to the following hypotheses:

• Null Hypothesis $\left(H_0\right)$: The sample data follow the specified theoretical distribution.
• Alternative Hypothesis $\left(H_1\right)$: The sample data do not follow the specified theoretical distribution.

A significance level of $0.05$ was used to determine whether to reject the null hypothesis. According to Table 1, the tests for all the models did not reject the null hypothesis, indicating that all of them are feasible modeling candidates according to the KS framework.

Overlaying the empirical cumulative distribution function (ECDF) with the theoretical fitted CDF provides a visual comparison between the observed and modeled cumulative distribution, allowing for an assessment of the goodness of fit. Notably, Figure 11 and Figure 12 demonstrate that, except for $H_1$ and $H_5$, the other distributions exhibit visually good fits. The same overall behavior was also observed for the variable Y.

Although some models exhibited similar visual modeling capabilities, the information criteria (Akaike information criterion—AIC; Bayesian information criterion—BIC; and efficient determination criterion—EDC) in

Table 2. Model selection using AIC, BIC, and EDC.

RV	Models	AIC	BIC	EDC
X	Fréchet	−14.61	−10.40	−16.93
X	Gumbel	−16.68	−13.88	−18.23
X	RWeibull	−20.01	−15.81	−22.34
X	GEV	−20.01	−15.81	−22.34
X	$H_1$	−8.99	−4.78	−11.31
X	$H_2$	−23.54	−19.34	−25.87
X	$H_5$	−11.22	−8.42	−12.77
X	Rayleigh	−20.05	−17.25	−21.60
Y	Fréchet	−2.36	−2.99	−6.61
Y	Gumbel	−4.38	−4.80	−7.22
Y	RWeibull	−14.13	−14.75	−18.38
Y	GEV	−14.13	−14.75	−18.38
Y	$H_1$	−1.42	−2.05	−5.67
Y	$H_2$	−26.98	−27.60	−31.23
Y	$H_5$	−3.55	−3.97	−6.39
Y	Rayleigh	−5.36	−5.78	−8.20

show that $H_2$ is the best model of the datasets, outperforming the Rayleigh model previously studied in [6]. In particular, EDC is a variation of both AIC and BIC such that

$$\begin{array}{ccc}\hfill AIC& =& 2k-2ln\left(\widehat{L}\right)\hfill \\ \hfill BIC& =& \ln \left(s\right)k-2ln\left(\widehat{L}\right)\hfill \\ \hfill EDC& =& \omega \ln \left(\ln \left(s\right)\right)k-2ln\left(\widehat{L}\right)\hfill \end{array}$$

(23)

where k is the number of estimated parameters in the model, $\widehat{L}$ is the maximized value of the likelihood function for the model, s is the sample size, and $\omega >0$ is a constant chosen to be 1 in the present paper.

Table 2 reveals that RWeibull and GEV had the same modeling capabilities, which is expected in cases of negative $\gamma$ values for GEV distributions. The top 3 distributions according to EDC would be $H_2$, GEV (RWeibull), and Rayleigh, but since this last distribution was already studied in a previous work [6], we included the Gumbel distribution as the third candidate model. Thus,

Table 3. Estimates of $R_{s,k}$ for some extreme-value models with $s=3$ and $k=5$.

Model	${\widehat{\mathit{R}}}_{3,5}$
Gumbel	0.30
GEV	0.33
$H_2$	0.31

shows the estimates for the MSSR of the monthly water capacity of the Shasta Reservoir in California in the analyzed period by considering $s=3$, $k=5$ and the Gumbel, GEV, and $H_2$ models in (7).

Table 3 shows that ${\widehat{R}}_{3,5}$ is not greater than $0.33$ for all models considered, which reveals the low reliability of the system (i.e., drought conditions are present). This should be used as a trigger to initiate contingency actions to mitigate the effects of possible droughts.

4. Conclusions

A multicomponent stress–strength analysis is a valuable tool for assessing hydrological scenarios. In particular, in the present paper, drought resilience in water systems was considered by evaluating the ability of resources to withstand varied stress factors indirectly associated with drought. This approach indirectly considers multiple parameters that influence water availability, such as precipitation patterns, reservoir levels, groundwater recharge rates, and usage demand by means of a proxy variable: water level. By comparing the combined “strength” (here taken as the water levels in particular months over years) of these factors to withstand stressors (also a particular water level), hydrologists can identify the likelihood of system failure (drought). Although it does not incorporate spatial data, this analysis is instrumental in modeling broad scenarios and quantifying drought-related risks. As a result, it aids in proactive decision-making for drought mitigation and resource allocation, ultimately enhancing the preparedness and sustainability of water systems during prolonged dry periods. In the present paper, the estimation of multicomponent stress–strength reliability when X and Y follow extreme value distributions was studied. New analytical expressions were obtained for several classes of extreme distributions, presenting a unified set of closed-form expressions in terms of the extreme-value $\mathbb{H}$-function. The expressions are easy to use, as was illustrated through their application to a real dataset. The good modeling performance of such expressions was shown, illustrating how these new results are useful to hydrologists. Specifically, the monthly water capacity data of the Shasta Reservoir in California (USA) during August and December from 1980 to 2015 were modeled, and the results obtained using the new expressions demonstrated a superior performance compared to previously published findings. In our simulation study, the primary objective was to conduct an exploratory analysis of the derived expressions and demonstrate their potential hydrological application. Therefore, exhaustive simulations were beyond the scope of this work. However, since our MSSR results are precise, an interesting direction for future simulation studies would involve assessing the impact on the MSSR of the parameter estimation procedures for each model. Thus, future research could benefit from developing bootstrap estimates for MSSR by incorporating the uncertainty associated with the maximum likelihood estimation of parameters. This approach could further refine reliability assessments in various hydrological scenarios and improve the robustness of MSSR estimates in practical applications.