1 Introduction

In the realm of reliability and life-testing experiments, a fundamental goal is to construct a theoretical framework for the empirical analysis of diverse physical phenomena. This framework is typically grounded in observational data, aiming to derive predictive inferences in a coherent manner. Numerous practical fields of study necessitate robust inferential analyses of various lifetime data. A key aspect of such investigations is to identify the underlying distribution that accurately describes the phenomenon of interest. The literature offers a plethora of distributions with flexible shape characteristics, suitable for modeling a wide range of lifetime data. For a comprehensive review of these distributions and their applications, readers may refer to Johnson et al. [1].

In statistical modeling and analysis, inverse distributions are essential, especially when working with data that are skewed or have large tails. These distributions are produced by reversing the original distribution’s order, which creates a new distribution with unique characteristics. To capture the distinctive qualities of some datasets, such as those found in industries like finance, insurance, and reliability engineering, this transformation is frequently required. Researchers can learn a great deal about the underlying patterns and behaviors of complicated data by understanding and employing inverted distributions. This can help them make more informed decisions and more accurate forecasts. Recently, Ghitany et al. [2] introduced and explored an inverse exponentiated class of distributions. This class of distributions has the following probability density function (PDF):

$${f_X}(x)=\frac{{{\upsilon _1}\omega }}{{{x^2}}}H^{\prime}\left( {\frac{1}{x}} \right)~~{\left( {1 - {e^{ - \omega H\left( {\frac{1}{x}} \right)}}} \right)^{{\upsilon _1} - 1}}{e^{ - \omega H\left( {\frac{1}{x}} \right)}};\,\,\,x,\omega,{\upsilon _1}>0,$$
(1)

where H(.) is an increasing function such that H(0) = 0 and H(∞) = ∞, while \(\omega,\) and \({\upsilon _1}\) are model parameters. This class encompasses several well-known models, including the inverted exponentiated exponential distribution, the inverted Burr X distribution, and the inverted exponentiated Pareto distribution (IEPD). Our interest here with the IEPD, presented by Ghitany et al. [2], by taking \(H\left( {\frac{1}{x}} \right)=\log \left( {1+\frac{1}{x}} \right)\) in Eq. (1), which yields the following PDF:

$${f_X}(x)={\upsilon _1}\omega {x^{\omega - 1}}~~{(1+x)^{ - (\omega +1)}}{\left[ {1 - {{\left( {1+\frac{1}{x}} \right)}^{ - \omega }}} \right]^{{\upsilon _1} - 1}};\,\,\,x>0,$$
(2)

where \(\omega>0,\) and \({\upsilon _1}>0\)are the shape parameters. The cumulative distribution function (CDF) and the hazard rate function (HF) associated with Eq. (2) are as follows:

$${F_X}(x)=1 - {\left[ {1 - {{\left( {1+\frac{1}{x}} \right)}^{ - \omega }}} \right]^{{\upsilon _1}}};\,\,\,x>0,$$
(3)

and,

\(h(x)=\omega {\upsilon _1}{x^{\omega - 1}}{(1+x)^{ - (\omega +1)}}{\left[ {1 - {{\left( {1+\frac{1}{x}} \right)}^{ - \omega }}} \right]^{ - 1}};\,\,\,x>0.\)

Figure 1 illustrates the HF of the IEPD for various parameter values. A visual examination reveals that the HF decreases rapidly over time when both shape parameters are less than one. Moreover, the function can exhibit non-monotonic behavior, suggesting that the IEPD may be suitable for modeling datasets with decreasing or non-monotonic HF. Consequently, potential applications of the IEPD include fields such as mortality data analysis, the study of mechanical or electrical component working conditions, degradation experiments, and fatigue failure. This distribution’s versatility positions it as a valuable tool in reliability and lifetime studies, comparable to other common models like the inverted exponential Rayleigh distribution, exponentiated moment exponential distribution, and generalized exponential distribution (see Pradhan and Kundu [3]). Studies concerning the IEPD have been discussed by several authors. Maurya et al. [4] investigated the estimation of the unknown parameters of an IEPD under progressive Type-II censoring. Shaikh and Patel [5] examined estimation of the IEPD parameters based on record values. Kumari et al. [6] focused on estimating the multicomponent stress-strength (MC-SS) reliability of the IEPD under Type II censoring.

Fig. 1
figure 1

The PDF and HF of the IEPD for selected parameter values

Full size image

In today’s highly competitive market, product reliability is a paramount concern. To ensure the quality and longevity of products, extensive life-testing and reliability experiments are conducted both before and during product launch. However, due to time constraints and other limitations, it is often impractical to observe the failure times of all units in a life test. Censoring schemes, which allow for the termination of the test before all units fail, are therefore essential tools in reliability analysis. Type-II censoring is the most commonly used test criterion and is advantageous for saving time and money. The progressive censoring (PC) scheme, a generalization of Type-II censoring, is particularly useful in scenarios where test units are removed during the experiment, not just at the final termination point. For a detailed understanding of the PC technique, readers are encouraged to consult works of Balakrishnan and Aggarwala [7]. Johnson [8] proposed a life test where the experimenter can group test units into sets and test each set simultaneously until the first failure occurs. Following a similar approach to Balasooriya [9], first-failure censoring can be employed for products with extended lifespans and limited testing capabilities. This scheme entails testing h × m units in m groups. The primary limitation of first-failure censoring is its inability to remove units at times other than the final termination. To address this issue, Wu and Kus [10] proposed progressively first-failure censored sampling (PFIF-CS), a life-testing strategy that combines first-failure and PC Type-II plans. The PFIF-CS allows for the removal of some sets of units from the life test before observing the first failures in all sets. Many authors, including Singh and Tripathi [11], Saini et al. [12], Chaturvedi et al. [13], Ramadan et al. [14], Mohammed et al. [15], Kayal et al. [16], Yousef and Almetwally [17], Cai and Gui [18], and Alsadat et al. [19], have examined this technique under various distributions. A description of the PFIF-CS is provided in Sect. 2.

The stress-strength (SS) model is a fundamental tool in mechanical engineering, used to evaluate the reliability and performance of system equipment. According to the model, the relationship between a system’s strength and the stress it experiences determines how reliable it is. A widely used formula in scientific texts is \(\gamma\)= P(Y < X), which represents the probability of the stress (Y) being less than the strength (X) of a material. Early research on SS reliability estimation was pioneered by Birnbaum [20]. Subsequent advancements in this field have been made by Nadar et al. [21], Najarzadegan et al. [22], Kundu and Gupta [23], and for more recent studies see ( [24,25,26,27,28,29,30]).

The MC-SS system is one that comprises multiple components. It can be viewed as an extension of the SS model, consisting of s independent and identically distributed (iid) strength variables, X₁,., Xs, and a common random stress Y that functions properly when at least r (≤ s) of the strength variables exceeds the random stress Y. This is referred to as the r-out-of-s: MC-SS model. Inference about multicomponent reliability is valuable in various practical applications, including bridge structures, communication systems, and military operations. For instance, a bridge with s vertical cables can be considered an MC-SS system. The bridge remains functional as long as the stress from factors like heavy traffic or wind load does not exceed the strength of at least r of its s cables. An automobile engine with eight cylinders can run if at least four of them are firing. This system is a 4-out-of-8: MC-SS system. Another example is an aircraft with two engines. The aircraft can still fly if only one engine is operational, which is a 1-out-of-2: MC-SS system. The reliability in the MC-SS system is defined by Bhattacharyya and Johnson [31] as follows:

$$\begin{gathered} {\eta _{r,s}}=P({\text{at}}\,{\text{least }}r{\text{ of the (}}{X_{\text{1}}},.,{X_s})\,{\text{exceed}}\,Y) \hfill \\ \,\,\,\,\,\,\,\,\,=\sum\limits_{{j=r}}^{s} {\left( {\begin{array}{*{20}{c}} s \\ j \end{array}} \right)} \int\limits_{{ - \infty }}^{\infty } {{{(1 - {F_X}(y))}^j}} {({F_X}(y))^{s - j}}d{F_Y}(y), \hfill \\ \end{gathered}$$
(4)

where X1,X2,…, Xs are iid random variables, where FX (.) and FY (.) are the CDF of strengths and stress random variables, respectively.

The MC-SS reliability models have gained significant attention in recent years, with numerous scholars investigating inference methods for various lifetime distributions and various sampling techniques. For example, Rao et al. [32] studied estimation of the MC-SS reliability for the Burr XII distribution. Dey et al. [33] conducted a comprehensive study of the MC-SS reliability for the Kumaraswamy distribution. References [34,35,36] have discussed the estimation of MC-SS reliability based on the exponentiated Pareto distribution. Kayal et al. [37] focused on the MC-SS for the Chen distribution. Several studies have investigated MC-SS reliability estimation via PC of Type-II for different distributions ( [38,39,40,41,42,43]). Kızılaslan [44] studied the MC-SS reliability for a family of inverse exponentiated distributions. The MC-SS reliability of Weibull and generalized inverse exponential distributions was extensively explored by Hassan et al. [45, 46] using record values. The estimation of the MC-SS reliability for the power Lomax distribution was thoroughly investigated by Ahmed et al. [47]. Little studies considered the estimation of MC-SS reliability based on PFIF-CS. These studies have considered certain underlying distributions, including Bur XII (Saini et al. [48]), the Kumaraswamy family of distributions (Saini and Garg [49]), and the generalized inverse exponential (Kumari et al. [50] and Saini [51]). For further insights into these inference methods, readers may refer to [52, 53].

This research explores statistical techniques to assess MC-SS reliability for the IEPD within the PFIF-CS framework. Our motivation for this research stems from the following:

  1. (i)

    The PFIF-CS offers significant time and cost savings compared to other censoring methods.

  2. (ii)

    There is a limited body of research, as mentioned above, based on MC-SS reliability models.

  3. (iii)

    The IEPD is a versatile and applicable distribution in various fields.

  4. (iv)

    This work can be considered an extension of the work presented by Kumari et al. [6], which focused on Type II censoring.

This paper is organized as follows: Sect. 2 presents the mathematical expression for the MC-SS reliability parameter, along with maximum likelihood estimates (MLEs) and asymptotic confidence intervals (As-CIs). Section 3 discusses Bayesian estimation, including Bayes estimates (BEs) and highest posterior density (HPD) credible intervals. Details of the numerical simulation results conducted to evaluate the proposed approach are given in Sect. 4. Section 5 demonstrates the applicability of the method using two real-world datasets. Section 6 concludes the paper by summarizing key findings and discussing potential future research directions.

2 Maximum Likelihood Estimator of ηr,s

This section gives an analytical expression of the reliability in the MC-SS model based on IEPD. The point and interval estimators of the model parameters and \({\eta _{r,s}}\)estimator are generated under the PFIF-CS in this section.

2.1 The MC-SS Reliability Expression

Here, the MC-SS reliability formula of the IEPD is determined. Suppose that X1, X2,…, Xs, and Y are independent random variables from IEPD\(({\upsilon _1},\omega )\) and IEPD\(({\upsilon _2},\omega ),\) respectively. Then the MC-SS expression is produced by using Eqs. (2), (3), and (4)

$${\eta _{r,s}}=\sum\limits_{{j=r}}^{s} {\left( {\begin{array}{*{20}{c}} s \\ j \end{array}} \right)} \int\limits_{0}^{\infty } {{{\left[ {1 - {{\left( {\frac{1}{y}+1} \right)}^{ - \omega }}} \right]}^{{\upsilon _1}j+{\upsilon _2} - 1}}{{\left[ {1 - {{\left[ {1 - {{\left( {\frac{1}{y}+1} \right)}^{ - \omega }}} \right]}^{{\upsilon _1}}}} \right]}^{s - j}}{\upsilon _2}\omega {y^{\omega - 1}}~~{{(1+y)}^{ - (\omega +1)}}} dy.$$
(5)

Using the binomial expansion in Eq. (5) gives

$${\eta _{r,s}}={\upsilon _2}\omega \sum\limits_{{j=r}}^{s} {\sum\limits_{{i=0}}^{{s - j}} {\left( {\begin{array}{*{20}{c}} s \\ j \end{array}} \right)} } \left( {\begin{array}{*{20}{c}} {s - j} \\ i \end{array}} \right){( - 1)^i}\int\limits_{0}^{\infty } {{{\left[ {1 - {{\left( {\frac{1}{y}+1} \right)}^{ - \omega }}} \right]}^{{\upsilon _1}(i+j)+{\upsilon _2} - 1}}{y^{\omega - 1}}~~{{(1+y)}^{ - (\omega +1)}}} dy.$$
(6)

Then after the simplified form of Eq. (6), we have

$${\eta _{r,s}}=\sum\limits_{{j=r}}^{s} {\sum\limits_{{i=0}}^{{s - j}} {\left( {\begin{array}{*{20}{c}} s \\ j \end{array}} \right)} } \left( {\begin{array}{*{20}{c}} {s - j} \\ i \end{array}} \right)\frac{{{{( - 1)}^i}{\upsilon _2}}}{{{\upsilon _1}(i+j)+{\upsilon _2}}}.$$
(7)

It’s important to note that expression (7) depends on the values of the parameters \({\upsilon _1},\)and \({\upsilon _2}.\) The 3D plots of the MC-SS reliability are shown in Fig. 2.

Fig. 2
figure 2

3D plots of \({\eta _{r,s}}\) at r = 2,3,4 and s = 5

Full size image

Figure 2 illustrates the impact of varying r values on the MC-SS reliability while s is held constant at 5. The three subplots represent the reliability values as functions of the parameters \({\upsilon _1}\) and\({\upsilon _2},\) with each subplot corresponding to different values of r (r = 2, 3, 4). The color gradient, ranging from red to yellow, indicates the reliability levels, where red color represents lower reliability values and yellow color represents higher values. The plots show that as \({\upsilon _1}\) and\({\upsilon _2}\)increase, the reliability tends to improve, as indicated by the transition from red to yellow. This suggests that higher values of the parameters \({\upsilon _1}\) and\({\upsilon _2}\) positively impact the system’s reliability. The surface plots effectively demonstrate how the reliability changes across different parameter settings, providing valuable insights for understanding the performance of the MC-SS model under varying conditions.

2.2 Maximum Likelihood Estimators

Here to obtain the MLE of the MC-SS reliability under PFIF-CS, we will begin by providing a concise description of the PFIF-CS.

Consider a life test where M ⋅ h items are divided into M independent groups of h items each. The test continues until a predetermined number, m, of failures occur. Upon observing the first failure at time X1:m: M:h,, Q1 live groups (including the one with the failure) are removed. This process is repeated for subsequent failures: at time X2:m: M:h, Q2 live groups are removed, and this pattern continues. Upon detecting the mth failure at time Xm: m:M: h, the remaining Qm live groups, along with the failed group, are discontinued from the study. The PFIF-CS data X1:m: M:h, X2:m: M:h,…, Xm: m:M: h, and a prefixed censoring scheme Q = (Q1, Q2,…,Qm), are obtained from this process, where M = m + Q1 + Q2+· · ·+Qm. The PFIF-CS plan encompasses various censoring schemes. For example, it can be adapted to:

  1. (i)

    At h = 1, the PFIF-CS reduces to PC of Type II.

  2. (ii)

    At Q = (0, 0,…,0), the PFIF-CS reduces to the first failure censoring plan.

  3. (iii)

    Under the conditions h = 1 and Q = (0, 0,…, M−m), the PFIF-CS simplifies to censoring of Type II.

  4. (iv)

    If h = 1 and Q = (0, 0,…,0), the PFIF-CS becomes equivalent to a complete sample.

Assume that the following PFIF-CS sample X1:m: M:h, X2:m: M:h,…, Xm: m:M: h is taken from a continuous population with CDF F(.) and PDF f(.). Then, according to Ref [10], the following represents the related likelihood function (LF):

$$f({x_1},.,{x_m})=d{h^m}\prod\limits_{{i=1}}^{m} {f({x_{i:m:M:h}})} {\left[ {1 - F({x_{i:m:M:h}})} \right]^{h({Q_i}+1) - 1}},$$
(8)

where \(d=\prod\limits_{{i=0}}^{m} {(M - i - \sum\limits_{{h=0}}^{i} {{Q_h}} } ),\)f(.) is the PDF, and F(.) is the CDF.

To derive the MLE of \({\eta _{r,s}},\) we will initially obtain MLEs for the parameters\({\upsilon _1},{\upsilon _2}\) and \(\omega\)using PFIF-CS. Let M × h1 systems, each system having S × h2 strength and 1 × h1 stress components, respectively, be put on a life testing experiment. After that, m systems, each system having s strength and one stress component, are observed under the PFIF-CS, and we get samples in the manner described below:

\(\quad \begin{array}{*{20}c} {{\text{Observed}}\,\,{\text{Strength}}\,{\text{Variables}}} \\ {\left( {\begin{array}{*{20}c} {X_{{11}} } & {X_{{12}} } & \cdots & {X_{{1s}} } \\ \vdots & \vdots & \ddots & \vdots \\ {X_{{m1}} } & {X_{{m2}} } & \cdots & {X_{{ms}} } \\ \end{array} } \right),} \\ \end{array} \quad {\text{and}}\quad \begin{array}{*{20}c} {{\text{Observed}}\,\,{\text{Stress}}\,{\text{Variables}}} \\ {\,\left( {\begin{array}{*{20}c} {Y_{1} } \\ \vdots \\ {Y_{m} } \\ \end{array} } \right)} \\ \end{array}\)

Let {Xi1, .., Xi s}, i = 1, .., m, be a PFIF-CS following the IEPD \(({\upsilon _1},\omega )\)with the PFIF-CS{S, s, h2, Q1, Q2, …, Qs } and {Y1, .., Ym} be an independent PFIF-CS following the IEPD \(({\upsilon _2},\omega )\) with the PFIF-CS {M, m, h2, \({Q^{\prime\prime}_1},....{Q^{\prime\prime}_m}\)}, the LF for both strengths and stress random variables, based on Eq. (8), is given by:

$$L={d_1}h_{1}^{m}\prod\limits_{{i=1}}^{m} {\left( {{d_2}h_{2}^{s}\prod\limits_{{j=1}}^{s} {f({x_{ij}}){{\left[ {1 - F({x_{ij}})} \right]}^{{h_1}({Q_j}+1) - 1}}} } \right)\,} f({y_i}){\left[ {1 - F({y_i})} \right]^{{h_2}({{Q^{\prime\prime}}_i}+1) - 1}},$$
(9)

The following is an expression for the LF of the observed data:

$$\begin{gathered} L\left( {data\left| \Theta \right.} \right)={d_2}h_{2}^{m}\prod\limits_{{i=1}}^{m} {\left( {{d_1}h_{1}^{s}\prod\limits_{{j=1}}^{s} {{\upsilon _1}\omega {x_{ij}}^{{\omega - 1}}~~{{(1+{x_{ij}})}^{ - (\omega +1)}}{{\left[ {A({x_{ij}},\omega )} \right]}^{{\upsilon _1}}}^{{{h_1}({Q_j}+1) - 1}}} } \right)\,} \hfill \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \times {\upsilon _2}\omega {y_i}^{{\omega - 1}}~~{(1+{y_i})^{ - (\omega +1)}}{\left[ {B({y_i},\omega )} \right]^{{\upsilon _2}}}^{{{h_2}({{Q^{\prime\prime}}_i}+1) - 1}}. \hfill \\ \end{gathered}$$
(10)

where \(\Theta ={({\upsilon _1},{\upsilon _2},\omega )^T}\)is the set of parameters, \(A({x_{ij}},\omega )=1 - {\left( {1+\frac{1}{{{x_{ij}}}}} \right)^{ - \omega }},\)and \(B({y_i},\omega )=1 - {\left( {1+\frac{1}{{{y_i}}}} \right)^{ - \omega }}.\)

Thus, the logarithm of the Eq. (10) is given by:

\(\begin{gathered} {L^ * } \propto ms\log (\omega {\upsilon _1})+\sum\limits_{{i=1}}^{m} {\sum\limits_{{j=1}}^{s} {\left[ {(\omega - 1)\log {x_{ij}} - (\omega +1)\log (1+{x_{ij}})} \right]+} } \sum\limits_{{i=1}}^{m} {\left[ {(\omega - 1)\log {y_i} - (\omega +1)\log (1+{y_i})} \right]} \hfill \\ \,\,\,\,\,\,\,+\sum\limits_{{i=1}}^{m} {\sum\limits_{{j=1}}^{s} {[{\upsilon _1}{h_1}({Q_j}+1) - 1]} } \,\log \left[ {A({x_{ij}},\omega )} \right]\,+m\log (\omega {\upsilon _2})+\sum\limits_{{i=1}}^{m} {[{h_2}{\upsilon _2}({{Q^{\prime\prime}}_i}+1) - 1]} \log \left[ {B({y_i},\omega )} \right]. \hfill \\ \end{gathered}\)

The following normal equations are employed to derive the MLEs of\({\upsilon _1},{\upsilon _2},\) and \(\omega\)

$$\frac{{\partial {L^ * }}}{{\partial {\upsilon _1}}}=\frac{{ms}}{{{\upsilon _1}}}+\sum\limits_{{i=1}}^{m} {\sum\limits_{{j=1}}^{s} {[{h_1}({Q_j}+1)]} } \,\log \left[ {A({x_{ij}},\omega )} \right]=0,$$
(11)
$$\frac{{\partial {L^ * }}}{{\partial {\upsilon _2}}}=\frac{m}{{{\upsilon _2}}}+\sum\limits_{{i=1}}^{m} {[{h_2}({{Q^{\prime\prime}}_i}+1)]} \log \left[ {B({y_i},\omega )} \right]=0,$$
(12)
$$\begin{aligned} \frac{{\partial L^{*} }}{{\partial \omega }} & = \frac{{m(s + 1)}}{\omega } + \sum\limits_{{i = 1}}^{m} {\sum\limits_{{j = 1}}^{s} {\left[ {\log x_{{ij}} - \log (1 + x_{{ij}} )} \right] + } } \sum\limits_{{i = 1}}^{m} {\left[ {\log y_{i} - \log (1 + y_{i} )} \right]} \\ & + \sum\limits_{{i = 1}}^{m} {\frac{{[h_{2} \upsilon _{2} (Q^{{\prime \prime }} _{i} + 1) - 1]B^{\prime } (y_{i},\omega )}}{{B(y_{i},\omega )}}} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} {\mkern 1mu} + \sum\limits_{{i = 1}}^{m} {\sum\limits_{{j = 1}}^{s} {\frac{{[\upsilon _{1} h_{1} (Q_{j} + 1) - 1]A^{\prime } (x_{{ij}},\omega )}}{{A(x_{{ij}},\omega )}}} }, \\ \end{aligned}$$
(13)

where, \({A^{\prime}_\omega }({x_{ij}},\omega )={\left( {1+\frac{1}{{{x_{ij}}}}} \right)^{ - \omega }}\log \left( {1+\frac{1}{{{x_{ij}}}}} \right),\) and \({B^{\prime}_\omega }({y_i},\omega )={\left( {1+\frac{1}{{{y_i}}}} \right)^{ - \omega }}\log \left( {1+\frac{1}{{{y_i}}}} \right).\)

Equations (11)–(13) can be used to create the MLEs \({\hat {\upsilon }_1},{\hat {\upsilon }_2},\) and \(\hat {\omega }\) by using the Newton-Raphson approach or any suitable iterative method. Consequently, the MLE \({\hat {\eta }_{r,s}}\)of \({\eta _{r,s}}\) may be created by utilizing the invariance characteristic of MLEs by merging \({\hat {\upsilon }_1},\) and \({\hat {\upsilon }_2}\) in Eq. (7) as follows:

$${\eta _{r,s}}=\sum\limits_{{j=r}}^{s} {\sum\limits_{{i=0}}^{{s - j}} {\left( {\begin{array}{*{20}{c}} {s - j} \\ i \end{array}} \right)} } \left( {\begin{array}{*{20}{c}} s \\ j \end{array}} \right)\frac{{{{( - 1)}^i}{{\hat {\upsilon }}_2}}}{{{{\hat {\upsilon }}_2}+{{\hat {\upsilon }}_1}(i+j)}}.$$
(14)

2.3 Asymptotic Confidence Intervals

Based on the asymptotic normality of the MLE, we construct the As-CIs of all parameters and MC-SS reliability \({\eta _{r,s}}.\) To build the As-CIs for \({\upsilon _1},{\upsilon _2},\) and \(\omega,\) we use the Fisher information matrix (FM)\(\hat {I}=E{\left[ {\frac{{ - \partial {L^ * }}}{{\partial {\tau _1}\partial {\tau _2}}}} \right]_{\Theta =\hat {\Theta }}},i,j{\text{ =1,2,3}}\) and \(\Theta =({\tau _1},{\tau _2},{\tau _3})=({\upsilon _1},{\upsilon _2},\omega )\)as described by Cohen [54]. Indeed, determining the exact formula for expectation can be challenging. After taking out the expectation operator E and replacing it with the observed variance–covariance matrix (the asymptotic FM) for the MLE \(\hat {\Theta }\) for \(\Theta\)

\({\hat {I}^{ - 1}}(\Theta )=\left[ {\begin{array}{*{20}{c}} {\frac{{ - {\partial ^2}{L^ * }}}{{\partial \upsilon _{1}^{2}}}}&{\frac{{ - {\partial ^2}{L^ * }}}{{\partial \upsilon _{1}^{{}}\partial {\upsilon _2}}}}&{\frac{{ - {\partial ^2}{L^ * }}}{{\partial \upsilon _{1}^{{}}\partial \omega }}} \\ {\frac{{ - {\partial ^2}{L^ * }}}{{\partial \upsilon _{2}^{{}}\partial {\upsilon _1}}}}&{\frac{{ - {\partial ^2}{L^ * }}}{{\partial \upsilon _{2}^{2}}}}&{\frac{{ - {\partial ^2}{L^ * }}}{{\partial \upsilon _{2}^{{}}\partial \omega }}} \\ {\frac{{ - {\partial ^2}{L^ * }}}{{\partial \omega \partial \upsilon _{1}^{{}}}}}&{\frac{{ - {\partial ^2}{L^ * }}}{{\partial \omega \partial \upsilon _{2}^{{}}}}}&{\frac{{ - {\partial ^2}{L^ * }}}{{\partial {\omega ^2}}}} \end{array}} \right]_{{_{{{\upsilon _1}={{\hat {\upsilon }}_1}}}{,_{{\upsilon _2}={{\hat {\upsilon }}_2},\omega =\hat {\omega }}}}}^{{ - 1}}=\left[ {\begin{array}{*{20}{c}} {{I_{11}}}&{{I_{12}}}&{{I_{13}}} \\ {{I_{21}}}&{{I_{22}}}&{{I_{23}}} \\ {{I_{31}}}&{{I_{32}}}&{{I_{33}}} \end{array}} \right]_{{}}^{{ - 1}}\)

where, the diagonal of \({\hat {I}^{ - 1}}(\Theta )\)gives the asymptotic variances of MLEs \({\hat {\upsilon }_1},{\hat {\upsilon }_2},\) and \(\hat {\omega }\) respectively. Note that the second partial derivatives of the unknown parameters are provided in the appendix at the end of the paper. The (1−\(\alpha\)) 100% As-CIs for\(\Theta =({\upsilon _1},{\upsilon _2},\omega )\) are determined as follows

\(\hat {\Theta } \pm Z_{{^{{{\alpha \mathord{\left/ {\vphantom {\alpha 2}} \right. \kern-0pt} 2}}}}}^{ * }\sqrt {\widehat {{\operatorname{var} }}(\hat {\Theta })},\,\)

where\(Z_{{^{{{\alpha \mathord{\left/ {\vphantom {\alpha 2}} \right. \kern-0pt} 2}}}}}^{ * }\) represents the top of \({\alpha \mathord{\left/ {\vphantom {\alpha 2}} \right. \kern-0pt} 2}\) percentile for the usual normal distribution.

To construct As-CIs for \({\eta _{r,s}}\), we need to calculate their variances. We approximate these variances using the delta method, as outlined in [55]. This approach allows us to estimate the variances of ​the MC − SS​ using the following formulas:

\(\operatorname{var} ({\hat {\eta }_{r,s}})={[\Delta \,]^T}[\hat {V}][\Delta ],\,\)

where,\(\Delta =\left( {\frac{{\partial {\eta _{r,s}}}}{{\partial {\upsilon _1}}},\frac{{\partial {\eta _{r,s}}}}{{\partial {\upsilon _2}}},\frac{{\partial {\eta _{r,s}}}}{{\partial \omega }}} \right),\)\(\frac{{\partial {\eta _{r,s}}}}{{\partial {\upsilon _1}}}=\sum\limits_{{j=r}}^{s} {\sum\limits_{{i=0}}^{{s - j}} {\left( {\begin{array}{*{20}{c}} s \\ j \end{array}} \right)} } \left( {\begin{array}{*{20}{c}} {s - j} \\ i \end{array}} \right)\frac{{{{( - 1)}^{i+1}}{\upsilon _2}(i+j)}}{{{{\left[ {{\upsilon _1}(i+j)+{\upsilon _2}} \right]}^2}}},\,\)\(\frac{{\partial {\eta _{r,s}}}}{{\partial \omega }}=0,\)\(\frac{{\partial {\eta _{r,s}}}}{{\partial {\upsilon _2}}}=\sum\limits_{{j=r}}^{s} {\sum\limits_{{i=0}}^{{s - j}} {\left( {\begin{array}{*{20}{c}} s \\ j \end{array}} \right)} } \left( {\begin{array}{*{20}{c}} {s - j} \\ i \end{array}} \right)\frac{{{{( - 1)}^i}{\upsilon _1}(i+j)}}{{{{\left[ {{\upsilon _1}(i+j)+{\upsilon _2}} \right]}^2}}}.\)

Therefore, 100\((1 - \alpha )\)% a two-sided As-CI of the MC-SS reliability can be generated as:

\({\hat {\eta }_{r,s}} \pm Z_{{^{{{\alpha \mathord{\left/ {\vphantom {\alpha 2}} \right. \kern-0pt} 2}}}}}^{ * }\sqrt {\widehat {{\operatorname{var} }}({{\hat {\eta }}_{r,s}})}.\,\)

3 Bayesian Estimation

This section deals with the BEs of the parameters \({\upsilon _1},{\upsilon _2},\) and \(\omega,\) as well as the MC-SS reliability estimate of \({\eta _{r,s}}.\) Given the substantial influence of the loss function, we consider both the squared error loss function (SELF) as symmetric and the linear exponential loss function (LLF) as asymmetric to derive the BEs. In the present study, we assume that the parameters \({\upsilon _1},{\upsilon _2},\) and \(\omega,\) are independent random variables following gamma distributions, i.e., \({\upsilon _1}\sim\) Gamma (a1, b1), \({\upsilon _2}\sim\) Gamma (a2, b2), and \(\omega \sim\) Gamma (a3, b3). In Bayesian statistics, the gamma distribution as an informative prior (INP) is often chosen as a prior distribution due to several reasons: (i) It can take on various shapes, including exponential and right-skewed distributions, making it suitable for modeling a wide range of positive-valued random variables, (ii) For certain likelihood functions, the gamma distribution is a conjugate prior, leading to analytical solutions for posterior distributions; (iii) The gamma distribution is defined for non-negative values, making it appropriate for modeling parameters like rates, scales, and counts, (iv) Its shape and scale parameters have clear interpretations, allowing for easy representation of prior beliefs; and (v) The gamma distribution’s connection to the gamma function simplifies Bayesian computations, such as calculating posterior distributions and credible intervals. The joint PDF of \({\upsilon _1},{\upsilon _2},\) and \(\omega\) is given by:

$$\pi (\Theta ) \propto \upsilon _{1}^{{^{{{a_1} - 1}}}}\upsilon _{2}^{{{a_2} - 1}}{\omega ^{_{{{a_3} - 1}}}}{e^{ - ({b_1}{\upsilon _1}+{b_2}{\upsilon _2}+{b_3}\omega )}},\,\,\,\,\,\,{a_i},{b_i}>0,$$
(15)

where ai, and bi, i = 1, 2, 3 are the hyper parameters. Notably, prior distributions are referred to as non-informative (N-INP) when \({a_i},\) and \({b_i}\) take values approaching zero. When using INPs, the gamma priors provided to ai, and bi have means and variances that match the MLEs of \({\hat {\upsilon }_1},{\hat {\upsilon }_2},\) and \(\hat {\omega }\). Consequently, we set the variance and mean of the MLEs of parameters with the gamma priors, which leads us to the conclusion presented by Dey et al. [56].

\(\frac{1}{d}\sum\limits_{{j=1}}^{d} {\hat {\Theta }_{i}^{j}} =\frac{{{a_i}}}{{{b_i}}},\,\,\,\,\,\,\,\,\frac{1}{{d - 1}}{\sum\limits_{{j=1}}^{d} {\left[ {\hat {\Theta }_{i}^{j} - \frac{1}{d}\sum\limits_{{j=1}}^{d} {\hat {\Theta }_{i}^{j}} } \right]} ^2}=\frac{{{a_i}}}{{b_{i}^{2}}},\,\,\,i=1,2,3,\,\,\,\hat {\Theta }=({\hat {\upsilon }_1}{\text{,}}{\hat {\upsilon }_{\text{2}}}{\text{,}}\hat {\omega }{\text{), }}\)

After solving the two previously mentioned equations, the estimated hyper-parameters are expressed as follows:

\({a_i}=\frac{{{{\left[ {{d^{ - 1}}\sum\limits_{{j=1}}^{d} {\hat {\Theta }_{i}^{j}} } \right]}^2}}}{{\,{{(d - 1)}^{ - 1}}{{\sum\limits_{{j=1}}^{d} {\left[ {\hat {\Theta }_{i}^{j} - \frac{1}{d}\sum\limits_{{j=1}}^{d} {\hat {\Theta }_{i}^{j}} } \right]} }^2}}},\,\,\,\,\,\,\,{b_i}=\frac{{{d^{ - 1}}\sum\limits_{{j=1}}^{d} {\hat {\Theta }_{i}^{j}} }}{{\,{{(d - 1)}^{ - 1}}{{\sum\limits_{{j=1}}^{d} {\left[ {\hat {\Theta }_{i}^{j} - \frac{1}{d}\sum\limits_{{j=1}}^{d} {\hat {\Theta }_{i}^{j}} } \right]} }^2}}}.\)

The joint posterior distribution of \({\upsilon _1},{\upsilon _2},\) and \(\omega,\) is produced by employing Equations (15), and (10), one can derive

$$\Pi (\Theta \left| {data} \right.)={M^{ - 1}}\,L\left( {data\left| \Theta \right.} \right)\pi (\Theta ),$$
(16)

where \(M=\int\limits_{0}^{\infty } {\int\limits_{0}^{\infty } {\int\limits_{0}^{\infty } {L\left( {data\left| \Theta \right.} \right)\pi (\Theta )d({\upsilon _1}{\upsilon _2}\omega )} } }\) is a normalizing constant. The conditional posteriors are given as:

$${\Pi _1}({\upsilon _1}\left| {\omega,data} \right.) \propto \upsilon _{1}^{{ms+{a_1} - 1}}{e^{ - {\upsilon _1}({b_1} - \sum\limits_{{i=1}}^{m} {\sum\limits_{{j=1}}^{s} {\left[ {{h_1}({Q_j}+1)} \right]\log \left[ {A({x_{ij}},\omega )} \right]} } }},\,$$
(17)
$${\Pi _2}({\upsilon _2}\left| \omega \right.,data) \propto \upsilon _{2}^{{m+{a_2} - 1}}{e^{ - {\upsilon _2}\left[ {{b_2} - \sum\limits_{{i=1}}^{m} {{h_2}({{Q^{\prime\prime}}_i}+1)\log \left[ {B({y_i},\omega )} \right]} } \right]}},$$
(18)
$$\begin{gathered} {\Pi _3}(\omega \left| {{\upsilon _1},{\upsilon _2},data} \right.) \propto \omega _{{}}^{{m(s+1)+{a_3} - 1}}{e^{ - \omega \left( {{b_3} - \sum\limits_{{i=1}}^{m} {\sum\limits_{{j=1}}^{s} {\left[ {\log (1+{x_{ij}}) - \log {x_{ij}}} \right] - \sum\limits_{{i=1}}^{m} {\left[ {\log (1+{y_i}) - \log {y_i}} \right]} } } } \right)}} \hfill \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{e^{\sum\limits_{{i=1}}^{m} {\sum\limits_{{j=1}}^{s} {\left( {[{\upsilon _1}{h_1}({Q_j}+1) - 1]\log \left[ {A({x_{ij}},\omega )} \right]} \right)+\sum\limits_{{i=1}}^{m} {\left( {[{\upsilon _2}{h_2}({{Q^{\prime}}_i}+1) - 1]\log \left[ {B({y_i},\omega )} \right]} \right)} } } }}. \hfill \\ \end{gathered}$$
(19)

The Bayesian estimator of \({\eta _{r,s}},\) denoted by \(\tilde {\eta }_{{r,s}}^{{SELF}},\)\(\tilde {\eta }_{{r,s}}^{{LLF}},\)under SELF and LLF, respectively, are given, as follows:

$$\left. \begin{gathered} {L_{SELF}}\left( {{\eta _{r,s}},\tilde {\eta }_{{r,s}}^{{SELF}}} \right)={\left( {{\eta _{r,s}} - \tilde {\eta }_{{r,s}}^{{SELF}}} \right)^2},\,\,\,\,\,\,\,\,\,\,{L_{LLF}}\left( {{\eta _{r,s}},\tilde {\eta }_{{r,s}}^{{LLF}}} \right)={e^{c({\eta _{r,s}} - \tilde {\eta }_{{r,s}}^{{LLF}})}} - c({\eta _{r,s}} - \tilde {\eta }_{{r,s}}^{{LLF}}) - 1 \hfill \\ \tilde {\eta }_{{_{{r,s}}}}^{{SELF}}=E\left( {{\eta _{r,s}}} \right),\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\tilde {\eta }_{{_{{r,s}}}}^{{LLF}}=\frac{{ - 1}}{c}\ln \left[ {E\left( {{e^{ - c{\eta _{r,s}}}}} \right)} \right] \hfill \\ \end{gathered} \right\},$$
(20)

where c is an LLF scale parameter. As seen from Eqs. (17) and (18), that the \(({\upsilon _1}\left| \omega \right.,data)\) has a distribution with parameters \(\Big( {ms+{a_1},\,{b_1} - \sum\limits_{{i=1}}^{m} {\sum\limits_{{j=1}}^{s} {\left[ {{h_1}({Q_j}+1)} \right]} } } \)\(\Big. {{{ \left[ {A({x_{ij}},\omega )} \right]} } } \Big)\) and \(({\upsilon _2}\left| \omega \right.,data)\) has a gamma distribution with parameters\(\left( {m+{a_2},\,{b_2} - \sum\limits_{{i=1}}^{m} {\left[ {{h_2}({{Q^{\prime\prime}}_i}+1)} \right]\log \left[ {B({y_i},\omega )} \right]} } \right).\) Therefore, any gamma-generating method can easily create samples of \({\upsilon _1}\) and \({\upsilon _2}.\) However, the posterior conditional distribution of \(\omega\) is not analytically tractable in our case, Markov Chain Monte Carlo (MCMC). Also, it is evident that a closed form cannot be derived from the BE of \({\eta _{r,s}}\) provided in Eq. (20). Gibbs sampling and Metropolis-Hastings (M-H) algorithms are prominent MCMC techniques. The M-H, combined with Gibbs sampling, is a particularly popular approach. Like acceptance-rejection sampling, the M-H proposes candidate values and accepts or rejects them based on a probability. The MCMC is also used to construct credible intervals. This is a detailed flowchart of the M-H algorithm:

  1. 1.

    Initialize the algorithm with parameter values \(\:{\upsilon\:}_{10},{\upsilon\:}_{20},\) and \(\:{\omega\:}_{0}\) (as likelihood estimators: \(\:{\upsilon\:}_{10}={\widehat{\upsilon\:}}_{1},{{\upsilon\:}_{20}=\widehat{\upsilon\:}}_{2},\) and \(\:{\omega\:}_{0}=\widehat{\omega\:}\)).

  2. 2.

    Generate candidate parameter values \(\:{{\upsilon\:}_{1}}^{*},{{\upsilon\:}_{2}}^{*},\) and \(\:{\omega\:}^{*}\) in each iteration t using a proposal distribution \(\:\left({{\upsilon\:}_{1}}^{*}\right|{{\upsilon\:}_{1}}^{t-1}),q({{\upsilon\:}_{2}}^{*}\left|{{\upsilon\:}_{2}}^{t-1}\right),\) and \(\:q\left({\omega\:}^{*}|{\omega\:}^{t-1}\right).\) A major factor influencing the convergence and efficiency of the method is the proposal distribution selection. One typical tactic is to employ a symmetric proposal distribution, such as the normal distribution.

  3. 3.

    Obtain the Following Acceptance Probability

\(\:{\text{{\rm\:A}}}_{{\upsilon\:}_{1}}=min\left(1,\frac{{\varPi\:}_{1}\left({{\upsilon\:}_{1}}^{*}\left|{\upsilon\:}_{20},\right.{\omega\:}_{0}\right)\:q\left({{\upsilon\:}_{1}}^{t-1}\right|{{\upsilon\:}_{1}}^{*})}{{\varPi\:}_{1}\left({{\upsilon\:}_{1}}^{t-1}\left|{\upsilon\:}_{20},\right.{\omega\:}_{0}\right)\:q\left({{\upsilon\:}_{1}}^{*}\right|{{\upsilon\:}_{1}}^{t-1}}\right),{\text{{\rm\:A}}}_{{\upsilon\:}_{2}}=min\left(1,\frac{{\varPi\:}_{2}\left({{\upsilon\:}_{2}}^{*}\left|{\upsilon\:}_{1t},\right.{\omega\:}_{t}\right)\:q\left({{\upsilon\:}_{2}}^{t-1}\right|{{\upsilon\:}_{2}}^{*})}{{\varPi\:}_{2}\left({{\upsilon\:}_{2}}^{t-1}\left|{\upsilon\:}_{1t},\right.{\omega\:}_{t}\right)\:q\left({{\upsilon\:}_{2}}^{*}\right|{{\upsilon\:}_{2}}^{t-1}}\right),\) & \(\:{\text{{\rm\:A}}}_{\omega\:}=min\left(1,\frac{{\varPi\:}_{3}\left({\omega\:}^{*}\left|{\upsilon\:}_{10},{\upsilon\:}_{20}\right.\right)\:q\left({\omega\:}^{t-1}\right|{\omega\:}^{*})}{{\varPi\:}_{3}\left({\omega\:}^{t-1}|{\upsilon\:}_{10},{\upsilon\:}_{20}\right)\:q\left({\omega\:}^{*}\right|{\omega\:}^{t-1}}\right)\)

as the ratio of the target distribution relative to the proposal distribution:

  1. 4.

    Create a uniform distribution between 0 and 1 and use it to generate the random numbers \(\:{u}_{{\upsilon\:}_{1}},{u}_{{\upsilon\:}_{2}}\text{a}\text{n}\text{d}\:{u}_{\omega\:}\)

If \(\:{u}_{{\upsilon\:}_{1}}\le\:{\text{{\rm\:A}}}_{{\upsilon\:}_{1}}\) accept the proposed value: \(\:{\stackrel{\sim}{\upsilon\:}}_{1t}={{\upsilon\:}_{1}}^{*}\). & If \(\:{u}_{{\upsilon\:}_{1}}>{\text{{\rm\:A}}}_{{\upsilon\:}_{1}}\) reject the proposed value: \(\:{\stackrel{\sim}{\upsilon\:}}_{1t}={{\upsilon\:}_{1}}^{t-1}\).

If \(\:{u}_{{\upsilon\:}_{2}}\le\:{\text{{\rm\:A}}}_{{\upsilon\:}_{2}}\) accept the proposed value: \(\:{\stackrel{\sim}{\upsilon\:}}_{2t}={{\upsilon\:}_{2}}^{*}\). & If \(\:{u}_{{\upsilon\:}_{2}}>{\text{{\rm\:A}}}_{{\upsilon\:}_{2}}\) reject the proposed value: \(\:{\stackrel{\sim}{\upsilon\:}}_{2t}={{\upsilon\:}_{2}}^{t-1}\).

If \(\:{u}_{\omega\:}\le\:{\text{{\rm\:A}}}_{\omega\:}\) accept the proposed value: \(\:{\stackrel{\sim}{\omega\:}}_{t}={\omega\:}^{*}\). & If \(\:{u}_{\omega\:}>{\text{{\rm\:A}}}_{\omega\:}\) reject the proposed value: \(\:{\stackrel{\sim}{\omega\:}}_{t}={\omega\:}^{t-1}\).

  1. 5.

    Repetition of stages 1–4 is required for a predetermined number of iterations or until convergence. Start the next iteration with the accepted or current parameter values, \(\:{\stackrel{\sim}{\upsilon\:}}_{1t},{\stackrel{\sim}{\upsilon\:}}_{2t},\) and \(\:{\stackrel{\sim}{\omega\:}}_{t}\)as the starting point for each iteration.

In Bayesian statistics, terminology like the credible interval and HPD interval are frequently employed. A technique for computing HPD intervals was published by Chen and Shao [57]. It is specifically made to handle multimodal distributions and comes in handy when the posterior distribution has complex forms. The algorithm developed by Chen and Shao [57] is as follows: Sort the parameter values that were sampled according to their posterior distribution. For the unknown parameters or any function of them, the 95% two-sided HPD credible interval is as follows:

\(\:\left[{{\upsilon\:}_{1}}_{\left[250\right]},{{\upsilon\:}_{1}}_{\left[9750\right]}\right],\left[{{\upsilon\:}_{2}}_{\left[250\right]},{{\upsilon\:}_{2}}_{\left[9750\right]}\right]\) and \(\:\left[{\omega\:}_{\left[250\right]},{\omega\:}_{\left[9750\right]}\right]\).

4 Simulation Study

In this section, we conduct a simulation study to assess the MC-SS reliability of the IEPD under the PFIF-CS. Using the R programming language and packages like “bbmle,” “coda,” “HDInterval,” and “rootSolve,” we generate random samples from the IEPD based on the specified model and scheme. To evaluate the performance of the proposed methods and estimate MC-SS reliability, we execute a comprehensive simulation experiment. The IEPD-based simulations under PFIF-CS offer researchers a valuable tool for overcoming challenges such as data scarcity, complex relationships, and scenario analysis. By generating artificial data and exploring various scenarios, these simulations can enhance our understanding of the model’s behavior and improve its robustness in diverse applications. Using Monte Carlo simulation, we assess the parameter performance and dependability of the MC-SS model with varying sample sizes. We vary r over 2, 3, and 4, keeping s fixed at 5. Table 1 lists the sample sizes and censored samples, and a total of 10,000 random full samples are generated.

Table 1 Selection of complete and censored sample sizes along with corresponding codes
Full size table

This section presents an empirical evaluation of the proposed methods’ performance across various sample sizes, based on Monte Carlo simulations. We choose the values of the true parameters based on several key considerations to ensure the robustness and relevance of our analysis. First, the parameter values were selected to reflect practical scenarios commonly observed in reliability studies, thus ensuring that our findings would be applicable in real-world situations. We also reviewed existing literature on the IEPD and other reliability models to identify typical parameter ranges, which helped us choose values that allow for meaningful comparisons with prior work. Additionally, we aimed to test the model’s performance across a diverse set of conditions, including different levels of stress and strength, to thoroughly evaluate the estimation methods under various settings. These values were carefully chosen to represent both typical and extreme cases, allowing us to assess the effectiveness and stability of the proposed methodology across a wide range of possible applications. The specific parameter values for \({\upsilon _1},{\upsilon _2}\) and \(\omega\) in the random variable generation process are chosen as follows:

In Table 2: \({\upsilon _1}=2,{\upsilon _2}=1.5,\omega =0.7\); In Table 3: \({\upsilon _1}=2,{\upsilon _2}=1.5,\omega =2\); In Table 4: \({\upsilon _1}=0.7,{\upsilon _2}=0.5,\)\(\omega =0.8\)and In Table 5: \({\upsilon _1}=0.7,{\upsilon _2}=2,\omega =0.8.\)

As mentioned in the preceding section, the hyperparameters for the priors \(\left( {{a_i},{b_i}} \right),\) where i ranges from 1 to 3, are found by applying the technique described by Dey et al. [56]. Two distinct schemes (Sc) are taken into consideration with respect to the censored samples during the removal procedure.

Sc 1: \({Q_i}=0;\,\,i=1,2,...,{m_j} - 1\)and \({Q_j}={n_j} - {m_j};\,\,j=1,...,5,y;\) while Sc 2: \({Q_{{m_j}}}={n_j} - {m_j};\,\,\)and \({Q_i}=0;\,\,i=2,...,{m_j}\)

For PFIF-CS, we choose h = 2 groups as the size for this mechanism. The value of the LLF scale parameter is fixed at 1 and 1.5. By computing the average bias and their accompanying mean squared errors (MSEs) and average of confidence lengths (ACLs), the parameter and the MC-SS reliability estimates for the IEPD based on PFIF-CS are assessed. The representation of ACLs in As-CI will be represented as LAs-CI with coverage probability (CP), whereas ACLs in HPD intervals will be LCCI.

The four tables provide a comprehensive analysis of parameter estimates for the MC-SS reliability model under the IEPD approach using PFIF-CS, focusing on different parameter combinations \({\upsilon _1},{\upsilon _2}\) and \(\omega\). All tables compare the MLE and the BE under different sample size settings. The metrics evaluated include bias, MSE, LAs-CI, LCCI, and CP, providing a detailed perspective on the efficiency, accuracy, and reliability of these estimators.

In tables, the effectiveness of the estimators is assessed under various scenarios, revealing differences in accuracy and reliability. For bias and MSE, lower values are ideal, as they indicate a more accurate estimation process. LAs-CI, LCCI reflects the uncertainty range, with shorter intervals typically being preferred, provided they still maintain adequate coverage. CP measures how often the interval captures the true value of the parameter, with values close to the desired confidence level (e.g., 95%) being favorable. By comparing these metrics for different parameter settings, the tables provide insights into the robustness of MLE and BE under varying conditions, offering valuable guidance for selecting appropriate estimation methods in reliability modeling.

Tables 2, 3 and 4 allow for the following important conclusions to be drawn:

  • Strong performance is demonstrated by the proposed MLEs and BEs for the IEPD’s unknown parameters in terms of bias, MSEs, and ACLs.

  • The suggested estimates in terms of bias, MSEs, and ACLs improve when sample size (i.e., the ratio of censored data to total observations) grows, which is consistent with predictions.

  • By including INP, the BEs that are obtained from gamma prior distributions perform better than the frequentist estimates.

  • The LLF (asymmetric) consistently outperforms the SELF (symmetric) when comparing symmetric and asymmetric loss functions.

  • The LLF is typically used by BEs to outperform other loss functions under comparison.

  • Components are individual elements of a system, each experiencing a unique stress level. The efficiency of the model grows with the number of components.

  • In comparison to conventional MLEs of\({\eta _{r,s}},\) BEs utilizing SELF and LLF perform better at estimating risks.

  • In the MC-SS model for the IEPD, the true reliability value is directly linked to the shape parameter \({\upsilon _1}\)of the strength distribution and the shape parameter of the stress distribution \({\upsilon _2}\). This relationship holds true for both (r, s) cases analyzed in the present work.

Table 2 The MC-SS reliability and parameter estimation under PFIF-CS: \({\upsilon _1}=2,{\upsilon _2}=1.5,\omega =0.7\)
Full size table
Table 3 The MC-SS reliability and parameter estimation under PFIF-CS: \({\upsilon _1}=2,{\upsilon _2}=1.5,\omega =2\)
Full size table
Table 4 The MC-SS reliability and parameter estimation under PFIF-CS: \({\upsilon _1}=0.7,{\upsilon _2}=0.5,\omega =0.8\)
Full size table
Table 5 The MC-SS reliability and parameter estimation under PFIF-CS: \({\upsilon _1}=0.7,{\upsilon _2}=2,\omega =0.8\)
Full size table

5 Data Analysis

We also investigate the suitability of the IEPD for analyzing the two datasets. This dataset was fitted using IEPD, and for comparative distributions, models as the Burr XII distribution (BXIID), the generalized inverted exponential distribution (GIED), and the exponentiated Pareto distribution (EPD) were also considered. The MLE and standard error (between brackets) for the unknown parameters of all competing models are shown in Tables 6 and 10. Additionally, these tables provide the Kolmogorov-Smirnov distance (KSD) statistic values along with their associated p-values (KSPV). In addition to Cramer-von Mises (CVM) and Anderson-Darling (AD) statistics have been considered and provided.

5.1 First Data Set

We examine a dataset that Maurya and Tripathi [42] have previously looked at. Failure times for a certain software model are included in this dataset. See “Link” for additional details. Variables \(\:{X}_{\left\{11\right\}},\:\dots\:,\:{X}_{\left\{15\right\}}\) represent failure times for units 120 to 124, while the variable \(\:{Y}_{1}\) indicates the failure time of the 125th unit. Similarly, \(\:{X}_{\left\{21\right\}},\:\dots\:,\:{X}_{\left\{25\right\}}\) denote failure times for units 126 to 130, and \(\:{Y}_{2}\) represents the failure time of the 131st unit. This pattern continues until the failure of unit 167. We are analyzing 8 systems, each with 5 components. The data set is as follows:

  • - For System 1, Y is 0.0037, while the values for \(\:{X}_{1},\:\:{X}_{2},\:{X}_{3},\:{X}_{4}\), and \(\:{X}_{5}\) are 0.0096, 0.0036, 0.0120, 1.0752, and 1.6488, respectively.

  • - For System 2, Y is 2.9736, and the predictor values are \(\:{X}_{1}\) = 5.4600, \(\:{X}_{2}\) = 1.1976, \(\:{X}_{3}\)= 0.0168, \(\:{X}_{4}\) = 0.1428, \(\:{X}_{5}\) = 1.5264.

  • - For System 3, Y is 1.4952, with the predictor values \(\:{X}_{1}\) = 4.7568, \(\:{X}_{2}\) = 0.0912, \(\:{X}_{3}\) = 0.0384, \(\:{X}_{4}\) = 0.3372, \(\:{X}_{5}\) = 1.5408.

  • - For System 4, Y is 0.5508, with \(\:{X}_{1}\)= 14.7636, \(\:{X}_{2}\)= 2.9400, \(\:{X}_{3}\)= 1.5336, \(\:{X}_{4}\)= 1.4160, \(\:{X}_{5}\)= 1.3236.

  • - For System 5, Y is 3.3420, and \(\:{X}_{1}\) = 1.1004, \(\:{X}_{2}\)= 0.0024, \(\:{X}_{3}\)= 5.93592, \(\:{X}_{4}\)= 1.8036, \(\:{X}_{5}\)= 14.4912.

  • - For System 6, Y is 1.9572, and the predictor values are \(\:{X}_{1}\)= 2.1240, \(\:{X}_{2}\)= 0.0096, \(\:{X}_{3}\)= 2.3472, \(\:{X}_{4}\)= 0.1296, \(\:{X}_{5}\)= 0.0012.

  • - For System 7, Y is 0.0108, and the predictor values are \(\:{X}_{1}\)= 7.9716, \(\:{X}_{2}\)= 7.8276, \(\:{X}_{3}\)= 0.0036, \(\:{X}_{4}\)= 0.0036, \(\:{X}_{5}\)= 0.0048.

  • - For System 8, Y is 0.0816, with \(\:{X}_{1}\)= 6.7380, \(\:{X}_{2}\)= 5.2896, \(\:{X}_{3}\)= 16.9416, \(\:{X}_{4}\)= 0.5292, \(\:{X}_{5}\)= 6.9864.

This dataset includes diverse values across systems for both the response and predictor variables, useful for evaluating multicomponent system performance.

Table 6 presents statistical inferences for stress and strength variables across four different models: IEPD, EPD, BXIID, and GIED. The parameter estimates and their standard errors (in parentheses) are provided for each of the variables (\(\:{x}_{1},\:{x}_{2},\:{x}_{3},\:{x}_{4},\:{x}_{5},\) and \(\:y\)). The goodness-of-fit measures include KSD and associated KSPV, CVM statistics, AD statistics. The IEPD generally provides better fits, as indicated by the lowest KSD and CVM values for many variables, suggesting a stronger alignment with empirical data. The KSPV values are also consistently higher for the IEPD, indicating significant results. Overall, the IEPD demonstrates the best performance across various metrics, suggesting it is the most suitable for fitting the stress and strength data compared to EPD, BXIID, and GIED.

Initially, we confirm that the IEPD model can be used to analyze the dataset represented in Figs. 3 and 4. Figure 3 shows the estimated CDF of the IEPD, BXIID, EPD, and GIED compared with the empirical CDF for each X-component. Figure 4 illustrates the estimated PDF of the IEPD along with histograms for each X-component. From these figures, we observe that the majority of the data fits well with the IEPD model.

Table 6 Statistical inference for each stress and strength variables for IEPD, EPD, BXIID, and GIED
Full size table
Fig. 3
figure 3

Calculated CDF for comparative models with each variable in the dataset

Full size image
Fig. 4
figure 4

Determined PDF for comparative models with each variable in the dataset

Full size image

Initially, it is assumed that \(\:(\:{X}_{1},\:\dots\:,\:{X}_{s}\sim\:\)IEPD (\(\omega,{\upsilon _1}\))) and (Y\(\:\sim\) IEPD (\(\omega,{\upsilon _2}))\). The MLE for the unknown parameters are: \(\hat {\omega }=\)0.4233, \({\hat {\upsilon }_1}=\) 0.6661 and \({\hat {\upsilon }_2}=\)0.8351, with a corresponding log-likelihood value of (L0 =−72.96877). Then, a second scenario assumes \(\:(\:{X}_{1},\:\dots\:,\:{X}_{s}\sim\) IEPD (\({\omega _1},{\upsilon _1}\))) and (Y\(\:\sim\) IEPD (\({\omega _2},{\upsilon _2}\))). The MLE for this case is \({\hat {\omega }_1}=\)0.4095, \({\hat {\upsilon }_1}=\) 0.6556, \({\hat {\omega }_2}=\)0.4966, and \({\hat {\upsilon }_2}=\)0.9191 with a log-likelihood value of (L1 =−72.90671). Subsequently, hypothesis tests are performed.

H0: \({\omega _1}={\omega _2}=\omega\) vs. H1: \({\omega _1} \ne {\omega _2}.\)

In this case, the value of the likelihood ratio test is calculated as \(\:-\)2(L0\(\:-\)L1) = 0.1241, with 1 degree of freedom. Additionally, the P-value from the Chi-squared test is 0.2753, indicating that the null hypothesis cannot be rejected. Therefore, H0 is considered valid in this scenario.

Table 7 Parameter and different reliability estimates of an MC-SS by ML estimation for comparative models under complete sample: data I
Full size table

Table 7 shows estimated values of the parameter and reliability for an MC-SS model using ML estimation for the IEPD, BXIID, EPD, and GIED based on Data I. The IEPD exhibits moderate parameter values and consistently provides the highest reliability estimates, particularly for \({\eta _{2,5}}\)and \({\eta _{3,5}}\), suggesting that it offers better performance compared to the other models. The BXIID and GIED also demonstrate relatively good reliability, while the EPD has the lowest reliability estimates across all metrics. Overall, the IEPD appears to be the most effective for modeling this MC-SS system.

Fig. 5
figure 5

Likelihood profile of parameters of IEPD for complete sample: Data I

Full size image

Figure 5 presents the likelihood profiles for the parameters \(\omega\), \({\upsilon _1}\), and \({\upsilon _2}\)of the IEPD based on the complete sample from Data I. Each plot shows the log-LF for one of the parameters, while keeping the others fixed, indicating the point where the likelihood is maximized. The blue dots mark the MLEs, representing the optimal parameter values that maximize the likelihood. The curves demonstrate a clear peak for each parameter, confirming well-defined MLEs and suggesting that the model parameters are identifiable and effectively estimated.

Table 8 Bayesian estimation of parameter and different MC-SS reliability for IEPD and EPD: complete sample
Full size table

Table 8 provides parameter and reliability estimate values for an MC-SS model using Bayesian estimation under N-INP for the IEPD and EPD. The IEPD shows consistently higher reliability values than the EPD, indicating that the IEPD may be more suitable for this data.

Table 9 presents the parameter and reliability estimates measures for an MC-SS model based on PFIF-CS using the IEPD, under two schemes, comparing MLE and BE with N-INP. The BEs generally yield higher parameter values and reliability estimates \({\eta _{2,5}}\), \({\eta _{3,5}}\), \({\eta _{4,5}}\), and \({\eta _{5,5}}\) compared to the MLE suggesting a more optimistic assessment of system reliability. Sc 1 shows slightly higher reliability measures compared to Sc 2, emphasizing that Bayesian approaches provide improved reliability estimates for the model compared to MLE.

Table 9 MC-SS Model: parameter and reliability based on PFIF-CS for first data
Full size table

Figure 6 displays trace and density plots for all parameters within the MCMC trace under Bayesian estimation. The posterior density of each parameter is shown, indicating a symmetric normal distribution that aligns with the expected distribution. These figures confirm the convergence of the MCMC results.

Fig. 6
figure 6

MCMC forms of MC-SS model parameters: Data I

Full size image

5.2 Second Data Set

In this example, one of our main objectives is to construct a scenario concerning excessive drought. The claim is that if, in at least one of the next five years, the water capacity of a reservoir in a particular region during May exceeds the water levels recorded in November of the previous year for studies, it can be concluded that no excessive drought will follow. To support this analysis, we examine monthly data from the Shasta Reservoir in California (refer to Kızılaslan and Nadar [58] and Kohansal [38]). This dataset can be accessed at “http://cdec.water.ca.gov/cgi-progs/queryMonthly?SHA.” The data provide insights into the amount of water discharged into rivers, the volume of water stored in reservoirs, precipitation accumulation, and snowpack water content, all from a flood management perspective. Typically, the reservoir’s maximum water level is observed in May, while the minimum is in November. Previous research by Kızılaslan and Nadar [58] and Nadar and Kızılaslan [59] explored this data in different contexts.

Our objective is to make inferences in order to take precautions against extreme drought situations under PFIF-CS. In the case of a complete sample with \(\:(\:r\:=\:5\:)\) and \(\:s\:=\:\text{2,3},\text{4,5}\), \(\:{Y}_{1}\) represents the capacity for November 1970, and \(\:{X}_{\left\{11\right\}},\:\dots\:,\:{X}_{\left\{15\right\}}\) represent the capacities for September from 1971 to 1975. Similarly, \(\:{Y}_{2}\) represents the capacity for November 1976, and \(\:{X}_{\left\{21\right\}},\:\dots\:,\:{X}_{\left\{25\right\}}\) represent the capacities for September from 1977 to 1981. This process continues until 2017, yielding \(\:N\:=\:8\) data points for \(\:Y\). For computational ease, we scaled each data point dividing by 4,552,000, which is the total capacity of the Shasta reservoir. The transformed data sets are provided below:

  • - First set: 0.719442, 0.138533, 0.794552, 0.460452, 0.689022, 0.483226, 0.412817, 0.418714.

  • - Second set: 0.717597, 0.753054, 0.711797, 0.359703, 0.678561, 0.561995, 0.304148, 0.254192.

  • - Third set: 0.728603, 0.690092, 0.434490, 0.294343, 0.507104, 0.694063, 0.389707, 0.352043.

  • - Fourth set: 0.803669, 0.729504, 0.705470, 0.369772, 0.755947, 0.479537, 0.729082, 0.617617.

  • - Fifth set: 0.784161, 0.544859, 0.463141, 0.681406, 0.730997, 0.666704, 0.733984, 0.742939.

  • - Y set: 0.767728, 0.343146, 0.724319, 0.400454, 0.443879, 0.633937, 0.686972, 0.563329.

Table 10 Statistical inference for each stress and strength variable for IEPD, EPD, BXIID, and GIED: data II
Full size table

Table 10 presents the parameter estimates of \(\omega\) and\(\upsilon\) for the IEPD, EPD, BXIID, and GIED across five variables (x1, x2, x3, x4, x5) and one dependent variable (y). Lower values of the KSD statistic, CVM, and AD with higher KSPV suggest a better fit. For results of Table 10: IEPD, EPD, and GIED consistently show lower KSD statistics and reasonable KSPV across the variables compared to other distributions, indicating that it fits the data well. The Y variable has changed results of this conclusion, as the best distribution fitting of the Y variable is IEPD.

Fig. 7
figure 7

CDF plots of four different models: IEPD, EPD, BXIID, and GIED for the variables \(\:{x}_{1},\:{x}_{2},\:{x}_{3},\:{x}_{4},\:{x}_{5}\) and y: Data II

Full size image
Fig. 8
figure 8

PDF plots of four different models: IEPD, EPD, BXIID, and GIED for the variables \(\:{x}_{1},\:{x}_{2},\:{x}_{3},\:{x}_{4},\:{x}_{5}\) and y: Data II

Full size image

By Figs. 7 and 8, compare four models (IEPD, EPD, BXIID, and GIED) against the empirical data for six variables.

Initially, it is assumed that \(\:(\:{X}_{1},\:\dots\:,\:{X}_{s}\sim\:\)IEPD (\(\omega,{\upsilon _1}\))) and (Y\(\:\sim\) IEPD (\(\omega,{\upsilon _2}\))). The MLE for the unknown parameters are: \(\hat {\omega }=\)6.1281, \({\hat {\upsilon }_1}=\) 333.4885 and \({\hat {\upsilon }_2}=\)373.4789, with a corresponding log-likelihood value of (L0 =\(\:-\)15.42049). Then, a second scenario assumes \(\:(\:{X}_{1},\:\dots\:,\:{X}_{s}\sim\) IEPD (\({\omega _1},{\upsilon _1}\))) and (Y\(\:\sim\) IEPD (\({\omega _2},{\upsilon _2}\))). The MLE for this case are \({\hat {\omega }_1}=\)5.9591, \({\hat {\upsilon }_1}=\) 286.6811, \({\hat {\omega }_2}=\)6.3495, and \({\hat {\upsilon }_2}=\)445.8943 with a log-likelihood value of (L1 =\(\:-\)15.48556). Subsequently, hypothesis tests are performed.

H0: \({\omega _1}={\omega _2}=\omega\) vs. H1: \({\omega _1} \ne {\omega _2}.\)

In this case, the value of the likelihood ratio test is calculated as \(\:-\)2(L0\(\:-\)L1) = 0.1301306 with one degree of freedom. Additionally, the P-value from the Chi-squared test is 0.2817, indicating that the null hypothesis cannot be rejected. Therefore, H0 is considered valid in this scenario.

Table 11 Parameter and different reliability estimates of an MC-SS by MLE for comparative models under a complete sample: data II
Full size table

Table 11 presents estimates of parameters and reliability based on MLE for four different models: IEPD, BXIID, EPD, and GIED, under a complete sample for Data II. The parameters \(\omega,\)\({\upsilon _1}\), and \({\upsilon _2}\) represent the model estimates, with each model showing significant variations in their values. The reliability measures \({\eta _{2,5}},\)\({\eta _{3,5}}\), \({\eta _{4,5}}\), and \({\eta _{5,5}}\) indicate the probability that the system’s strength exceeds the applied stress, with IEPD and BXIID demonstrating higher reliability across the columns, especially for \({\eta _{2,5}}\), and \({\eta _{5,5}}\). GIED and EPD show slightly lower reliability estimates, with EPD having the lowest reliability values. Overall, the results of this table highlight the effectiveness of the IEPD in achieving relatively higher reliability measures compared to the other three models. Coupled with the findings from the previous Table 11, it becomes evident that the IEPD plays a significant role in explaining the phenomenon and yielding the highest reliability values. This suggests that the IEPD is particularly well-suited for analyzing and interpreting the data in this context.

Figure 9 shows the log-LF plotted against the parameters \(\omega\), \({\upsilon _1}\), and \({\upsilon _2}\) with the blue dot indicating the MLEs. The curves peak at the MLE values, highlighting the optimal parameter estimates for the model.

Fig. 9
figure 9

Likelihood profile of parameters of IEPD for complete sample: Data II

Full size image

Figure 10 presents trace and density plots for BEs of all parameters within the MCMC analysis. The symmetric normal distribution shown in the posterior density plots for each parameter supports our earlier hypotheses and indicates that the MCMC technique has successfully converged. This convergence is further supported by the trace plots, which display the path of the MCMC chain across iterations and indicate that the sampler has sufficiently explored the parameter space because no patterns or trends are visible.

Fig. 10
figure 10

MCMC forms of MC-SS model parameters: Data II

Full size image

Table 12 presents MC-SS reliability and parameter estimates using Bayesian estimation under N-NIP for the IEPD in the case of a complete sample. The table shows that the BE results in higher parameter estimates for \(\omega\), \({\upsilon _1}\), and \({\upsilon _2}\).

Table 12 Bayesian estimation and different reliability of an MC-SS for IEPD: Data II under complete sample
Full size table

Additionally, BEs of \({\eta _{2,5}},\)\({\eta _{3,5}},\)\({\eta _{4,5}},\) and \({\eta _{5,5}}\) are generally higher, indicating a more optimistic assessment of reliability with this approach.

Table 13 MC-SS Model: parameter and reliability based on PFIF-CS for second data
Full size table

Table 13 compares the MLE and BE using N-INP for an MC-SS reliability model across two schemes. In both schemes, the BE consistently produces higher parameter estimates and reliability values \({\eta _{2,5}},\)\({\eta _{3,5}},\)\({\eta _{4,5}},\) and \({\eta _{5,5}}\)compared to MLE. The results suggest that the BE yields the most optimistic estimates, while MLE provides the most conservative reliability measures.

6 Concluding Remarks

The primary objective of this study is to investigate reliability inference for the MC-SS model under PFIFC-S. We assume that both stress and strength are independent random variables that follow the IEPD with a shared second shape parameter. We employ both Bayesian and non-Bayesian procedures to obtain parameter and MC-SS reliability estimators. Additionally, we construct credible intervals for the highest posterior density and asymptotic confidence intervals. Using the MCMC approach, we obtain Bayesian estimates of the MC-SS reliability under both symmetric (squared error) and asymmetric (linear exponential) loss functions. To evaluate the effectiveness of our proposed methodology, we conduct comprehensive simulation analyses. We compare IEPD’s performance to other models using goodness-of-fit metrics along with graphical representation in order to assess the model’s suitability for analyzing two datasets. We evaluate the competing models’ MC-SS reliability and estimate unknown parameters using Bayesian and non- Bayesian techniques. The analysis demonstrates that the IEPD consistently generates the most accurate estimates and performs better than the other models. One of the study’s limitations is that it solely focuses on the premise that strength and stress are independent random variables. An additional restriction is taken into account only for the two loss functions (SELF, LLF) when calculating the Bayesian estimators for utilizing the MCMC method. Further studies can be considered in the following scenarios:

  • Develop the MC-SS model that accounts for the dependence between stress and strength variables [60, 61].

  • Explore alternative Bayesian estimation methods, such as variational Bayes.

  • Investigating the impact of other loss functions on the Bayesian estimates of reliability.