1.
Introduction
With the development and progress of manufacturing technology, modern products are designed with complex structures and have high reliability. However, for some high reliablity products, it is hard to obtain their failure data through traditional life tests within a short period of time. However, in many cases, degradation measurements can provide valuable information related to the failure mechanism of the product. Therefore, the product's reliability can be inferred and estimated through the degradation data of quality characteristics obtained [1].
In recent years, various kinds of degradation models have emerged and been studied. Some of these models are stochastic process models, mixed-effect models [2,3], and so on. Typical degradation models include the general degradation path models and stochastic process models. Meeker et al.[4] used a nonlinear regression model with mixed-effects to analyze constant-stress accelerated degradation test (CSADT) data. Shi and Meeker [5] discussed the accelerated destructive degradation test planning of a nonlinear regression model through the Bayesian method. The stochastic process models include the Wiener process model [6,7,8,9,10,11,12], the Gamma process model [13,14,15,16,17], and the inverse Gaussian (IG) process model [18,19,20,21,22]. Although Ye and Xie [23] have made a comprehensive study on degradation analysis of products with single quality characteristics (QC), reliability analysis of complex systems with two or more competing failure modes (e.g. sudden failure, degradation failure) is still a challenge.
The modeling and statistical analysis of competing risk data has increasingly become a hot issue in the field of reliability, and there are many literatures on the statistical analysis of competing risk data (or competing risk model), such as Nassr et al. [24], Ramadan et al. [25], Mohamed et al. [26] and Mohamed et al. [27]. Huang et al. [28] studied the optimal maintenance scheme of multi-dependent competitive degradation and shock processes. Fan et al. [29] used degradation-shock dependence to model the dependent competitive failure process and used Monte Carlo techniques to calculate the reliability of the system. Xu et al. [30] modeled competing failure with the bivariate Wiener degradation process. Wang et al. [31] proposed two semiparametric additive mean models for clustered panel count data, and estimated the regression parameters of interest by constructing the estimation equations. Mutairi et al. [32] studied the inverse Weibull model based on jointly type-II hybrid censoring samples through the Bayesian or non-Bayesian methods. Bhat et al. [33] discussed the properties and Bayesian estimation of the odd lindley power rayleigh distribution.
A motivating example of this study is provided by Huang and Askin [34]. Units in the system may fail when the solder/pad interface breaks due to fatigue [35], or when the electrical/optical signal drops to unacceptable levels due to aging degradation [36]. In this example, an electronic device failed caused by two independent failure elements: the light intensity degradation (soft failure), which is considered a degradation phenomenon, because at some common inspection times to observe and measure the light intensity of the device, the solder/bond pad interface breaks, which is regarded as a sudden failure (hard failure). The original data given in Tables 1 and 2 was measured under the same conditions. The degradation data is the ratio of the current brightness to the startup brightness. When this ratio is reduced by 60%, the product is assumed to fail. These two failure processes are both competitive and independent of each other. In a competitive failure model, the lifetime of the system is the least of many random lifetimes.
To analyze the above data, Huang and Askin [34] assumed that both the sudden and degradation failures are modeled by a Weibull distribution and they discussed reliability analysis of this competing failure model. In their paper, they assumed that the population is homogeneous and describe the degradation process by assuming that the light intensity level at each inspection time follows a Weibull distribution whose shape and scale parameters are time dependent. Both of the shape and scale parameters are estimated by the degradation levels observed at each time. Zhao and Elsayed [37] assumed that the sudden failure time follows a Weibull distribution and the degradation failure process is modeled by a Brownian motion, and they used the maximum likelihood estimate (MLE) method to obtain the estimates of the model. Studies based on Huang and Askin [34] and Cha et al. [38] assume that the large heterogeneity observed during degradation is described in part by considering two distinct subpopulations and using least squares estimation to obtain the main reliability features.
Reliability is often closely related to system security. Hence, reliable inference procedures for competing failure model studies with small sample cases have become an important issue in reliability analysis. The challenge of providing reliable inference procedures based on small samples inspires us to explore interval estimation approaches for competing failure models. In this paper, we propose a Wiener-weibull competing failure model and develop the generalized pivotal quantity (GPQ) method to explore the interval estimation of system's reliability metrics under small sample case, and use the proposed model and method to analyze the data in the motivated example.
The rest of the paper is arranged as follows. In Section 2, we outline the Wiener-weibull competing failure model. In Section 3, the MLEs and inverse estimates (IEs) of model parameters are derived and the exact confidence intervals (ECIs), generalized confidence intervals (GCIs) for model parameters, and some important reliability metrics such as the pth quantile of lifetime, the reliability function, and the mean time to failure (MTTF) of system are developed. In Section 4, Monte Carlo techniques are used to examine the performance of the proposed GCIs in terms of the coverage percentage (CP) and average interval length (AL). In Section 5, an illustrative example is given to apply the proposed method. Finally, we summarize the article in Section 6.
2.
Wiener-weibull competing failure model
Supposed that the system is equipped with two groups of components: The first group contains a component, whose degradation process of quality characteristic is described as a stochastic process; the second group contains a component, whose lifetime is modeled by sudden failure. Moreover, the two components are operating independently. In this paper, we assume that the degradation process of quality characteristic for the first component is modeled by a Wiener process, and the lifetime of the second component due to sudden failure follows a Weibull distribution.
2.1. Wiener degradation process
It is assumed that the degradation path of the quality characteristics of the first component can be fitted using the Wiener process {X(t),t≥0}, denoted by
where μ and σ>0 are the drift and diffusion parameters, respectively, μ reflects the degradation rate, and B(⋅) denotes a standard Brownian motion. The Wiener process X(t) has the following properties:
● X(0)=0 is true with probability one.
● X(t|t≥0) has independent increments, that is, the increments X(t1)−X(t0),…,X(tn)−X(tn−1) are independent random variables for ∀0<t0<t1<⋯<tn−1<tn.
● Each increment, ΔX(t)=X(t+Δt)−X(t), follows a normal distribution N(0,σ2Δt).
The lifetime T1 of the first component is defined as the first hitting time of X(t) to a degradation threshold L. As is known to all, T1 follows the IG distribution IG(L/μ,L2/σ2). Therefore, the cumulative distribution function (CDF) of T1 is presented as
where Φ(⋅) is the CDF of N(0,1) distribution.
2.2. Weibull sudden failure mode
Suppose that the lifetime T2 of the second component due to sudden failure follows a Weibull distribution, denoted by Weibull(β,η). The probability density function (PDF) of T2 is
and the CDF of T2 is
where η>0 and β>0 are the scale and shape parameters, respectively.
Therefore, the lifetime of the system can be defined as T=min(T1,T2). The CDF of T and the reliability function of the system at time t are presented as
respectively.
The MTTF of the system can be obtained by
2.3. Data of competing failure model
Suppose that n systems are tested. Let ri denote the number of measurements for the first component of the ith system. The measurement times for the first component of the ith system ti,1,ti,2,…,ti,ri,1≤i≤n, are usually predetermined. Therefore, the degradation data is X={X(ti,j);i=1,2,…,n,j=1,2,…,ri}. Let Δti,jˆ=ti,j−ti,j−1, ΔXi,jˆ=X(ti,j)−X(ti,j−1) denote the degradation increment between ti,j−1 and ti,j, for i=1,2,…,n,j=1,2,…,ri. For convenience of expression, let T=∑ni=1∑rij=1Δti,j and M=∑ni=1ri denote the total test duration and the total number of measurements for the whole test, respectively. The sudden failure time of the second component for the ith system is Ti,2, i=1,2,…,n. Hence, the sudden failure times of the second component for n systems are T=(T1,2,T2,2,…,Tn,2).
The degradation data refers to 10 electronic devices whose degradation level (brightness) was measured at the same inspection times, with equal inspection time interval of 500 hours and the test duration up to 4,000 hours. Suppose that the degradation process {X(t);t≥0} is a Wiener process with drift parameter μ and diffusion parameter σ2. As the degradation level reaches to (or exceeds) the threshold level L, the device is considered as a fail. Where the degradation level is X(ti,j)=100−Yi(ti,j), Yi(ti,j) denotes the light intensity (in percentage relative to the original measurement) of the test unit i at time ti,j. The sudden failure data and the transformed degradation data for test units are given in Tables 2 and 3.
3.
Estimation for competing failure model
In this section, we first give the MLE of parameters μ and σ2 for the Wiener degradation process. On basis of the MLEs of μ and σ2, the ECIs of μ and σ2 are obtained. Unfortunately, to get the confidence interval of the scale parameter η as intractable, we develop the GCIs of parameter η for the sudden failure model. It is well known that the pth quantile of system lifetime, the reliability function, and the MTTF of a system are three important characteristics in reliability analysis. However, it is intractable to obtain the ECIs of these three reliability characteristics, so we consider getting the GCIs of them.
3.1. Estimation for Wiener degradation model
Notice that the degradation increments of quality characteristic ΔXi,j are mutually independent, and ΔXi,j∼N(μΔti,j,σ2Δti,j) for i=1,2,…,n,j=1,2,…,ri. Hence, on basis of the degradation data X, the likelihood function is expressed as
Therefore, the MLEs of parameters μ and σ2 are obtained by
respectively.
Next, we will develop the ECIs for parameters μ and σ2. To derive the ECIs of μ and σ2, the following Theorem 1 is needed.
Theorem 3.1. Suppose that the degradation increments D={ΔXi,j;i=1,2,…,n,j=1,2,…,ri} are from the Wiener degradation process {X(t);t≥0} above. Let ˆμ=∑ni=1∑rij=1ΔXi,j/T, S2=1M−1∑ni=1∑rij=1(ΔXi,j−ˆμΔti,j)2Δti,j, then
1) ˆμ is an unbiased estimator of μ, and ˆμ∼N(μ,σ2/T);
2) S2 is an unbiased estimator of σ2, and (M−1)S2/σ2∼χ2(M−1);
3) S2 is independent of ˆμ.
Proof Notice that ∑ni=1∑rij=1ΔXi,j∼N(μT,σ2T), so ˆμ∼N(μ,σ2/T) is obvious. By telescoping ΔXi,j−μΔti,j as (ΔXi,j−ˆμΔti,j)+(ˆμΔti,j−μΔti,j), we have the following factorization:
According to Cochran [39], ˆμ and S2 are independent and (M−1)S2/σ2∼χ2(M−1).
Using the results of Theorem 1, the 100(1−γ)% ECIs of μ and σ2 are obtained by
respectively, where tγ(n) and χ2γ(n) are the lower γ percentiles of t and χ2 distributions with free degrees n, respectively.
3.2. Estimation for Weibull sudden failure model
In this subsection, for the Weibull sudden failure model, we will give the ECI of shape parameter β. Moreover, for point estimation, the IEs of parameters η and β are obtained. To construct the ECI for parameter β, the following Lemmas 1 and 2 tend to be useful.
Lemma 3.1. Suppose that Y1,Y2,…,Yn are independent identically distributed (i.i.d) random variables from Weibull distribution (2). Let Zi=(Yiη)β,i=1,2,…,n, then the Z1,Z2,…,Zn are independent standard exponential variables.
Lemma 1 is obvious, so here we neglect the detailed proof. □
Lemma 3.2. Given that Z1,Z2,…,Zn are standard exponential random variables and Z(1),Z(2),…,Z(n) are their order statistics, let ξ1=nZ(1),ξi=(n−i+1)(Z(i)−Z(i−1)),i=2,3,…,n; Si=∑ij=1ξj,U(i)=Si/Sn,i=1,2,…,n−1, and Sn=∑ni=1ξi, then
1) ξ1,ξ2,…,ξn are independent standard exponential random variables;
2) U(1)<U(2)<⋯<U(n−1) are the corresponding order statistics of uniform distribution U(0,1) with sample size n−1;
3) 2Sn follows the distribution χ2(2n).
Proof 1) As is known to all, the joint probability density function (JPDF) of (Z(1),Z(2),…,Z(n)) is
Notice that ∑ni=1ξi=∑ni=1Z(i) and the Jacobian determinant is J=|∂(Z(1),Z(2),…,Z(n))∂(ξ1,ξ2,…,ξn)|=1n!, so the JPDF of (ξ1,ξ2,…,ξn) is obtained by
That is, ξ1,ξ2,…,ξn are independent standard exponential random variables.
2) From U(i)=Si/Sn,i=1,2,…,n−1, we know that ξ1=U(1)Sn,ξn=Sn−U(n−1)Sn and ξi=U(i)Sn−U(i−1)Sn,i=2,…,n−1. As the Jacobian determinant J=|∂(ξ1,ξ2,…,ξn)∂(U(1),…,U(n−1),Sn)|=Sn−1n, the JPDF of (U(1),…,U(n−1),Sn) is given by
By marginal integral, the JPDF of (U(1),…,U(n−1)) is obtained by
Hence, U(1)<U(2)<⋯<U(n−1) are the corresponding order statistics of uniform distribution U(0,1) with sample size n−1.
3) Notice that Sn∼Ga(n,1), then we have 2Sn∼Ga(n,1/2)=χ2(2n). □
Next, we will construct pivotal quantities (PQs) for parameters β and η. Since the sudden failure data T is a sequence from the Weibull distribution (2), the corresponding order failure data is denoted by {T(1),2,T(2),2,…,T(n),2}. Based on Lemma 1, we know that the transformation {(Ti,2/η)β,i=1,2,…,n} is a sequence of standard exponential random variables. Thus, from Lemma 2, we have that
are order statistics of the uniform distribution U(0,1).
For shape parameter β, consider the following PQ
From Eq (3.1), we find that for Weibull distribution (2), W1 is a function with respect to the shape parameter β and free of the scale parameter η.
It is obvious that W1 is nonnegative. Notice that ∑n−1i=1logU(i)=∑n−1i=1logUi and Ui,i=1,2,…,n−1 are i.i.d random variables from the uniform distribution U(0,1). Moreover, we can prove the fact that W1∼χ2(2n−2).
Next, we will prove that W1 is strictly monotonic with respect to parameter β. Let Q(j,i)=(T(j),2/T(i),2)β. Note that
It can be observed from Eq (3.2) that W1 is strictly increasing with respect to parameter β, because Q(j,i) is strictly increasing (decreasing) for j>i (j<i). Hence, given a realization W1 from χ2(2n−2), there exists a unique solution g(W1,T) of β for Eq (3.1), then the PQ for parameter β is given as P1=g(W1,T). Therefore, an ECI of β for the Weibull distribution can be derived by the following Theorem 2.
Theorem 3.2. If T1,2,T2,2,…,Tn,2 are i.i.d random variables from Weibull distribution (2),
T(1),2,T(2),2,…,T(n),2 are the corresponding order statistics, then for any 0<γ<1,
is a 1−γ level confidence interval of the shape parameter β. Here χ2γ(n) denotes the lower γ percentile of the χ2 distribution with freedom degrees n, and for t>0, W−11(t) is the solution of β for the equation W1(β)=t.
Notice that W1∼χ2(2n−2) and E(W1)=2(n−1). So, W1 converges to 2(n−1) with probability one. Let W1=2(n−1). Based on the following Eq (3.3), we can get the point estimator ˆβ of the shape parameter β
similar to the discussion above. Eq (3.3) also has a unique solution for parameter β.
Denote An=∑ni=1(Ti,2/η)β, so An∼Ga(n,1) and E(An)=n. Similarly, let An=n, and the corresponding point estimator ˆη of η is obtained by
The estimators obtained from Eqs (3.3) and (3.4) are named as IEs of parameters β and η, which was proposed in Wang [40].
3.3. GCIs for η,Tp,R(t0), and MTTF
In practical applications, some reliability metrics of a system, such as the pth quantile of lifetime, the reliability function R(t0), and the MTTF of system, may be of more importance than the model parameters. However, since these reliability metrics involve multiple parameters, it is intractable to obtain their exact confidence intervals. Therefore, we develop the GCIs for these reliability metrics.
Now we will construct the GPQ for the scale parameter η. Based on Lemmas 1 and 2, we know the quantity
then the scale parameter η can be expressed as η=(2∑ni=1Tβ(i),2/W2)1/β. Recall that the PQ of β is P1=g(W1,T). Using the substitution method given by Weerahandi [41], we replace β by P1 in the expression of η and obtain the GPQ of parameter η
It can be observed from Eq (2.3) that for Weibull sudden failure model, the reliability is R2(t0)=1−F2(t0|η,β). Using the substitution method, the GPQ of reliability R2(t0) is obtained by
To derive the GPQ for reliability R1(t0) of the Wiener degradation model, we first construct the PQs of μ and σ.
Let
Obviously, U∼N(0,1) and V∼χ2(M−1) and they are mutually independent. Thus, μ and σ can be formulated as
respectively, so the GPQs of μ and σ are obtained by
It should be pointed out that ˆμ and S2 are treated as known quantities in generalized inference [41]. Using the substitution method given in [41], the GPQ of R1(t0) is given by
Based on Eqs (2.4)–(2.6), the GPQs for pth quantile of lifetime T, the reliability function, and the MTTF of a system can be obtained by
respectively.
Let Pi,γ denote the γ percentile of Pi, then [Pi,γ/2,Pi,1−γ/2],i=2,5,6,7 are the 1−γ level GCIs of η,Tp,R(t0), and MTTF, respectively. The percentiles of Pi,i=2,5,6,7 can be acquired through the following Monte Carlo Algorithm.
Algorithm : The percentiles of η,Tp,R(t0), and MTTF.
Step 1 Given degradation data X and sudden failure data T, compute ˆμ,S2, and T.
Step 2 Generate W1∼χ2(2n−2),W2∼χ2(2n),U∼N(0,1), and V∼χ2(M−1), respectively.
Based on Eqs (3.1), (3.5), and (3.7), compute P1,P2,P3, and P4.
Step 3 Based on P1,P2,P3, and P4, using Eqs (3.8)–(3.10) to compute P5,P6, and P7
Step 4 Repeat steps (2) and (3) K times, then K values of Pi,i=2,5,6,7 are obtained, respectively.
Step 5 Sorting all Pi values in ascending order: Pi,(1)<Pi,(2)<…<Pi,(K),i=2,5,6,7, then the γ percentile of Pi is estimated by Pi,(γK).
4.
Simulation study
The Monte Carlo simulation technique is used to evaluate the performance of the proposed GCIs in the aspect of the CP and AL. Table 4 lists the different combinations of the model parameters μ,σ2,η,β and the threshold L for simulation study. Moreover, we take n=10,15,20, riˆ=r=8,10,12, Δti,jˆ=Δt=10, and K=10,000 in the simulation study. Based on 5000 replications, all the simulation results are provided in Tables 5–9.
To examine the performance of the point estimates of model parameters (μ,σ2,β,η), simulation studies were carried out in terms of relative-bias (R-Bias) and relative-mean square error (R-MSE) under the parameter setting II, III, and IV for (n,r)=(10,8),(15,10),(20,12). Motivated by Luo et. al. [42], the R-Bias and R-MSE are defined as:
where θ and ˆθi are the true value and the estimate of a parameter, respectively.
Based on 5000 replications, the simulation results about the R-Bias and R-MSE of the proposed estimates for model parameters (μ,σ2,β,η) are provided in Table 5. It can be seen from Table 5 that both the R-Bias and R-MSE are small compared with the true values, and for given parameter settings as n and r increase, the R-MSEs decrease as expected. The simulation results show that these estimates perform well in aspect of both R-Bias and R-MSE. Hence, we recommend the proposed point estimates for model parameters (μ,σ2,β,η).
It is well known that the parametric bootstrap method is a classic approach to get confidence intervals for model parameters. In order to fully assess the performances of the GCIs, we also considered the bootstrap CIs for the proposed competing failure model. A comparative analysis is conducted between the CIs obtained by the GPQ method and the bootstrap-p method. Based on 5,000 bootstrap samples, the bootstrap-p CIs are obtained and provided in Tables 6–9.
It is observed from Tables 6–9 that the CPs of the proposed GCIs are quite close to the nominal levels, even in the small sample case. In many cases, the differences between the real CP and the nominal level are small, ranging in 1%, but the CPs obtained by the bootstrap-p method are away from the nominal levels for some parameters and quantities. In particular, we find that the bootstrap-p CIs of the scale parameter η are far below the nominal levels for all cases. In addition, from Tables 6 and 8 we also find that the bootstrap-p CIs perform bad for some quantities. For example, the CPs of lifetime quantile T0.1 and reliability function R(5) deviate from the nominal levels.
When the sample size n turns large, the CPs of the bootstrap-p CIs also near the nominal levels. Tables 6–9 report that, for fixed parameter settings, when n and r increase, the ALs become shorter for both GPQ and bootstrap-p CIs as expected. The comparison shows that the GCIs perform better than the corresponding bootstrap-p CIs in terms of CP. Hence, we recommend the GCIs for model parameter η and some quantities, such as Tp,R(t0),MTTF, particularly in the case of small sample.
5.
An illustrative example
In this section, we use the proposed Wiener-Weibull competing failure model and the GPQ method to analyze the real example provided by Huang and Askin [34]. The product is treated as a fail when the luminance ratio decreases by 60%. For convenience, we assumed that the original brightness is 100, so the threshold is L=60. The transformed degradation data of luminance ratio is presented in Table 3 and the weld interface fracture data (sudden failure data) is given in Table 2. Figure 1 shows the luminance ratio degradation paths of 10 test units. In this study, we use the Wiener process to model the degradation data in Table 3. Similar to Huang and Askin [34], the sudden failure data is also fitted by a Weibull distribution.
For point estimation, the MLEs of μ and σ2 are given by ˆμ=0.0310,ˆσ2=0.0076, respectively. The IEs of η and β are given by ˆη=1.0350×103,ˆβ=4.7684, respectively. Figure 2 shows the degradation paths, the sample average degradation path, and the fitted mean degradation path by model. It is clear that the estimates of the mean degradation paths provide good fitting for the sample average degradation paths. This means that it is reasonable to use the Wiener process to fit the luminance ratio degradation data. Given p=0.1,L=60, and t0=800(days), the pth percentile of lifetime, the reliability function at time t0, and the MTTF of the system are obtained by T0.1=642.8386(days), R(800)=0.7416, MTTF=943.5986(days), respectively. The point estimate of ^MTTF=22646(hours) is near to the estimate of MTTF (22,765 hours) provided by Huang and Askin [34].
For interval estimation, we use the GPQ method proposed in Section 3 to analyze the real dataset. As is known to all, some reliability metrics such as the pth percentile of lifetime T, the reliability function, and the MTTF of a system are more important than the model parameters in reliability analysis. However, as these metrics are very complicated, getting their interval estimations is usually difficult, so we developed the GCIs for them. Based on Eqs (3.8)–(3.10) and using the GPQ method, the GCIs of Tp,R(t0), and MTTF can be obtained. Take K=10000; the results are given in Table 10. According to the methods in Huang and Askin [34] and Cha et al. [38], they can only provide the point estimation for MTTF and not give its interval estimation.
6.
Conclusions
In this paper, a competing failure model involving both degradation failure and sudden failure was studied by modeling degradation failure as a Wiener process and sudden failure as a Weibull distribution. For model parameters, the MLEs of μ,σ2, and the IEs of η,β were derived and the ECIs of parameters μ, σ2, and β are obtained.
In this study, the GPQ method was proposed to investigate the scale parameter and some reliability metrics of the competing failure model. By constructing the GPQs, the GCI of parameter η was developed. Using the substitution method, the GCIs for the reliability function, the pth percentile of the lifetime, and the MTTF of a system were also developed. Simulation studies were carried out to assess the performances of the proposed intervals. Simulation results reported that the proposed GCIs work well in aspect of the CP and are better than the corresponding bootstrap CIs. Finally, we applied the proposed model and the GPQ method to a real data example and found the ECIs and GCIs of model parameters and some reliability metrics.
Use of AI tools declaration
The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.
Acknowledgements
The authors thank the Editor and the reviewers for their detailed comments and valuable suggestions, which helped improve the manuscript significantly. The work was supported in part by the National Social Science Foundation of China (21CTJ005), in part by the Talent Cultivation and Research Start-up Foundation of Anhui Polytechnic University (S022022014), in part by the Pre-research Project of the National Science Foundation of Anhui Polytechnic University (Xjky08201903), in part by the program for outstanding young talents in colleges and universities of Anhui Province (gxyqZD2022046), in part by the Excellent and Top-notch Personnel Cultivation Project of Universities (gxyq2022083), and in part by the Natural Science Research Project of Anhui Educational Committee (2023AH050927).
Conflict of interest
The authors declare that there is no conflict of interest.