Stress–Strength Reliability of the Type P(X < Y) for Birnbaum–Saunders Components: A General Result, Simulations and Real Data Set Applications

Abstract

An exact expression for $R=P(X<Y)$ has been obtained when X and Y are independent and follow Birnbaum–Saunders (BS) distributions. Using some special functions, it was possible to express R analytically with minimal parameter restrictions. Monte Carlo simulations and two applications considering real datasets were carried out to show the performance of the BS models in reliability scenarios. The new expressions are accurate and easy to use, making the results of interest to practitioners using BS models.

1. Introduction

The Birnbaum–Saunders (BS) fatigue-life model with two parameters was originally introduced in [1] and stems from the principles of renewal theory, focusing on the count of cycles required to surpass a critical threshold and induce fatigue crack propagation. In the same year, those same authors [2] provided maximum likelihood estimates (MLE) for BS model parameters. Later, Desmond [3] presented an alternative derivation of the distribution using a biological model and relaxing several of the assumptions made in [1], and Desmond [4] further studied how the Birnbaum–Saunders distribution could be related to inverse Gaussian distributions. The cumulative distribution function (CDF) of the Birnbaum–Saunders ($BS(\alpha ,\beta )$) model is given by:

$$F(t;\alpha ,\beta )=\Phi \left(\frac{1}{\alpha }\xi \left(\frac{t}{\beta }\right)\right),t>0,$$

(1)

where $\Phi (·)$ is the CDF of a standard normal distribution, $\alpha >0$ and $\beta >0$ are the shape and scale parameters, respectively, and

$$\xi \left(x\right)=\sqrt{x}-\frac{1}{\sqrt{x}},x>0.$$

(2)

In this work, the focus is on the stress–strength reliability metric $P(X<Y)$ (cf. [5]) when both X and Y are distributed as $BS(\alpha ,\beta )$ random variables. In broad terms, one may be interested in the likelihood of a system or component failing based on a comparison between the applied stress and the system’s strength. Let us consider stress, denoted as variable Y, and strength, denoted as variable X. Assuming these are independent continuous random variables (RVs) with probability density function (PDF) $f_Y$ and CDF $F_X$, respectively, the stress–strength probability (also referred to as reliability) is defined as follows:

$$R:=P(X<Y)={\int }_{-\infty }^{\infty }{F}_X\left(x\right)f_Y\left(x\right)dx.$$

This theory finds numerous applications, including statistics, economics, engineering and decision theory, and the reader may refer to the seminal work by Kotz et al. [5] for additional details.

For instance, when considering the reliability of a particular system, such a metric can indicate the likelihood of one system performing better than another concerning reliability or failure rates. Thus, specialists might utilise this metric to determine which system is more reliable for a specific application. Also, financial analysts might employ this metric to evaluate the relative risk linked with various investment strategies or assets. In contrast, manufacturing sectors could utilise this measure to appraise the effectiveness of quality control methods or production processes, thereby pinpointing areas for enhancement and optimising operations to enhance product quality and customer satisfaction.

In essence, the significance of reliability metrics such as R lies in their capacity to enable quantitative assessments for comparison, evaluation and informed decision making in intricate and uncertain circumstances. This, in turn, facilitates improved risk management, allocation of resources and optimisation across a broad spectrum of domains.

The stress–strength reliability metric when X and Y are independent RVs with distributions, respectively, $BS({\alpha }_x,{\beta }_x)$ and $BS({\alpha }_y,{\beta }_y)$, was studied by Xiuyun et al. [6], who derived an expression for the stress–strength probability R based on progressively Type II-censored samples. They presented an approximated expression for R when ${\alpha }_x={\alpha }_y=\alpha$. MLE and Bayesian estimators were also studied in [6]. In the present paper, it is of interest to explore cases without such parameter restrictions and in an exact and compact framework through the use of generalised hypergeometric special functions like the ones proposed by Rathie et al. [7].

Special functions are widely studied in the literature, such as the cases of the generalised hypergeometric function, Meijer’s G-function, Fox’s H-function, modified Bessel $K_{\nu }$-function of the third kind (see Definitions 1.1, 1.5, 1.6 and 1.11 in [8]), $\widehat{I}$-function ([7]). Thus, the results are hereby obtained using such functions and their properties. In the context of stress–strength reliability, the use of special functions is a way to provide expressions for R with fewer parameter restrictions (see [9] for example).

The paper is organised as follows: in Section 2, some preliminary concepts are presented, especially the definition of some useful special functions and of the Birnbaum–Saunders distribution (its density and properties of interest). Section 3 deals with the derivation of R when X and Y are independent Birnbaum–Saunders RVs as well as the estimation of R. In Section 4, Monte-Carlo simulations are presented to show the correctness of the analytical expressions derived. Besides, in Section 4, the modelling of two real datasets which account for the strength of different-length carbon fibres and daily wind speeds in two Atlantic coastal cities are presented. Then, conclusions are presented in Section 5.

2. Preliminaries

In this section, some definitions and results, which will be used subsequently, are presented.

2.1. Special Functions and Mellin Transforms

The Fox’s H-function (cf. Definition 1.1 in [8]) is defined by:

$$\begin{array}{c}\hfill H_{p,q}^{m,n}\left[z|\begin{array}{c}(a_1,A_1),\cdots ,(a_p,A_p)\\ (b_1,B_1),\cdots ,(b_q,B_q)\end{array}\right]=\\ \hfill \frac{1}{2\pi i}{\int }_L\frac{{\prod }_{k=1}^m\Gamma (b_j+B_{j}s){\prod }_{j=1}^n\Gamma (1-a_j-A_{j}s)}{{\prod }_{k=m+1}^q\Gamma (1-b_j-B_{j}s){\prod }_{j=n+1}^p\Gamma (a_j+A_{j}s)}{z}^{-s}ds,\end{array}$$

(3)

where $0\le m\le q$, $0\le n\le p$ (not both m and n simultaneously zero), $\Gamma (·)$ is the gamma function and i stands for square root of $(-1)$, $A_j>0$ ($j=1,\cdots ,p$), $B_k>0$ ($k=1,\cdots ,q$), $a_j$ and $b_k$ are complex numbers such that no poles of $\Gamma (b_k+B_{k}s)$ ($k=1,\cdots ,m)$ coincide with poles of $\Gamma (1-a_j-A_{j}s)$ ($j=1,\cdots ,n)$. L is a suitable contour $w-i\infty$ to $w+i\infty$, $w\in \mathbb{R}$, separating the poles of the two types mentioned above. For more details, see [8].

The Meijer’s G-function is a particular case of the H-function and can be defined (Definition 1.5 in 1.5 [8]) by:

$$G_{p,q}^{m,n}\left[z|\begin{array}{c}{a}_1,\cdots ,a_p\\ b_1,\cdots ,b_q\end{array}\right]=\frac{1}{2\pi i}{\int }_L\frac{{\prod }_{k=1}^m\Gamma (b_k-s){\prod }_{j=1}^n\Gamma (1-a_j+s)}{{\prod }_{k=m+1}^q\Gamma (1-b_k+s){\prod }_{j=n+1}^p\Gamma (a_j-s)}{z}^{-s}ds,$$

where $x\ne 0$, an empty product is interpreted as unity, $0\le m\le q$, $0\le n\le p$ (not both m and n simultaneously zero). The parameters $a_j$ ($j=1,\cdots ,n)$ and $b_k$ ($k=1,\cdots ,m)$ are such that no poles of ${\prod }_{j=1}^n\Gamma (b_k-s)$ coincide with poles of ${\prod }_{j=1}^n\Gamma (1-a_j+s)$. Refer to the study by [10] for details about the contour L and about the convergence conditions of the integral.

Let $\mathcal{M}\left[f\right]\left(s\right)$ denote the Mellin transform of a function f, which can be mathematically defined as (cf. [8]):

$$\mathcal{M}\left[f\right]\left(s\right)={\int }_0^{\infty }{t}^{s-1}f\left(t\right)dt,$$

(4)

provided that the integral converges. The inverse Mellin transform is obtained by the contour integral:

$$f\left(x\right)={\mathcal{M}}^{-1}\left[\mathcal{M}\left[f\right]\right]\left(x\right)=\frac{1}{2\pi i}{\int }_{c-i\infty }^{c+i\infty }{x}^{-s}\mathcal{M}\left[f\right]\left(s\right)ds.$$

Rathie et al. [7] proposed a generalisation of the H-function, thereby named the $\widehat{I}$ function, defined as a contour complex integral which contains H-functions in their integrands. The function is given by:

$${\widehat{I}}_m\left[z\left|\begin{array}{c}({\overline{a}}_{m,1},{\widehat{a}}_{m,1},{\overline{A}}_{m,1}),({\overline{a}}_{m,2},{\widehat{a}}_{m,2},{\overline{A}}_{m,2}),({\overline{a}}_{m,3},{\widehat{a}}_{m,3},{\overline{A}}_{m,3})\\ ({\overline{b}}_{m,1},{\widehat{b}}_{m,1},{\overline{B}}_{m,1}),({\overline{b}}_{m,2},{\widehat{b}}_{m,2},{\overline{B}}_{m,2})\\ ({\overline{\gamma }}_m,{\overline{\Gamma }}_m,{\overline{\pi }}_m,{\overline{\Pi }}_m,{\overline{\rho }}_m,{\overline{\sigma }}_m)\\ ({\overline{\alpha }}_m,{\overline{\beta }}_m,{\overline{\Lambda }}_m,{\overline{\Theta }}_m,{\overline{\zeta }}_m,{\overline{\eta }}_m)\end{array}\right.\right]=\frac{1}{2\pi i}{\int }_L\Psi \left(s\right)z^{-s}ds,$$

(5)

in which $\Psi \left(s\right)$ is the Mellin transform of the new function and can be explicitly given as:

$$\begin{array}{c}\hfill \Psi \left(s\right)=\prod _{j=1}^m{({\overline{\alpha }}_{j}s+{\overline{\beta }}_j)}^{{\overline{\Lambda }}_{j}s+{\overline{\Theta }}_j}{e}^{{\overline{\zeta }}_{j}s+{\overline{\eta }}_j}\times \\ \hfill \times {\left(H_{3,2}^{2,1}\left[\left.{({\overline{\gamma }}_j+{\overline{\Gamma }}_{j}s)}^{{\overline{\pi }}_j+{\overline{\Pi }}_{j}s}\right|\begin{array}{c}({\overline{a}}_{j,1}+{\widehat{a}}_{j,1}s,{\overline{A}}_{j,1}),({\overline{a}}_{j,2}+{\widehat{a}}_{j,2}s,{\overline{A}}_{j,2}),({\overline{a}}_{j,3}+{\widehat{a}}_{j,3}s,{\overline{A}}_{j,3})\\ ({\overline{b}}_{j,1}+{\widehat{b}}_{j,1}s,{\overline{B}}_{j,1}),({\overline{b}}_{j,2}+{\widehat{b}}_{j,2}s,{\overline{B}}_{j,2})\end{array}\right]\right)}^{{\overline{\rho }}_{j}s+{\overline{\sigma }}_j},\end{array}$$

(6)

where ${\overline{A}}_{j,k}$ and ${\overline{B}}_{j,k}$ are assumed to be positive real quantities, the ${\overline{a}}_{j,k}$, ${\widehat{a}}_{j,k}$, ${\overline{b}}_{j,k}$, ${\widehat{b}}_{j,k}$, ${\overline{\Pi }}_j$, ${\overline{\pi }}_j$, ${\overline{\Gamma }}_j$, ${\overline{\gamma }}_j$, ${\overline{\rho }}_j$, ${\overline{\sigma }}_j$, ${\overline{\alpha }}_j$, ${\overline{\beta }}_j$, ${\overline{\Lambda }}_j$, ${\overline{\Theta }}_j$, ${\overline{\zeta }}_j$, ${\overline{\eta }}_j$, $j=1,\cdots ,m$ are real numbers. The contour L runs from $c-i\infty$ to $c+i\infty$, where c is a real number, and exists in accordance with Mellin inversion theorem, taking into account all the singularities.

2.2. The BS Model

The PDF of the BS model is given by:

$$f\left(t\right)=\varphi \left(\frac{1}{\alpha }\xi \left(\frac{t}{\beta }\right)\right)\frac{1}{\alpha }{\xi }^{\prime }\left(\frac{t}{\beta }\right)\frac{1}{\beta },t>0,$$

(7)

which can also be explicitly expressed as:

$$f(t;\alpha ,\beta )=\frac{1}{2\sqrt{2\pi }\alpha \beta }\left[{\left(\frac{\beta }{t}\right)}^{1/2}+{\left(\frac{\beta }{t}\right)}^{3/2}\right]exp\left\{-\frac{1}{2{\alpha }^2}\left[\frac{t}{\beta }+\frac{\beta }{t}-2\right]\right\},t>0,$$

(8)

since $\varphi (·)$ denotes the PDF of a standard normal distribution. Figure 1 shows the PDF given in (7) and the CDF given in (1) for some choices of parameters.

From [11,12], it follows that:

$$E\left[T^p\right]={\beta }^{p}I(p,\alpha ),$$

(9)

where $T\sim BS(\alpha ,\beta )$ and

$$I(p,\alpha )=\frac{K_{p+1/2}\left({\alpha }^{-2}\right)+K_{p-1/2}\left({\alpha }^{-2}\right)}{2K_{1/2}\left({\alpha }^{-2}\right)},$$

(10)

and $K_{\nu }\left(z\right)$ is the modified Bessel function of the third kind (see Definition 1.11 in [8]). Also, literature [8] (p. 24) indicates that:

$$\begin{array}{ccc}\hfill K_{\nu }\left(z\right)& =& \frac{1}{2}{H}_{0,2}^{2,0}\left[\left.\frac{z^2}{4}\right|\begin{array}{cc}& \\ (\nu /2,1),& (-\nu /2,1)\end{array}\right]\hfill \\ \hfill & =& \frac{1}{2}{G}_{0,2}^{2,0}\left[\left.\frac{z^2}{4}\right|\frac{\nu }{2},-\frac{\nu }{2}\right].\hfill \end{array}$$

(11)

Equation (11) will be used to obtain an expression for the stress–strength probability R in terms of $\widehat{I}$ functions.

3. Stress–Strength Probability for BS Models

Let X and Y be independent random variables with $X\sim BS({\alpha }_x,{\beta }_x)$, and $Y\sim BS({\alpha }_y,{\beta }_y)$, ${\alpha }_j,{\beta }_j\in {\mathbb{R}}_{+}$ ($j\in \{x,y\}$). The stress–strength probability is given by:

$$R=P(X<Y)={\int }_0^{\infty }{F}_X(t;{\alpha }_x,{\beta }_x)f_Y(t;{\alpha }_y,{\beta }_y)dt.$$

(12)

As described, Xiuyun et al. [6] presented an approximate expression for (12) when ${\alpha }_1={\alpha }_2=\alpha$. The results hereby obtained, on the other hand, will not consider such restrictions.

Thus, it follows from (7) and (12) that:

$$\begin{array}{ccc}\hfill R& =& \hfill {\int }_0^{\infty }\Phi \left(\frac{1}{{\alpha }_x}\xi \left(\frac{t}{{\beta }_x}\right)\right)\varphi \left(\frac{1}{{\alpha }_y}\xi \left(\frac{t}{{\beta }_y}\right)\right)\frac{1}{{\alpha }_y}{\xi }^{\prime }\left(\frac{t}{{\beta }_y}\right)\frac{1}{{\beta }_y}dt\hfill \end{array}$$

At first, let us consider a special case where $\xi \left(\frac{t}{{\beta }_x}\right)=\xi \left(\frac{t}{{\beta }_y}\right)=u$. Then, ${\beta }_x={\beta }_y$, which implies:

$$R={\int }_{-\infty }^{\infty }\Phi \left(\frac{1}{{\alpha }_x}u\right)\varphi \left(\frac{1}{{\alpha }_y}u\right)\frac{1}{{\alpha }_y}du.$$

In this case, it is worth noticing that:

• $\Phi \left(\frac{1}{{\alpha }_x}u\right)$ is the CDF of $V_x\sim N(0,{\alpha }_x)$;
• $\varphi \left(\frac{1}{{\alpha }_y}u\right)$ is the PDF of $V_y\sim N(0,{\alpha }_y)$.

Thus, since $V_x$ and $V_y$ are independent:

$$R=P(V_x<V_y)=P(V_x-V_y<0)=\frac{1}{2},$$

which suggests the following remark:

Remark 1.

If ${\beta }_x={\beta }_x=\beta$, the stress–strength reliability metric in (12) is given by $R=0.5$.

Coming back to the definition of R given in (12), then:

$$R=P(X<Y)=P\left(\frac{X}{Y}<1\right),$$

since Y is a positive RV. One may look for the CDF of $W=\frac{X}{Y}$ to investigate R. The Mellin transform, as defined in (4), of the PDF $f_W$ is given by:

$$\mathcal{M}\left[f_W\right]\left(s\right)=\mathcal{M}\left[f_X\right]\left(s\right)\mathcal{M}\left[f_Y\right](2-s).$$

One may refer to [13] for further details on the properties of Mellin transforms in the context of the algebra of random variables.

From (9):

$$\mathcal{M}\left[f_X\right]\left(s\right)=E\left[X^{s-1}\right]={\beta }_x^{s-1}I(s-1,{\alpha }_x)$$

and, similarly,

$$\mathcal{M}\left[f_Y\right]\left(s\right)=E\left[Y^{s-1}\right]={\beta }_y^{s-1}I(s-1,{\alpha }_y).$$

Then, by inverting the Mellin transform:

$$f_W\left(z\right)=\frac{{\beta }_y}{{\beta }_x}\frac{1}{2\pi i}{\int }_{L}I(s-1,{\alpha }_x)I(1-s,{\alpha }_y){\left(\frac{{\beta }_{y}z}{{\beta }_x}\right)}^{-s}ds.$$

The direct integration of the PDF $f_W\left(z\right)$ leads to:

$$\begin{array}{ccc}\hfill F_W\left(z\right)& =& {\int }_0^z{f}_W\left(w\right)dw\hfill \\ \hfill & =& \frac{{\beta }_y}{{\beta }_x}\frac{1}{2\pi i}{\int }_{L}I(s-1,{\alpha }_x)I(1-s,{\alpha }_y){\left(\frac{{\beta }_y}{{\beta }_x}\right)}^{-s}{\int }_0^z{w}^{-s}dwds\hfill \\ \hfill & =& \frac{{\beta }_y}{{\beta }_x}\frac{1}{2\pi i}{\int }_{L}I(s-1,{\alpha }_x)I(1-s,{\alpha }_y){\left(\frac{{\beta }_y}{{\beta }_x}\right)}^{-s}\frac{z^{1-s}}{1-s}ds,\hfill \end{array}$$

provided that $Re(1-s)>0$ ($Re\left(z\right)$ is the real part of z). Then:

$$F_W\left(z\right)=z\frac{{\beta }_y}{{\beta }_x}\frac{1}{2\pi i}{\int }_L\frac{I(s-1,{\alpha }_x)I(1-s,{\alpha }_y)}{1-s}{\left(\frac{{\beta }_{y}z}{{\beta }_x}\right)}^{-s}ds,$$

which implies:

$$R=P\left(\frac{X_1}{X_2}<1\right)=F_W\left(1\right)=\frac{{\beta }_y}{{\beta }_x}\frac{1}{2\pi i}{\int }_L\frac{I(s-1,{\alpha }_x)I(1-s,{\alpha }_y)}{1-s}{\left(\frac{{\beta }_y}{{\beta }_x}\right)}^{-s}ds.$$

(13)

Using (10) and (11), the inner part of the integral (13) becomes:

$$\begin{array}{c}\hfill I(s-1,{\alpha }_x)I(1-s,{\alpha }_y)=\frac{1}{16}\frac{1}{K_{1/2}\left({\alpha }_x^{-2}\right)K_{1/2}\left({\alpha }_y^{-2}\right)}\\ \hfill \times \left\{G_{0,2}^{2,0}\left[\left.\frac{{\alpha }_x^{-4}}{4}\right|\frac{s-1/2}{2},-\frac{s-1/2}{2}\right]+G_{0,2}^{2,0}\left[\left.\frac{{\alpha }_x^{-4}}{4}\right|\frac{s-3/2}{2},-\frac{s-3/2}{2}\right]\right\}\\ \hfill \times \left\{G_{0,2}^{2,0}\left[\left.\frac{{\alpha }_y^{-4}}{4}\right|\frac{s-1/2}{2},-\frac{s-1/2}{2}\right]+G_{0,2}^{2,0}\left[\left.\frac{{\alpha }_y^{-4}}{4}\right|\frac{s-3/2}{2},-\frac{s-3/2}{2}\right]\right\}.\end{array}$$

This implies that R is the sum of four $\widehat{I}$ functions (5), proving the following result:

Theorem 1.

Let $X\sim BS({\alpha }_x,{\beta }_x)$ and $Y\sim BS({\alpha }_y,{\beta }_y)$ be independent RVs, ${\alpha }_x,{\alpha }_x,{\beta }_x,{\beta }_y\in {\mathbb{R}}_{+}$. Then:

$$\begin{array}{ccc}\hfill R& =& \hfill P(X<Y)=\frac{1}{16}\frac{{\beta }_y}{{\beta }_x}\frac{1}{K_{1/2}\left({\alpha }_x^{-2}\right)K_{1/2}\left({\alpha }_y^{-2}\right)}\left\{I_1+I_2+I_3+I_4\right\},\hfill \end{array}$$

(14)

where, for $k=1,\cdots ,4$:

$$I_k={\widehat{I}}_2\left[\frac{{\beta }_y}{{\beta }_x}\left|\begin{array}{c}\left\{\right(0,0,0),(1,0,0),(1,0,0\left)\right\},\left\{\right(0,0,0),(1,0,0),(1,0,0\left)\right\}\\ \{({\nu }_{1,k}/2,1/2,1),(-{\nu }_{1,k}/2,-1/2,1)\},\{({\nu }_{2,k}/2,1/2,1),(-{\nu }_{2,k}/2,-1/2,1)\}\\ \left\{({\alpha }_x^{-4}/4,0,1,0,0,1)\right\},\left\{({\alpha }_y^{-4}/4,0,1,0,0,1)\right\}\\ \left\{\right(-1,1,0,-1,0,0\left)\right\},\left\{\right(0,1,0,1,0,0\left)\right\}\end{array}\right.\right],$$

in which:

$${\nu }_{1,k}=\left\{\begin{array}{cc}-\frac{1}{2},& k\in \{1,2\},\\ -\frac{3}{2},& k\in \{3,4\},\end{array}\right.$$

and

$${\nu }_{2,k}=\left\{\begin{array}{cc}-\frac{1}{2},& k\in \{1,3\},\\ -\frac{3}{2},& k\in \{2,4\}.\end{array}\right.$$

Estimation

Let $X\sim BS({\alpha }_x,{\beta }_x)$ and $Y\sim BS({\alpha }_y,{\beta }_y)$ be independent RVs. It is of interest to estimate $R=P(X<Y)$ based on random samples of X and Y. For this, estimates of $\mathit{\theta }=({\alpha }_x,{\alpha }_y,{\beta }_x,{\beta }_y)$ are required. For simplicity of notation, let $\mathbf{x}=(x_1,\cdots ,x_n)$ and $\mathbf{y}=(y_1,\cdots ,y_m)$ denote observed random samples of X and Y, respectively.

The log-likelihood function, denoted by $l\left(\mathit{\theta }\right):=l(\mathit{\theta };\mathbf{x},\mathbf{y})$, is given by

$$\begin{array}{ccc}\hfill l\left(\mathit{\theta }\right)& =& -n(log{\alpha }_x+log{\beta }_x)-m(log{\alpha }_y+log{\beta }_y)\hfill \\ \hfill & +& \sum _{i=1}^n\left\{-\frac{1}{2{\alpha }_x^2}{\xi }^2\left(\frac{t_i}{{\beta }_x}\right)+log{\xi }^{\prime }\left(\frac{t_i}{{\beta }_x}\right)\right\}\hfill \\ & +& \sum _{j=1}^m\left\{-\frac{1}{2{\alpha }_y^2}{\xi }^2\left(\frac{u_j}{{\beta }_y}\right)+log{\xi }^{\prime }\left(\frac{u_j}{{\beta }_y}\right)\right\}.\hfill \end{array}$$

(15)

The maximum likelihood estimator (MLE) of $(\alpha ,\beta )$ was presented in [2]. In the present paper, two independent samples are considered with no parameter equality constraint (as ${\alpha }_x={\alpha }_y$ considered by [6]), such that $\widehat{\mathit{\theta }}$ is obtained directly by the MLEs $({\widehat{\alpha }}_x,{\widehat{\beta }}_x)$ and $({\widehat{\alpha }}_y,{\widehat{\beta }}_y)$.

Remark 2.

The MLE of R is obtained using the invariance property of MLE. This is due to the Theorem 1 that describes R in terms of $\widehat{I}$ functions (which are integrals, hence continuous and measurable functions), that is, using the estimates ${\widehat{\alpha }}_x$, ${\widehat{\alpha }}_y$${\widehat{\beta }}_x$ and ${\widehat{\beta }}_y$. Then, from (14):

$${\widehat{R}}_{MLE}=R({\widehat{\alpha }}_x,{\widehat{\alpha }}_y,{\widehat{\beta }}_x,{\widehat{\beta }}_y).$$

Algorithm 1 describes the approach used in the next section to obtain bootstrap confidence intervals of R.    

Algorithm 1: Let $\mathbf{x}$ and $\mathbf{y}$ be samples of sizes n and m, respectively, and a positive integer M.

Step 1  Generate independent bootstrap samples ${\mathbf{x}}_i$ and ${\mathbf{y}}_i$ from $\mathbf{x}$ and $\mathbf{y}$ (or directly sample from distributions if true parameters are known). Step 2  Compute the estimate ${\widehat{\mathit{\theta }}}_i=({\widehat{\alpha }}_{x,i},{\widehat{\alpha }}_{y,i},{\widehat{\beta }}_{x,i},{\widehat{\beta }}_{y,i})$ based on ${\mathbf{x}}_i$ and ${\mathbf{y}}_i$. Step 3  Obtain ${\widehat{R}}_i=\widehat{R}\left({\widehat{\mathit{\theta }}}_i\right)$ by (14). Step 4  Repeat steps $\mathbf{1}$ to $\mathbf{3}$ M times. Step 5  The approximate $100(1-\alpha )\%$ confidence interval of R is given by $[{\widehat{R}}_M(\alpha /2),{\widehat{R}}_M(1-\alpha /2)]$, where ${\widehat{R}}_M\left(\alpha \right)\approx {\widehat{G}}^{-1}\left(\alpha \right)$ and  $\widehat{G}$ is the cumulative distribution function of $\widehat{R}$ considering all $i=1,...,M$ samples.

It is also possible to calculate the mean of the estimates obtained in Algorithm 1, denoting this mean as ${\widehat{R}}_{boot}=M^{-1}{\sum }_{i=1}^M{\widehat{R}}_i$. It is clear that ${\widehat{R}}_{MLE}$ is the value of ${\widehat{R}}_i$ when $\mathbf{x}={\mathbf{x}}_i$ and $\mathbf{y}={\mathbf{y}}_i$, i.e., when the full sample is used to estimate the parameters of the distribution. In the following section, the estimates for both ${\widehat{R}}_{MLE}$ and ${\widehat{R}}_{boot}$ are obtained. When matched data are present ($m=n$), these results can be compared with a nonparametric estimator denoted as ${\widehat{R}}_{NP}$, which is defined as:

$${\widehat{R}}_{NP}=\frac{1}{n}\sum _{j=1}^n{\mathbb{1}}_{x_j\le y_j},$$

where ${\mathbb{1}}_A$ denotes the indicator function on the set A.

4. Applications

Monte Carlo simulations are presented to assess the estimation of R for the BS model as well as the modelling of two real data sets involving different-length carbon fibres and daily wind speeds in two Atlantic coastal cities. All the applications have been devised to show the full capabilities of the expressions hereby developed, and readers may inquire about the codes used directly to the corresponding author.

4.1. Monte Carlo Simulations

In this subsection, a Monte Carlo simulation study is carried out to evaluate the Bootstrap MLE ${\widehat{R}}_{boot}$, described in Algorithm 1. Additionally, the assessment of the estimates of ${\alpha }_x,{\beta }_x,{\alpha }_y$ and ${\beta }_y$ are also considered. The following fixed parameters are considered for the simulations:

• True population parameters $({\alpha }_x,{\beta }_x,{\alpha }_y)=(0.3,35,0.5)$ and ${\beta }_y\in \{27.5,30,32.5,37.5\}$;
• $M=100$ Monte Carlo replications of samples $\mathbf{x}$ and $\mathbf{y}$ directly sampled from true distributions since the parameters are known;
• sample sizes $m=n\in \{25,100,500,1000\}$.

Aiming to study the Monte Carlo MLEs, the bias and root mean squared error (RMSE) of each estimated value with respect to their true population’s values are computed. The simulation study was programmed in Python. The simulation results of ${\widehat{R}}_{boot}$ are presented in

Table 1. Monte Carlo simulation results for R estimation.

n	${\mathit{\alpha }}_{\mathit{x}}$	${\mathit{\beta }}_{\mathit{x}}$	${\mathit{\alpha }}_{\mathit{y}}$	${\mathit{\beta }}_{\mathit{y}}$	R	${\widehat{\mathit{R}}}_{\mathit{boot}}$	$\mathit{Bias}\left({\widehat{\mathit{R}}}_{\mathit{boot}}\right)$	$\mathit{RMSE}\left({\widehat{\mathit{R}}}_{\mathit{boot}}\right)$
25	0.3	35	0.5	27.5	0.6618	0.6694	0.0076	0.0060
25	0.3	35	0.5	30	0.6051	0.6092	0.0041	0.0071
25	0.3	35	0.5	32.5	0.5510	0.5599	0.0089	0.0075
25	0.3	35	0.5	37.5	0.4525	0.4484	−0.0042	0.0068
100	0.3	35	0.5	27.5	0.6618	0.6629	0.0011	0.0016
100	0.3	35	0.5	30	0.6051	0.6063	0.0012	0.0016
100	0.3	35	0.5	32.5	0.5510	0.5531	0.0021	0.0019
100	0.3	35	0.5	37.5	0.4525	0.4523	−0.0003	0.0017
500	0.3	35	0.5	27.5	0.6618	0.6619	0.0001	0.0003
500	0.3	35	0.5	30	0.6051	0.6058	0.0007	0.0003
500	0.3	35	0.5	32.5	0.5510	0.5516	0.0006	0.0004
500	0.3	35	0.5	37.5	0.4525	0.4522	−0.0003	0.0004
1000	0.3	35	0.5	27.5	0.6618	0.6614	−0.0004	0.0001
1000	0.3	35	0.5	30	0.6051	0.6048	−0.0003	0.0002
1000	0.3	35	0.5	32.5	0.5510	0.5516	0.0006	0.0002
1000	0.3	35	0.5	37.5	0.4525	0.4527	0.0002	0.0002

. Figure 2 and Figure 3 present the bias and RMSE for each estimated parameter. In short, each time Step 2 in Algorithm 1 is run, one ends up having estimates of the parameters ${\alpha }_x,{\beta }_x,{\alpha }_y$ and ${\beta }_y$. After $M=100$ runs, there are 100 values which are compared to the true population’s values mentioned before. Observe that by increasing the sample size n, the bias and RMSE of the Monte Carlo simulations decrease.

4.2. Real Data Set

Taking advantage of real data previously analysed in the literature, two applications are presented. The validity of the model for both datasets is studied and it is shown that the BS distribution presents a good fit in both cases.

4.2.1. Carbon Fibres

In the work of Bader and Priest [14], single carbon fibres were tested under tension at gauge lengths of 20 mm and 10 mm, having their strength measured in GPa. The data for both lengths are compared to assess which is statistically higher than the other. These data sets were modelled previously in several works (see, for example, [15]). For the convenience of the reader, data sets of length 20 mm ($\mathbf{x}$) and length 10 mm ($\mathbf{y}$) are presented below:

$$\begin{array}{ccc}\hfill \mathbf{x}& =& (1.312,1.314,1.479,1.552,1.700,1.803,1.861,1.865,1.944,1.958,1.966,1.977,2.006,\hfill \\ \hfill & & 2.021,2.027,2.055,2.063,2.098,2.140,2.179,2.224,2.240,2.253,2.270,2.272,2.274,\hfill \\ \hfill & & 2.301,2.301,2.359,2.382,2.382,2.426,2.434,2.435,2.478,2.490,2.511,2.514,2.535,\hfill \\ \hfill & & 2.554,2.566,2.570,2.586,2.629,2.633,2.642,2.648,2.684,2.697,2.726,2.770,2.773,\hfill \\ \hfill & & 2.800,2.809,2.818,2.821,2.848,2.880,2.954,3.012,3.067,3.084,3.090,3.096,3.128,\hfill \\ \hfill & & 3.233,3.433,3.585,3.585)\hfill \end{array}$$

and

$$\begin{array}{ccc}\hfill \mathbf{y}& =& (1.901,2.132,2.203,2.228,2.257,2.350,2.361,2.396,2.397,2.445,2.454,2.474,2.518,\hfill \\ \hfill & & 2.522,2.525,2.532,2.575,2.614,2.616,2.618,2.624,2.659,2.675,2.738,2.740,2.856,\hfill \\ \hfill & & 2.917,2.928,2.937,2.937,2.977,2.996,3.030,3.125,3.139,3.145,3.220,3.223,3.235,\hfill \\ \hfill & & 3.243,3.264,3.272,3.294,3.332,3.346,3.377,3.408,3.435,3.493,3.501,3.537,3.554,\hfill \\ \hfill & & 3.562,3.628,3.852,3.871,3.886,3.971,4.024,4.027,4.225,4.395,5.020).\hfill \end{array}$$

Descriptive statistics for $\mathbf{x}$ and $\mathbf{y}$ are presented in

Table 2. Descriptive statistics for breaking strength (in GPa) of carbon fibre at gauge lengths 20 mm ($\mathbf{x}$) and 10 mm ($\mathbf{y}$).

Data Set	Min.	1st Qu.	Median	Mean	3rd Qu.	Max.	Std. dv.
$\mathbf{x}$	1.31	2.10	2.48	2.45	2.77	3.58	0.50
$\mathbf{y}$	1.90	2.55	3.00	3.06	3.42	5.02	0.62

. The boxplot in Figure 4 shows that the data in $\mathbf{x}$ are a bit skewed to the right while the data in $\mathbf{y}$ are closer to being symmetric. Computing the value of the statistic R is important to quantitatively assess the differences observed in the data sets.

As the data sets have positive support, the BS distribution is a candidate to model such data sets. Let us assume that $\mathbf{x}$ and $\mathbf{y}$ are observed samples of $X\sim BS({\alpha }_x,{\beta }_x)$ and $Y\sim BS({\alpha }_y,{\beta }_y)$, respectively. Estimates of parameters ${\alpha }_x,{\alpha }_y$, ${\beta }_x$ and ${\beta }_y$ are presented in

Table 3. Estimated parameters, p-value of Kolmogorov–Smirnov (KS) test.

Data Set	$\widehat{\mathit{\alpha }}$	$\widehat{\mathit{\beta }}$	KS p-Value
$\mathbf{x}$	0.2138	2.3965	0.8496
$\mathbf{y}$	0.1984	3.0003	0.7678

as well as the p-value of the Kolmogorov–Smirnov (KS) test. The results indicate that it is not possible to reject the null hypothesis that BS models the CDF. Figure 5 shows the fit of distributions to data sets.

Using Theorem 1 and Algorithm 1, the values of ${\widehat{R}}_{MLE}=0.7802$, ${\widehat{R}}_{boots}=0.7808$ and the $95\%$ bootstrap confidence interval $(0.7092,0.8472)$ were obtained. The spread of the bootstrap estimates can be assessed in Figure 6. From the results, it is possible to say that since $P(X<Y)>0.5$ and $0.5$ is not within the CI of the Bootstrap estimates, carbon fibres with a length of 20 mm (sampled from X) show less strength when compared to carbon fibres with a length of 10 mm (sampled from Y).

As mentioned, this dataset has already been analysed in the stress–strength context previously. Valiollahi et al. [15] estimated R after transforming the original data (so that the transformed data had the same scale parameter) and modelling it using the Weibull distribution. They considered MLE, approximate MLE (AML) and Bayes estimator obtaining the values $0.5002$, $0.5172$ and $0.5221$, respectively. They also obtained the 95% Boot-p and Boot-t confidence intervals as $(0.4758,0.5789)$ and $(0.4726,0.5823)$. Their results cannot provide conclusive evidence of the true problem, as the transformation severely impairs the reliability calculations (now $P(X<Y)=0.5$ is within the CI). In the present analysis, such a transformation is not needed since the BS distribution was a good candidate for modelling the data and Theorem 1 does not require equality of parameters between $\mathbf{x}$ and $\mathbf{y}$. To compare the findings of the present study with their transformed datasets, using Tables 4 and 5 of [15], it was estimated that $X\sim BS(0.2138,0.9040)$, $Y\sim BS(0.1983,0.9051)$ and ${\widehat{R}}_{MLE}=0.4984$, which also sits on the reported confidence interval even though $BS$ distributions are considered.

4.2.2. Daily Wind Speeds

An application of stress–strength probability $R=P(X<Y)$ in the modelling and comparison of daily wind speeds (in 0.1 m/s units) in two Atlantic coastal cities, Coruña (Spain)—X—and Bergen (Norway)—Y—from 1 January 2010 till 31 December 2019, is presented. The data are presented below and were first analysed and studied in [16]:

$$\begin{array}{ccc}\hfill \mathbf{x}& =& (81,33,39,78,28,22,53,25,25,28,17,44,31,28,39,\hfill \\ \hfill & & 22,22,42,39,31,36,39,44,33,36,25,33,36,28,44)\hfill \end{array}$$

and

$$\begin{array}{ccc}\hfill \mathbf{y}& =& (36,57,26,52,29,93,50,72,53,11,31,27,37,15,28,\hfill \\ \hfill & & 38,28,26,48,17,34,28,35,30,45,22,100,126,21,39).\hfill \end{array}$$

Descriptive statistics for $\mathbf{x}$ (Coruña) and $\mathbf{y}$ (Bergen) are presented in

Table 4. Descriptive statistics for daily wind speeds (in 0.1 m/s units) in Coruña ($\mathbf{x}$) and Bergen ($\mathbf{y}$).

City	Min.	1st Qu.	Median	Mean	3rd Qu.	Max.	Std. dv.
Coruña	17.00	28.00	33.00	36.03	39.00	81.00	14.37
Bergen	11.00	27.25	34.50	41.80	49.50	126.00	25.99

. The boxplot in Figure 7 shows that the 50% highest $\mathbf{y}$ values appear to be larger and more dispersed than the 50% highest $\mathbf{x}$ values. Computing R can be useful to assess the suitability of building a wind power plant in either city.

Figure 8 shows the fit of BS models to the datasets; meanwhile, the estimated parameters are shown in

Table 5. Estimated parameters, p-value of Kolgomorov–Smirnov (KS) test for daily winter speeds in Coruña and Bergen.

City	Data Set	$\widehat{\mathit{\alpha }}$	$\widehat{\mathit{\beta }}$	KS p-Value
Coruña	$\mathbf{x}$	0.3463	34.0008	0.7125
Bergen	$\mathbf{y}$	0.5559	36.2300	0.8309

. Results of the Kolgomorov–Smirnov (KS) test indicate that BS random variables provide a reasonable representation of the datasets.

Using Theorem 1 with the parameters of Table 5, it was possible to obtain ${\widehat{R}}_{MLE}=0.5390$ from the fitted distributions, ${\widehat{R}}_{NP}=0.5333$ and ${\widehat{R}}_{boots}=0.5394$ ($M={10}^4$) whose confidence interval is $(0.3921,0.6825)$ at a $95\%$ significance level. It is important to highlight that the size of the bootstrap confidence interval was large, similar to the results obtained in [16] ($\widehat{R}=0.57$ and $CI=(0.45,0.69)$). One can see the distribution of the bootstrap estimates in Figure 9.

Since $P(X<Y)=0.5$ is within the CI, it is not possible to indicate which city would have higher coastal wind speeds, statistically.

5. Conclusions

Assessing reliability measures of the type $R=P(X<Y)$ has received much attention in several fields. When both X and Y follow Birnbaum–Saunders distributions, i.e., $X\sim BS({\alpha }_x,{\beta }_x)$ and $Y\sim BS({\alpha }_y,{\beta }_y)$, previous works have addressed the exact and approximate evaluation of R assuming several restrictions on the parameters of each distribution. In the present paper, R is obtained analytically in an exact and compact form in terms of generalised hypergeometric functions. This allows one to provide maximum likelihood estimates for R in a straightforward way. The results obtained indicate that the case when ${\beta }_x={\beta }_y=\beta$ is not interesting because it implies $R=0.5$. A Monte Carlo simulation indicated the correctness of the estimators proposed. Additionally, two real datasets were modelled using the new expression hereby derived, showing they are easy to use and accurate. Accurate and precise calculations of such metrics are essential for offering quantitative evaluations that aid in comparing, evaluating and making informed decisions in complex and uncertain situations. Consequently, this paper contributes to enhanced risk management, resource allocation and optimisation across various domains which consider Birnbaum–Saunders distributions to model the analysed random variables.