Many real-world systems encounter extremes in operation, which frequently result in malfunctions. Recently, there has been a lot of focus on the phenomena of systems failing to carry out their intended functions when they reach their lowest, highest, or both extreme operating conditions. This phenomenon is modeled by multi-stress-strength reliability \(R = \Pr \left( {W < Q < Z} \right),\) where system strength \((Q)\) is influenced by two random stresses ( \(W\) and \(Z\) ). Outliers, data points significantly deviating from the dataset's overall pattern, pose challenges in this context. This study address’s reliability estimation for the model \(R = \Pr \left( {W < Q < Z} \right),\) assuming that \(Q,W,\) and \(Z\) are independent inverted Kumaraswamy distributions. In the existence of outliers, maximum likelihood estimation for reliability \(R\) is computed, and the Bayesian approach employing independent gamma priors is also explored. It is possible to derive Bayesian estimators with symmetric or asymmetric loss functions. Markov chain Monte Carlo techniques are typically used to handle the complex computations required for this kind of computation. The outcomes of the numerical study show that both techniques' reliability estimates were enhanced by larger sample numbers in scenarios with and without outliers. In contrast, precision measures decreased in both scenarios with and without outliers as sample sizes increased. Remarkably, Bayesian estimates under the precautionary loss function frequently surpass those under other loss functions both in scenarios with and without outliers. Three real data sets are applied to the proposed methodology.
