Model Selection for Reliability Estimation in Series Systems

The statistical treatment of masked system failure data has a substantial history. Usher and Hodgson [12] introduced maximum likelihood methods for component reliability estimation from masked system life-test data, establishing the foundational framework for this field. Lin, Usher, and Guess [13] extended this work to derive exact maximum likelihood estimators, while Lin, Usher, and Guess [6] developed Bayesian approaches for component reliability estimation from masked data.

For Weibull-distributed components specifically, Usher [14] addressed component reliability prediction in the presence of masked data. Guess and Usher [3] proposed an iterative approach for estimating component reliability that handles the computational challenges of masked data likelihood functions.

Sarhan [7] examined reliability estimation from masked system life data under various distributional assumptions, later extending this work to linear failure rate models [8]. Tan [9, 10] contributed methods for exponential component reliability estimation from masked binomial system testing data and uncertain life data in series and parallel systems. More recently, Guo, Niu, and Szidarovszky [5] studied estimating component reliabilities from incomplete system failure data, providing the baseline system configuration used in our simulation studies.

Building on this literature, Towell [11] developed a comprehensive likelihood model that incorporates both right-censoring and candidate sets for series systems with Weibull components, along with extensive simulation studies validating the maximum likelihood approach. The present paper extends that work by investigating model selection—specifically, when a reduced model assuming homogeneous component shapes is appropriate.

This paper addresses a fundamental question in reliability engineering: when is a simpler model adequate? The key contributions are:

1. 1.

   A striking robustness result: For well-designed series systems, the reduced homogeneous-shape model cannot be rejected even with sample sizes approaching 30,000—far larger than typically available in practice. This provides strong justification for using the simpler model.

2. 2.

   Sharp sensitivity boundaries: We quantify exactly how much component heterogeneity is needed before the likelihood ratio test rejects the reduced model, across a range of sample sizes and shape parameter values.

3. 3.

   Practical model selection guidance: Our results translate directly into actionable recommendations for practitioners facing the bias-variance tradeoff in reliability assessment.

The remainder of this paper is organized as follows. Section 2 presents mathematical preliminaries including formal notation, series system definitions, Weibull distribution properties, and foundational theorems. Section 3 summarizes the likelihood model for masked and censored data. Section 4 presents simulation studies assessing estimator sensitivity to system design variations. Section 5 introduces the reduced homogeneous shape model and evaluates its appropriateness using likelihood ratio tests. Section 6 concludes with practical guidance on model selection.

Consider a series system composed of $m$ components, where the system fails when any single component fails. We observe $n$ independent and identically distributed system lifetimes. For the $i$-th system ($i=1,\mathrm{\dots },n$), let $T_{ij}$ denote the lifetime of component $j$ ($j=1,\mathrm{\dots },m$). The system lifetime is given by

$$T_i=\min \{T_{i1},T_{i2},\mathrm{\dots },T_{im}\}.$$

(1)

We denote the complete parameter vector by $𝜽=(k_1,{\lambda }_1,k_2,{\lambda }_2,\mathrm{\dots },k_m,{\lambda }_m)$, where $k_j>0$ is the shape parameter and ${\lambda }_j>0$ is the scale parameter for component $j$. Bold symbols represent vectors throughout this paper.

Each component lifetime follows a two-parameter Weibull distribution with probability density function

$$f_j(t;{\lambda }_j,k_j)=\frac{k_j}{{\lambda }_j}{\left(\frac{t}{{\lambda }_j}\right)}^{k_j-1}\exp \left\{-{\left(\frac{t}{{\lambda }_j}\right)}^{k_j}\right\},t>0,$$

(2)

reliability function

$$R_j(t;{\lambda }_j,k_j)=\exp \left\{-{\left(\frac{t}{{\lambda }_j}\right)}^{k_j}\right\},$$

(3)

and hazard function

$$h_j(t;{\lambda }_j,k_j)=\frac{k_j}{{\lambda }_j}{\left(\frac{t}{{\lambda }_j}\right)}^{k_j-1}.$$

(4)

The mean time to failure (MTTF) for a Weibull-distributed component is

$${\text{MTTF}}_j={\lambda }_j\mathrm{\Gamma }\left(1+\frac{1}{k_j}\right),$$

(5)

where $\mathrm{\Gamma }(\cdot )$ is the gamma function. The shape parameter $k_j$ characterizes the failure mode:

• •

  If , the hazard function decreases with time, indicating infant mortality or early-life failures.

• •

  If $k_j=1$, the hazard function is constant, corresponding to an exponential distribution with memoryless failures.

• •

  If $k_j>1$, the hazard function increases with time, indicating wear-out or aging failures.

The observed data for each system $i$ consists of:

• •

  Observed system lifetime $t_i$: The time at which system $i$ fails or is censored.

• •

  Censoring indicator ${\delta }_i$: Equals 1 if system $i$ failed and 0 if right-censored.

• •

  Candidate set $C_i\subseteq \{1,2,\mathrm{\dots },m\}$: For failed systems (${\delta }_i=1$), the set of components that could have caused the failure. For censored systems, $C_i$ is undefined.

This data structure is termed masked data because the true component cause of failure is not directly observed but is known to belong to the candidate set $C_i$. The masking mechanism assumes that:

1. 1.

   The candidate set always contains the true failed component.

2. 2.

   Given the system failure time and the true failed component, the masking mechanism is non-informative, meaning the process generating candidate sets is independent of the parameter vector $𝜽$. For a detailed treatment of these conditions, see [11].

A well-designed series system is characterized by components having similar but not necessarily identical failure characteristics. Operationally, we define a well-designed system as one where:

1. 1.

   Component MTTFs are of similar magnitude (within a factor of 2-3).

2. 2.

   Component shape parameters are reasonably aligned (within approximately 20-30% of each other).

3. 3.

   No single component dominates as a weak point (i.e., component failure probabilities are relatively balanced).

This concept is important for assessing the appropriateness of reduced models that assume parameter homogeneity.

The following properties characterize the lifetime distribution of series systems with independent component lifetimes.

For series systems with Weibull components, these properties yield specific forms presented in the following subsection.

Applying Properties 2.2, 2.3, and 2.1 to series systems with Weibull-distributed components yields specific analytical forms for the system lifetime distribution.

The lifetime of the series system composed of $m$ Weibull components has a reliability function given by

$$R(t;𝜽)=\exp \left\{-\sum _{j=1}^m{\left(\frac{t}{{\lambda }_j}\right)}^{k_j}\right\}.$$

(6)

The Weibull series system’s hazard function is given by

$$h(t;𝜽)=\sum _{j=1}^m\frac{k_j}{{\lambda }_j}{\left(\frac{t}{{\lambda }_j}\right)}^{k_j-1},$$

(7)

whose proof follows from Property 2.3.

The pdf of the series system is given by

$$f(t;𝜽)=\left\{\sum _{j=1}^m\frac{k_j}{{\lambda }_j}{\left(\frac{t}{{\lambda }_j}\right)}^{k_j-1}\right\}\exp \left\{-\sum _{j=1}^m{\left(\frac{t}{{\lambda }_j}\right)}^{k_j}\right\}.$$

(8)

When components have heterogeneous shape parameters, the series system hazard function can exhibit complex behavior, including both infant mortality (initial decrease) and aging (eventual increase) phases. This bathtub-shaped hazard is commonly observed in engineered systems where early failures due to defects give way to a period of stable operation, eventually followed by wear-out failures.

Throughout this paper, we use a baseline 5-component series system configuration for simulation studies. The component parameters are specified in Table 1.

This baseline system represents a well-designed series system configuration. All components have similar MTTFs ranging from approximately 799 to 913 time units, ensuring no single component dominates as a weak point. The shape parameters are tightly clustered between 1.13 and 1.26, all indicating slight aging behavior ($k_j>1$) with similar failure characteristics. Component 3, with shape parameter $k_3=1.1308$, has the smallest shape value and serves as a reference point for sensitivity analyses. This configuration is ideal for assessing the appropriateness of the reduced homogeneous-shape model, as the components already exhibit substantial similarity in their failure modes.

In series systems, a critical quantity for understanding system behavior and estimator performance is the probability that a particular component causes system failure. Let $K_i$ denote the index of the component that causes the $i$-th system to fail. The probability that component $j$ is the cause of failure is given by:

$$P_j=\mathrm{Pr}\{K_i=j\}={\int }_0^{\mathrm{\infty }}{f}_{T_i,K_i}(t,j;𝜽)𝑑t,$$

(9)

$$P_j={\int }_0^{\mathrm{\infty }}{f}_j(t;{\theta }_j)R_{\setminus j}(t;{𝜽}_{\setminus j})𝑑t=E_{𝜽}\left\{\frac{h_j(T_i;{\theta }_j)}{h(T_i;𝜽)}\right\},$$

(10)

For Weibull components in series, the component failure probability depends on both shape and scale parameters in a complex, non-linear manner. A key insight is that MTTF alone is insufficient for determining failure probabilities in series systems with heterogeneous shape parameters.

For the $i$-th system ($i=1,\mathrm{\dots },n$), the observed data consists of:

• •

  System failure or censoring time $t_i$

• •

  Censoring indicator ${\delta }_i\in \{0,1\}$

• •

  Candidate set $C_i\subseteq \{1,2,\mathrm{\dots },m\}$ (for failed systems only)

Let $K_i\in \{1,\mathrm{\dots },m\}$ denote the (unobserved) component cause of failure for system $i$. The likelihood contribution for a single observation depends on whether the system failed or was censored.

For a censored observation (${\delta }_i=0$), the likelihood contribution is simply the system reliability function evaluated at the censoring time:

$$L_i(𝜽|t_i,{\delta }_i=0)=R(t_i;𝜽)=\prod _{j=1}^m\exp \left\{-{\left(\frac{t_i}{{\lambda }_j}\right)}^{k_j}\right\}.$$

(11)

For a failed system (${\delta }_i=1$) with candidate set $C_i$, the likelihood contribution accounts for the fact that the failed component is known to be in $C_i$ but is otherwise unknown:

$$L_i(𝜽|t_i,{\delta }_i=1,C_i)=\sum _{j\in C_i}\mathrm{Pr}\{K_i=j|t_i,𝜽\}\cdot f(t_i;𝜽),$$

(12)

where the conditional probability that component $j$ failed given system failure time $t_i$ is

$$\mathrm{Pr}\{K_i=j|t_i,𝜽\}=\frac{h_j(t_i;{\lambda }_j,k_j)}{h(t_i;𝜽)}=\frac{h_j(t_i;{\lambda }_j,k_j)}{{\sum }_{\mathrm{\ell }=1}^m{h}_{\mathrm{\ell }}(t_i;{\lambda }_{\mathrm{\ell }},k_{\mathrm{\ell }})}.$$

(13)

The complete log-likelihood for the sample of $n$ systems is

$$\mathrm{\ell }(𝜽|D)=\sum _{i=1}^n\left[(1-{\delta }_i)\log R(t_i;𝜽)+{\delta }_i\log \left(\sum _{j\in C_i}\mathrm{Pr}\{K_i=j|t_i,𝜽\}\cdot f(t_i;𝜽)\right)\right].$$

(14)

The likelihood model relies on the following assumptions:

1. 1.

   Independence: System lifetimes are independent and identically distributed.

2. 2.

   Series structure: The system fails when the first component fails.

3. 3.

   Weibull components: Each component lifetime follows a two-parameter Weibull distribution.

4. 4.

   Candidate set validity: For failed systems, the candidate set $C_i$ always contains the true failed component, i.e., $K_i\in C_i$.

5. 5.

   Non-informative masking: Given the system failure time $t_i$ and the true failed component $K_i$, the masking mechanism that determines $C_i$ is non-informative about the parameters $𝜽$.

6. 6.

   Identifiability: We assume standard regularity conditions ensuring that the parameter vector $𝜽$ is identifiable from the observed data. For series systems with masked data, identifiability requires sufficient variation in component failure times and candidate sets; see [11] for a detailed treatment.

Maximum likelihood estimates are obtained by maximizing the log-likelihood function $\mathrm{\ell }(𝜽|D)$ with respect to $𝜽$. Due to the complexity of the likelihood surface with $2m$ parameters, numerical optimization is required. The optimization is performed using the L-BFGS-B algorithm [1], which handles box constraints on the parameters (all shape and scale parameters must be positive).

Previous simulation studies [11] demonstrated that maximum likelihood estimation produces accurate results despite small samples and significant masking and censoring. However, shape parameters exhibit greater variability and are more challenging to estimate precisely than scale parameters, particularly when candidate sets are large or when certain components rarely fail.

For each simulation scenario, we employ the following methodology:

• •

  Number of replications: 1000 Monte Carlo replications per parameter configuration

• •

  Sample size: $n=100$ systems per replication (unless otherwise specified)

• •

  Masking probability: $p=0.215$ (moderate masking, unless otherwise specified)

• •

  Censoring mechanism: Right-censoring at the $q=0.825$ quantile of the system lifetime distribution (moderate censoring, unless otherwise specified)

• •

  Candidate set generation: For each failed system, components are independently included in the candidate set with probability $p$, ensuring the true failed component is always included

• •

  Optimization: L-BFGS-B algorithm [1] with box constraints requiring all parameters to be positive

• •

  Confidence intervals: Bias-corrected and accelerated (BCa) bootstrap confidence intervals [2] with 2000 bootstrap samples per replication, targeting 95% nominal coverage

• •

  Performance metrics: Bias, dispersion (interquartile range of point estimates), coverage probability, and confidence interval width

The term moderate refers to levels that are substantial enough to pose estimation challenges but not so extreme as to make inference infeasible. Specifically, a masking probability of 0.215 results in candidate sets containing approximately 2-3 components on average, and a censoring quantile of 0.825 censors approximately 17.5% of observations.

By Equation 5, we see that MTTF_j is proportional to the scale parameter ${\lambda }_j$, which means when we decrease the scale parameter of a component, we proportionally decrease the MTTF. In this scenario, we start with the well-designed series system described in Table 1, and we will manipulate the MTTF of component 3, MTTF₃, by changing its scale parameter, ${\lambda }_3$, and observing the effect this has on the MLE. Since the other components had a similar MTTF, we will arbitrarily choose component 1 to represent the other components. The bottom plot shows the coverage probabilities for all parameters.

In Figure 1, we show the effect of changing the scale parameter of component $3$, ${\lambda }_3$, but map ${\lambda }_3$ to MTTF₃ to make it more intuitive to reason about. We vary the MTTF of component 3 from $300$ to $1500$ and the other components have their MTTFs fixed at around $900$, as shown in Table 1. We fix the masking probability to $p=0.215$ (moderate masking), the right-censoring quantile to $q=0.825$ (moderate censoring), and the sample size to $n=100$ (moderate sample size).

The shape parameter determines the failure characteristics. We vary the shape parameter of component 3 from $0.1$ to $3.5$ and observe the effect it has on the MLE. When , this indicates infant mortality, and when $k_3>1$, this indicates wear-out failures.

We analyze the effect of component 3’s shape parameter on the MLE and the bootstrapped confidence intervals for the shape and scale parameters of components 1 and 3 (the component we are varying). First, we look at the effect on the scale parameter.

The homogeneous shape model (also called the reduced model) assumes that all components share a common shape parameter $k$ while retaining individual scale parameters ${\lambda }_j$. The parameter vector for the reduced model is

$${𝜽}_R=(k,{\lambda }_1,{\lambda }_2,\mathrm{\dots },{\lambda }_m),$$

(15)

reducing the number of parameters from $2m$ in the full model to $m+1$ in the reduced model.

Under this assumption, each component has a Weibull lifetime with

$$T_{ij}\sim \text{Weibull}(k,{\lambda }_j),j=1,\mathrm{\dots },m.$$

(16)

A key property of the homogeneous shape model is that the series system lifetime distribution is itself Weibull. Specifically, for the system lifetime $T_i=\min \{T_{i1},\mathrm{\dots },T_{im}\}$, we have

$$T_i\sim \text{Weibull}(k,{\lambda }_s),$$

(17)

where the system scale parameter is given by

$${\lambda }_s={\left(\sum _{j=1}^m{\lambda }_j^{-k}\right)}^{-1/k}.$$

(18)

This Weibull property of the system lifetime provides several advantages:

1. 1.

   Analytical tractability: System reliability metrics (MTTF, reliability function, hazard function) have closed-form expressions.

2. 2.

   Interpretability: The system exhibits a single failure mode characterized by the common shape parameter $k$.

3. 3.

   Reduced variance: Fewer parameters to estimate results in lower estimator variance, particularly for the shape parameter.

4. 4.

   Computational efficiency: Optimization is faster with $m+1$ parameters instead of $2m$ parameters.

However, the reduced model is appropriate only when components genuinely have similar failure characteristics. When component shape parameters differ substantially, forcing homogeneity can lead to poor model fit and biased parameter estimates.

In order to determine if a reduced model (e.g., Weibull series system in which all of the shape parameters are homogeneous) is appropriate, a hypothesis test may be conducted to determine if there is statistically significant evidence against the null hypothesis $H_0$ that all of the shape parameters are homogeneous. If the test fails to reject $H_0$, the reduced model is supported; if it rejects $H_0$, the full model with heterogeneous shapes is preferred.

The likelihood function of the reduced model is related to the likelihood function of the full model. We denote the full model likelihood function by $L_F$ and the reduced model likelihood by $L_R$. The reduced model is obtained by setting the shape parameter of each component to be the same, i.e., $k_1=\mathrm{\cdots }=k_m=k$. Thus, the reduced model likelihood function is given by

$$L_R(k,{\lambda }_1,{\lambda }_2,\mathrm{\cdots },{\lambda }_m|D)=L_F(k,{\lambda }_1,k,{\lambda }_2,\mathrm{\dots },k,{\lambda }_m|D),$$

The same may be done for the score and hessian of the log-likelihood functions.

Given that we employ a well-defined likelihood model, the likelihood ratio test (LRT) is a good choice. The LRT statistic is given by

$$\mathrm{\Lambda }=-2(\log L_R({\widehat{\theta }}_R|D)-\log L_F(\widehat{\theta }|D))$$

where $L_R$ is the likelihood of the reduced (null) model evaluated at its MLE ${\widehat{\theta }}_R$ given a random sample $D$ of masked data and $L_F$ is the likelihood of the full model evaluated at its MLE $\widehat{\theta }$ given the same set of data $D$. Under the null model, the LRT statistic is asymptotically distributed chi-squared with $m-1$ degrees of freedom, where $m$ is the number of components in the series system,

$$\mathrm{\Lambda }\sim {\chi }_{m-1}^2.$$

If the LRT statistic is greater than the critical value of the chi-squared distribution with $m-1$ degrees of freedom, ${\chi }_{m-1,1-\alpha }^2$, where $\alpha$ denotes the significance level, then we find the data to be incompatible with the null hypothesis $H_0$.

To assess how far a system deviates from the homogeneous-shape assumption, we introduce a divergence metric based on the coefficient of variation (CV) of the shape parameters:

$${\text{CV}}_k=\frac{\text{sd}(k_1,k_2,\mathrm{\dots },k_m)}{\text{mean}(k_1,k_2,\mathrm{\dots },k_m)}.$$

(19)

A system with ${\text{CV}}_k=0$ has perfectly homogeneous shape parameters, while larger values indicate greater heterogeneity. An alternative metric is the max/min ratio, ${\max }_j(k_j)/{\min }_j(k_j)$, which equals 1 for homogeneous systems.

Table 2 shows how varying $k_3$ in our baseline system translates to different divergence levels. The well-designed baseline ($k_3=1.1308$) has CV $\approx 4\%$ and max/min $\approx 1.11$, representing very mild heterogeneity.

We aim to assess the appropriateness of the reduced model under varying divergence levels (controlled by $k_3$) and sample sizes. We employ a simulation study using the likelihood ratio test (LRT) for this purpose, where the null hypothesis, $H_0$, assumes homogeneous shape parameters.

We take the well-designed series system described in Table 1, and manipulate the shape parameter of the third component ($k_3$) to cause the components to have different failure characteristics. As shown in Table 2, $k_3=1.1308$ corresponds to a well-designed series system with CV $\approx 4\%$, while extreme values like $k_3=0.25$ or $k_3=3.0$ produce systems with CV $>40\%$. We also vary the sample size $n$ to assess the impact of sample size on the appropriateness of the reduced model.

The previous analysis focused on p-values from the LRT under varying divergence levels. To provide a more complete understanding of LRT behavior, we conducted additional simulations examining: (1) Type I error control under perfect homogeneity, (2) power curves across divergence levels, and (3) the effects of masking probability, censoring, and system complexity on test performance.

The extended simulation study reveals that the LRT has correct Type I error control and predictable power characteristics. The power of the LRT for well-designed series systems (CV ) is remarkably low, requiring many thousands of observations before the test has sufficient power to reject the null hypothesis. This is not necessarily problematic—it indicates that the reduced model genuinely fits well-designed systems.

Our analysis suggests practical guidelines based on divergence levels:

1. 1.

   Low divergence (CV ): Use the reduced model confidently. The LRT will not reject it even with very large samples, and it provides the benefits of simplicity, interpretability, and reduced variance.

2. 2.

   Moderate divergence (CV 10–20%): The choice depends on sample size. For , the reduced model is unlikely to be rejected and may be preferred. For larger samples, consider the full model.

3. 3.

   High divergence (CV $>25\%$): Use the full model. Even modest sample sizes will reject the reduced model, indicating genuine heterogeneity that should be captured.

For systems believed to be well-designed, employing the reduced model is supported both statistically and practically due to its simplicity, reduced estimator variability, and analytical tractability. When prior information about component homogeneity is unavailable, engineers should estimate the shape CV from initial data to guide model selection.