# Reliability Modeling for Systems with Multiple Degradation Processes Using Inverse Gaussian Process and Copulas

Authors: Zhenyu Liu; Xiaobing Ma; Jun Yang; Yu Zhao
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

We develop a reliability model for systems with s-dependent degradation processes using copulas. The proposed model accommodates assumptions of s-dependence among degradation processes and allows for different marginal distributions. This flexibility makes the model more attractive compared with the multivariate distribution model, which lay on the limitation of the homogeneous marginal distribution and can only describe linear correlation. Marginal degradation process is modeled by the inverse Gaussian (IG) process with time scale transformation. Furthermore, we incorporate random drift to account for the possible heterogeneity in population. This paper also develops the statistical inference method using EM algorithm with two-stage procedure. The comparison results of the reliability estimation under both s-dependent and s-independent assumptions are illustrated in the illustrative example to demonstrate the applicability of the proposed method.


---

## Body

## 1. Introduction

For systems with high reliability, it is challenging to do reliability assessment in a timely manner because usually minor failures occur during such short time. Thus it is difficult to assess reliability based only on limited lifetime data. Compared with traditional lifetime data, degradation data contains more information, which records the accumulation of damage over time. Degradation analysis aims to characterize the underlying failure process and uses more information than lifetime data analysis. In reality, a system may have multiple degradation processes or consists of multiple degrading components. The reliability of a system is usually calculated under the assumption of s-independence. However, assuming s-independence among degradation processes may underestimate system reliability. It is more realistic to assume some sort of dependence among degradation processes.In practice, the continuous degradation process is often modeled by a stochastic process. Wiener process and the Gamma process are two classes of stochastic processes that have been widely used in degradation modeling [1]. A distinct feature of the Wiener process is that its degradation measures are not necessarily monotone [2], which is not applicable in many cases. As an alternative, the Gamma process is often used when monotonicity is required [3]. Although the Wiener process and Gamma process have wide applications in degradation modeling, only these two classes of models cannot fit all degradation data well. Wang and Xu [4] showed an application in which both Wiener and Gamma processes do not fit the data well. They try to fit the GaAs Laser data by using Wiener and Gamma process models, but probability plots show that both of them fit the data poorly. They suggest inverse Gaussian (IG) process as another good choice for degradation modeling which provides a monotone degradation path. The IG process was proposed by Wasan [5] in 1968. Notwithstanding wide applications of the IG distribution, the IG process was scarcely used in degradation modeling. Ye and Chen [6] attributed the reason of its unclear physical meanings to reliability engineers and prove that IG process is the limit of a compound Poisson process whose arrival rate goes to infinity while the jump size goes to zero in a certain way. Their research provides a physical meaning for the IG process as a degradation model.For multiple degradation processes, Wang and Coit [7] introduced a general modeling approach to estimating the reliability of a system using multivariate normal distribution under the condition that degradation mechanism is unknown. Pan and Balakrishnan [8] introduced a reliability model for systems with bivariate degradation processes by utilizing bivariate Birnbaum-Saunders distribution. In their model, the degradation process is described by Gamma process and distribution of the first passage time is approximated by Birnbaum-Saunders distribution. In recent years, much attention has been placed on modeling the s-dependence behavior between degradation variables by copulas, which allows us to couple degradation variables from different distribution families. Sari et al. [9] introduced a two-stage model for bivariate processes using the copula function. Experiment data from LED lighting systems are analyzed to illustrate the application of the proposed model. Pan et al. [10] proposed a bivariate Wiener degradation process to describe the dependence between degradation characteristics, where marginal distribution and multivariate dependence can be modeled separately, and marginal distribution functions of degradation measures may follow different distributions families.Based on the bivariate Wiener degradation model, Wang et al. [11] put forward an adaptive method for residual life estimation. In their method, the dependence of degradation variables is characterized by the Frank copula function. Hao and Su [12] presented a bivariate nonlinear Wiener degradation model and use Markov Chain Monte Carlo (MCMC) method to estimate the model parameters. Wang and Pham [13] established a s-dependent competing risk model to handle the dependence among multiple degradation processes and random shocks using copulas. Two types of dependence are considered in their model: dependence among degradation processes and dependence between degradation and shocks.Previous research has been focused on utilizing multivariate distributions. This approach has some drawbacks. One major drawback is that it limits the distribution of each degradation variable. Degradation variables should have the same distribution. Meanwhile, such approaches are applicable only when degradation variables are linearly correlated. Some researchers use copulas to model the dependence among degradation processes, which allows for different marginal distributions. A more interesting problem that has not been addressed in the literature is the effect of copulas on s-dependence modeling for systems with IG processes. This paper makes a further study on s-dependence modeling using IG process and copulas. We first use the IG process with random drift and time scale transformation to model the monotonic degradation process; then we employ the copula method to fit the joint distribution of multiple degradation processes. We also develop the statistical inference method using EM algorithm with two-stage procedure to estimate model parameters.The remainder of the paper is organized as follows. Section2 introduces an IG process with random drift and time scale transformation. Section 3 proposes a system reliability model under s-dependent assumption using copulas. The statistical inference for the model parameters is discussed in Section 4. In Section 5, a simulation study is conducted to demonstrate the quality of the estimators. Section 6 gives an example of crack growth problem to demonstrate the proposed method. Section 7 concludes the paper.


## 2. IG Process with Random Drift and Time Scale Transformation

In practice, many systems or components have their performance characteristics whose degradation over time can be measured. The degradation path over time can be modeled by a stochastic process{X(t);t≥0}. In our study, we use the inverse Gaussian (IG) process to describe monotonic degradation process.Definition 1.
The inverse Gaussian process is defined as the stochastic process satisfying the following.(i)
X(t+s)-X(s) follows an inverse Gaussian distribution with mean t and scale σt2.(ii)
X(t) has independent increments on nonoverlapping intervals; that is, X(t4)-X(t3) and X(t2)-X(t1) are independent, ∀t1<t2≤t3<t4.Consider that the IG process{X(t);t≥0} with X(t)-X(s) follows IG(μ(τ(t)-τ(s)),σ(τ(t)-τ(s))2), where τ(·) is the transformed time scale; the transformation depends on particular mechanism. General discussions about time transformation can be found in Whitmore and Schenkelberg [15]. μ is the drift parameter; σ is the diffusion parameter. Let τ(0)=0 and X(0)=0; thus X(t) follows IG(μτ(t),στ2(t)) with mean μτ(t) and variance μ3τ(t)/σ, where IG(μτ(t),στ2(t)) denotes the IG distribution with probability density function (PDF):
(1)fIG(x)=στ2(t)2πx3exp⁡[-σ(x-μτ(t))22μ2x],
and cumulative distribution function (CDF)
(2)FIG(x)=Φ[σx(τ(t)-xμ)]-exp⁡(2στ(t)μ)Φ[-σx(τ(t)+xμ)].In many engineering applications, the failure of a system is defined when the degradation processX(t) crosses a critical threshold level D for the first time. The corresponding time TD is
(3)TD=inf⁡{t:X(t)>D}.Because IG process has monotonic property, the reliability function of the first passage timeTD can be obtained from the relation P(TD>t)=P(X(t)<D)=FIG(D) as
(4)R(t)=Φ[-σD(τ(t)-Dμ)]+exp⁡(2στ(t)μ)Φ[-σD(τ(t)+Dμ)].Suppose that a system hasK degradation processes and each of them is governed by an IG process. We assume that all samples have common inspection times. The inference can be easily extended to the case where inspection times are different. In a degradation test, N units are tested and the degradation measures of all units are observed at ordered inspection times t1<t2<⋯<tM with observations xp(i)(tj), p=0,1,2,…,K, i=1,2,…,N, j=0,1,2,…,M, where K is the total number of degradation processes, N is the total number of test units, and M is the total number of inspection times. In general, the degradation data can be presented in the form
(5)xKN×M=[x1⋮xK]=[x1(1)(t1)x1(1)(t2)⋯x1(1)(tM)⋮⋮⋱⋮x1(N)(t1)x1(N)(t2)⋯x1(N)(tM)⋮⋮⋱⋮xK(1)(t1)xK(1)(t2)⋯xK(1)(tM)⋮⋮⋱⋮xK(N)(t1)xK(N)(t2)⋯xK(N)(tM)].Because of the independent property of IG process, increments for nonoverlapping intervals are independent. Therefore, we focus on the degradation increment. Definet0=0, xp(i)(t0)=0, xpij=xp(i)(tj)-xp(i)(tj-1), and τpj=τp(tj)-τp(tj-1) for p=1,2,…,K, i=1,2,…,N, j=1,2,…,M. We have
(6)xpij~IG(μpτpj,σpτpj2),
where μp and σp are parameters for pth degradation process.The joint PDF ofxpi={xpi1,xpi2,…,xpiM} is given by
(7)f(xpi)=∏j=1Mσpτpj22πxpij3exp⁡[-σp(xpij-μpτpj)22μp2xpij].Considering heterogeneity within the population, random effects are often used to account for heterogeneous degradation rates across units [16, 17]. Normality assumption for μp has been adopted by many studies [18]. To avoid the negative values of μp and simplify the deviation, we assume that μp-1 conditional on 0≤μp-1<∞ follows a truncated normal distribution TN(ωp,δp-2). This means that μp-1 is fixed but unknown for each sample. Its probability density function (PDF) for 0≤μp-1<∞ is given by
(8)g(μp-1)=δp·ϕ[δ(μp-1-ωp)]1-Φ(-δpωp),
where ϕ(·) is the PDF of standard normal distribution, Φ(·) is the CDF of standard normal distribution, and p=1,2,…,K. Note that, conditional on μp-1, the unconditional PDF of xpi can be obtained by integrating μp-1 out of (7), which yields
(9)f(xpi)=∫0∞(∏j=1Mσpτpj22πxpij2exp⁡[-σp(xpij-μpτpj)22μp2xpij]hhhhhhhh·g(μp-1)∏j=1Mσpτpj22πxpij2exp⁡[-σp(xpij-μpτpj)22μp2xpij])dμp-1.The random drift model covers simple model without random effects whenδp→∞. The time scale transformation can also be eliminated if we let τp(t)=t. When both δp→∞ and τp(t)=t, the random drift model with time scale transformation degenerates to the regular IG process. This means that the random drift model with time scale transformation can be viewed as an extended model which includes several IG process models as special cases. In degradation data analysis, an extended model is a much flexible choice to fit a specific dataset.


## 3. Modeling for Systems with Multiple Degradation Processes

Although, in many studies, multiple degradation processes are assumed to have independent property, it may be more realistic to assume some sort of dependence among degradation processes. According to Sklar’s theorem [19], if H is a joint distribution function with continuous margins F1,F2,…,Fn, then there exists a unique copula C such that, for all x1,x2,…,xn in R,
(10)H(x1,x2,…,xn)=C(F1(x1),F2(x2),…,Fn(xn)).This theorem provides the theoretical foundation for the application of copulas. Based on Sklar’s theorem, any multivariate distribution can be decomposed into a copula and its marginal. Thus, copula functions provide a much more flexible and realistic way to study multivariate distributions. Table1 lists some well-known one-parameter Archimedean copulas, in which u and v are marginal distributions.Table 1
Summary of some one-parameter Archimedean copulas.
Copulas
C
(
u
,
v
)
Parameter
Gumbel copula
exp
⁡
{
-
[
(
-
ln
⁡
u
)
α
+
(
-
ln
⁡
v
)
α
]
1
/
α
}
α
∈
[
1
,
∞
)
Frank copula
-
1
α
ln
⁡
(
1
+
(
e
-
α
u
-
1
)
(
e
-
α
v
-
1
)
e
-
α
-
1
)
α
∈
(
-
∞
,
∞
)
∖
{
0
}
Clayton copula
[
max
⁡
(
u
-
α
+
v
-
α
-
1,0
)
]
-
1
/
α
α
∈
[
-
1
,
∞
)
∖
{
0
}
Joe copula
1
-
[
(
1
-
u
)
α
+
(
1
-
v
)
α
-
(
1
-
u
)
α
(
1
-
v
)
α
]
1
/
α
α
∈
[
1
,
∞
)
Amh copula
u
v
1
-
α
(
1
-
u
)
(
1
-
v
)
α
∈
[
-
1,1
]Properties of the copula have been explored in a number of studies; for example, see Jia et al. [20].Consider a system subject to multiple s-dependent degradation processes. The degradation processes are denoted as{X1(t),X2(t),…,XK(t)}, where K is the total number of degradation processes and Xp(t) are degradation process for p=1,2,…,K. The joint distribution of {X1(t),X2(t),…,XK(t)} is
(11)H(X1(t),X2(t),…,XK(t))=C(F1(X1(t)),F1(X1(t)),…,FK(XK(t))),
where Fp(Xp(t)) is the CDF of Xp(t) for p=1,2,…,K.The system is considered to be failed if either degradation processXp(t) reaches its corresponding failure threshold, which is known as Dp for p=1,2,…,K. Therefore, the reliability of a system with s-dependent degradation processes is
(12)R(t)=P((X1(t),X2(t),…,XK(t))∈S)=P(X1(t)<D1,X2(t)<D2,…,XK(t)<DK)=C(F1(D1),F2(D2),…,FK(DK)),
where S means safe region, within which the system is working and Fp(Dp) is the CDF of Xp(t) for p=1,2,…,K.


## 4. Statistical Inference

### 4.1. Time Scale Determination

In a few situations we understand well the parametric form for time scaleτ(t) based on physics mechanism or prior knowledge. When there is no such failure mechanism to describe the time scale, we use a statistical method to help us determine the parametric form. From (6), we can see that the mean degradation path is E[μ]τ(t), where E[μ] can be obtained by using the property of truncated normal distribution. We may average the samples to obtain the mean degradation curve and specify a parametric form for τ(t) by fitting the mean path. Linear relationship, power law relationship, and exponential law relationship are the most used parametric forms.


### 4.2. Parameter Estimation

Suppose that the multivariate distribution is specified byK margins with CDF Fp and PDF fp, p=1,2,…,K, and a copula C. Let α be the parameter in copula and βp=(σp,γp,ωp,δp) the parameters in the pth degradation process. We will remark that γp is the parameter in time scale τp(t). Now, consider the joint distribution of degradation increments xpij in the interval [tj-1,tj]. The joint PDF can be computed by
(13)f(x1ij,x2ij,…,xKij)=c(F1(x1ij;β1),F2(x2ij;β2),…,FK(xKij;βK);α)×∏p=1Kfp(xpij;βp),
where Fp(·;βp) and fp(·;βp) are the pth marginal CDF and PDF with parameter βp; c(F1(·;β1),F2(·;β2),…,FK(·;βK);α)=(∂KC(F1(·;β1),F2(·;β2),…,FK(·;βK);α)/∂F1(·;β1)∂F2(·;β2),…, ∂FK(·;βK)) is the density of copula C(·;α); and xpij is degradation increment in the interval [tj-1,tj].Using MLE, the full log-likelihood function of the copula with marginal function can be expressed as(14)ln⁡L(α,β)=∑i=1N∑j=1M[∑p=1Kln⁡fp(xpij;βp)ln⁡c(F1(x1ij;β1),F2(x2ij;β2),…,hhhhhhhhhFK(xKij;βK);α)+∑p=1Kln⁡fp(xpij;βp)].The parameters can be estimated by maximizing the above log-likelihood function. However, the log-likelihood is computationally difficult to work with, or even infeasible. Here we use a two-stage procedure of firstly estimating the marginal parameters and then estimating the copula parameter from the joint likelihood with the marginal parameters estimated in the first stage. Computational difficulty is greatly reduced since each stage has a very small number of parameters. Joe [21] has studied the efficiency properties of this two-stage estimation procedure. He points out that this method has good efficiency and can sometimes be equivalent to MLE method.Stage 1 (marginal parameters).
Each marginal distributionFp has its own parameter set βp=(σp,γp,ωp,δp). In the first stage, the log-likelihood of the K marginals is separately maximized to get estimates β=(β1,β2,…,βK)T.Based on (9), the log-likelihood function of each degradation process can be expressed as
(15)ln⁡L(βp)=∑i=1N∑j=1Mln⁡fp(xpij;βp),p=1,2,…,K,
where K is the total number of degradation processes, N is the total number of test units, and M is the total number of inspection times.The log-likelihood function is very flat aroundβ^p. Direct maximization of the likelihood function often yields a solution far away from the MLE. We use Expectation-Maximization (EM) algorithm to solve this problem as it is free of the rounding error. We denote μp=[μp1,μp2,…,μpN] as the unobserved data and xp as the observed data. Given the complete data including μp and xp, the log-likelihood of unknown parameters βp=(σp,γp,ωp,δp), up to a constant, can be expressed as
(16)ln⁡L(σp,γp,ωp,δp)=∑i=1N∑j=1M[12ln⁡σp+ln⁡τpj-σp(xpij-μpiτpj)22μpi2xpij]+∑i=1N[ln⁡δp-δp2(μpi-1-ωp)22-ln⁡(1-Φ(-δpωp))].For theE-Step, we need to estimate E(L(σp,γp,ωp,δp)) given the observed data and current estimated parameters. Ye and Chen [6] have proven that the conditional distribution of μpi-1 follows TN((σp∑j=1Mtij+ωpδp2)/(σp∑j=1Mxij+δp2),σp∑j=1Mxij+δp2). By using the property of truncated normal distribution, we can easily obtain api=E[μpi-1∣xpi] and bpi=E[μpi-2∣xpi] for each unit. After substituting them into (16), we have
(17)E(ln⁡L(σp,γp,ωp,δp))=l1(σp,γp)+l2(ωp,δp),
where
(18)l1(σp,γp)=∑i=1N∑j=1M[12ln⁡σp+ln⁡τpj-σp(xpijbpi2-τpjapi+τpj22xpij)],(19)l2(ωp,δp)=∑i=1N[ln⁡δp-δp2(12bpi-apiωp+12ωp2)hhhhhhhhhhhhh-ln⁡(1-Φ(-δpωp))δp2(12bpi-apiωp+12ωp2)].It is noted thatl1(σp,γp) only depends on σp, γp, and l2(ωp,δp) only depends on ωp, δp. So, in the M-Step, we can compute the updated parameter estimates separately. The estimation equations for σp, γp can be easily obtained by maximizing (18); that is,
(20)σp-1=1MN∑i=1N∑j=1M(xpijbpi-2τpjapi+τpji2xpij)∑i=1N∑j=1M[τpj′τpj+σpτpj′(api-τpj′xpij)]=0.We can also obtainωp, δp by solving
(21)Nϕ(-ωpδp)[1-Φ(-ωpδp)]=(∑i=1Napi-Nωp)δp,δp-2=1N∑i=1N(bpi-apiωp).After obtaining updated parameter estimates in theM-Step, we can iterate to obtain new estimates until convergence. It is worth noting that both (20) and (21) contain only two parameters. It is much easier to find the optimal solution, as opposed to the direct maximization of the original log-likelihood function.Stage 2 (copula parameters).
With givenβ^ in Stage 1, the estimation of the parameter α from the copula function can be performed as
(22)α^=argmax⁡α∈Θ∑i=1N∑j=1Mln⁡c(F1(x1ij;β^1),hhhhhhhhhhhhhhhhF2(x2ij;β^2),…,FK(xKij;β^K);α),
where Θ is the parameter space. c(F1(·;β1),F2(·;β2),…,FK(·;βK);α) is the density function derived from copula C(·;α).The copula parameterα can be estimated by maximizing the full log-likelihood function in (22).


### 4.3. Goodness-of-Fit and Model Selection

When there are several potential copula functions available, the most straightforward way to quantitatively select the model is to use the Akaike Information Criterion (AIC). The AIC is defined as(23)AIC=2q-2ln⁡(L),
where q is the number of unknown parameters in the statistical model and ln⁡(L) is the maximized value of the log-likelihood function for the estimated model.Based on this criterion, the model with the smallest AIC value is selected as the best fitting model. However, the model with the smallest AIC does not necessarily be a correct model. We should also check the correlation coefficient between variables. The most widely known correlation coefficients are Pearson’sr, Kendall’s tau, and Spearman’s rho. Because Pearson’s r only measures the linear correlation between random variables, it is not a proper measure of association in this study. We only focus on Kendall’s tau and Spearman’s rho, which does not rely on linear correlated restriction. When applying these correlation coefficients to the bivariate degradation processes, Kendall’s tau can be stated as
(24)τ=4∬[0,1]2C(u,v)dC(u,v)-1,
and Spearman’s rho can be stated as
(25)ρ=12∬[0,1]2C(u,v)dudv-3,
where u and v are the CDFs of marginal random variables.The model that fits the data well is expected to yield the maximum Kendall’s tau and Spearman’s rho. More details could be found in Nelsen [19].


## 4.1. Time Scale Determination

In a few situations we understand well the parametric form for time scaleτ(t) based on physics mechanism or prior knowledge. When there is no such failure mechanism to describe the time scale, we use a statistical method to help us determine the parametric form. From (6), we can see that the mean degradation path is E[μ]τ(t), where E[μ] can be obtained by using the property of truncated normal distribution. We may average the samples to obtain the mean degradation curve and specify a parametric form for τ(t) by fitting the mean path. Linear relationship, power law relationship, and exponential law relationship are the most used parametric forms.


## 4.2. Parameter Estimation

Suppose that the multivariate distribution is specified byK margins with CDF Fp and PDF fp, p=1,2,…,K, and a copula C. Let α be the parameter in copula and βp=(σp,γp,ωp,δp) the parameters in the pth degradation process. We will remark that γp is the parameter in time scale τp(t). Now, consider the joint distribution of degradation increments xpij in the interval [tj-1,tj]. The joint PDF can be computed by
(13)f(x1ij,x2ij,…,xKij)=c(F1(x1ij;β1),F2(x2ij;β2),…,FK(xKij;βK);α)×∏p=1Kfp(xpij;βp),
where Fp(·;βp) and fp(·;βp) are the pth marginal CDF and PDF with parameter βp; c(F1(·;β1),F2(·;β2),…,FK(·;βK);α)=(∂KC(F1(·;β1),F2(·;β2),…,FK(·;βK);α)/∂F1(·;β1)∂F2(·;β2),…, ∂FK(·;βK)) is the density of copula C(·;α); and xpij is degradation increment in the interval [tj-1,tj].Using MLE, the full log-likelihood function of the copula with marginal function can be expressed as(14)ln⁡L(α,β)=∑i=1N∑j=1M[∑p=1Kln⁡fp(xpij;βp)ln⁡c(F1(x1ij;β1),F2(x2ij;β2),…,hhhhhhhhhFK(xKij;βK);α)+∑p=1Kln⁡fp(xpij;βp)].The parameters can be estimated by maximizing the above log-likelihood function. However, the log-likelihood is computationally difficult to work with, or even infeasible. Here we use a two-stage procedure of firstly estimating the marginal parameters and then estimating the copula parameter from the joint likelihood with the marginal parameters estimated in the first stage. Computational difficulty is greatly reduced since each stage has a very small number of parameters. Joe [21] has studied the efficiency properties of this two-stage estimation procedure. He points out that this method has good efficiency and can sometimes be equivalent to MLE method.Stage 1 (marginal parameters).
Each marginal distributionFp has its own parameter set βp=(σp,γp,ωp,δp). In the first stage, the log-likelihood of the K marginals is separately maximized to get estimates β=(β1,β2,…,βK)T.Based on (9), the log-likelihood function of each degradation process can be expressed as
(15)ln⁡L(βp)=∑i=1N∑j=1Mln⁡fp(xpij;βp),p=1,2,…,K,
where K is the total number of degradation processes, N is the total number of test units, and M is the total number of inspection times.The log-likelihood function is very flat aroundβ^p. Direct maximization of the likelihood function often yields a solution far away from the MLE. We use Expectation-Maximization (EM) algorithm to solve this problem as it is free of the rounding error. We denote μp=[μp1,μp2,…,μpN] as the unobserved data and xp as the observed data. Given the complete data including μp and xp, the log-likelihood of unknown parameters βp=(σp,γp,ωp,δp), up to a constant, can be expressed as
(16)ln⁡L(σp,γp,ωp,δp)=∑i=1N∑j=1M[12ln⁡σp+ln⁡τpj-σp(xpij-μpiτpj)22μpi2xpij]+∑i=1N[ln⁡δp-δp2(μpi-1-ωp)22-ln⁡(1-Φ(-δpωp))].For theE-Step, we need to estimate E(L(σp,γp,ωp,δp)) given the observed data and current estimated parameters. Ye and Chen [6] have proven that the conditional distribution of μpi-1 follows TN((σp∑j=1Mtij+ωpδp2)/(σp∑j=1Mxij+δp2),σp∑j=1Mxij+δp2). By using the property of truncated normal distribution, we can easily obtain api=E[μpi-1∣xpi] and bpi=E[μpi-2∣xpi] for each unit. After substituting them into (16), we have
(17)E(ln⁡L(σp,γp,ωp,δp))=l1(σp,γp)+l2(ωp,δp),
where
(18)l1(σp,γp)=∑i=1N∑j=1M[12ln⁡σp+ln⁡τpj-σp(xpijbpi2-τpjapi+τpj22xpij)],(19)l2(ωp,δp)=∑i=1N[ln⁡δp-δp2(12bpi-apiωp+12ωp2)hhhhhhhhhhhhh-ln⁡(1-Φ(-δpωp))δp2(12bpi-apiωp+12ωp2)].It is noted thatl1(σp,γp) only depends on σp, γp, and l2(ωp,δp) only depends on ωp, δp. So, in the M-Step, we can compute the updated parameter estimates separately. The estimation equations for σp, γp can be easily obtained by maximizing (18); that is,
(20)σp-1=1MN∑i=1N∑j=1M(xpijbpi-2τpjapi+τpji2xpij)∑i=1N∑j=1M[τpj′τpj+σpτpj′(api-τpj′xpij)]=0.We can also obtainωp, δp by solving
(21)Nϕ(-ωpδp)[1-Φ(-ωpδp)]=(∑i=1Napi-Nωp)δp,δp-2=1N∑i=1N(bpi-apiωp).After obtaining updated parameter estimates in theM-Step, we can iterate to obtain new estimates until convergence. It is worth noting that both (20) and (21) contain only two parameters. It is much easier to find the optimal solution, as opposed to the direct maximization of the original log-likelihood function.Stage 2 (copula parameters).
With givenβ^ in Stage 1, the estimation of the parameter α from the copula function can be performed as
(22)α^=argmax⁡α∈Θ∑i=1N∑j=1Mln⁡c(F1(x1ij;β^1),hhhhhhhhhhhhhhhhF2(x2ij;β^2),…,FK(xKij;β^K);α),
where Θ is the parameter space. c(F1(·;β1),F2(·;β2),…,FK(·;βK);α) is the density function derived from copula C(·;α).The copula parameterα can be estimated by maximizing the full log-likelihood function in (22).


## 4.3. Goodness-of-Fit and Model Selection

When there are several potential copula functions available, the most straightforward way to quantitatively select the model is to use the Akaike Information Criterion (AIC). The AIC is defined as(23)AIC=2q-2ln⁡(L),
where q is the number of unknown parameters in the statistical model and ln⁡(L) is the maximized value of the log-likelihood function for the estimated model.Based on this criterion, the model with the smallest AIC value is selected as the best fitting model. However, the model with the smallest AIC does not necessarily be a correct model. We should also check the correlation coefficient between variables. The most widely known correlation coefficients are Pearson’sr, Kendall’s tau, and Spearman’s rho. Because Pearson’s r only measures the linear correlation between random variables, it is not a proper measure of association in this study. We only focus on Kendall’s tau and Spearman’s rho, which does not rely on linear correlated restriction. When applying these correlation coefficients to the bivariate degradation processes, Kendall’s tau can be stated as
(24)τ=4∬[0,1]2C(u,v)dC(u,v)-1,
and Spearman’s rho can be stated as
(25)ρ=12∬[0,1]2C(u,v)dudv-3,
where u and v are the CDFs of marginal random variables.The model that fits the data well is expected to yield the maximum Kendall’s tau and Spearman’s rho. More details could be found in Nelsen [19].


## 5. Simulation Study

To demonstrate the quality of the EM estimators proposed in last section, a simulation study has been carried out. Assumeτ(t)=tγ as time scale. Without loss of generality, we fix σ=10, ω=0.5, δ=1, and γ=1.5. We choose the number of test units to be N=10, 15, and 30 and the number of inspections times to be M=5, 10, 15, and 20. Each unit is inspected at ti=i, for i=1,2,…,M. Under each parameter setting (N, M), the bias and the mean square error (MSE) of the MLE are computed based on 1000 Monte Carlo replications. The results are displayed in Table 2.Table 2
Biases and standard errors of the estimates based on the proposed EM approach.
(N, M)
σ
ω
δ
γ
Bias
MSE
Bias
MSE
Bias
MSE
Bias
MSE
(10, 5)
1.04
10.37
0
.
11
3
0.20
0.37
0.51
-
5
e
-
3
1.4
e
-
3
(10, 10)
0.57
4.62
0.0
86
0.21
0.34
0.47
-
3.8
e
-
3
4.9
e
-
4
(10, 15)
0.31
2.82
0.0
92
0.20
0.33
0.42
-
2.6
e
-
3
2.7
e
-
4
(10, 20)
0.29
2.04
0.0
78
0.19
0.33
0.39
-
1.6
e
-
3
1.6
e
-
4
(15, 5)
0.67
5.13
0.0
82
0.17
0.26
0.25
-
4
e
-
3
8.1
e
-
4
(15, 10)
0.44
2.93
0.0
72
0.17
0.24
0.25
-
3.7
e
-
3
3.1
e
-
4
(15, 15)
0.30
1.97
0.07
9
0.17
0.23
0.21
-
2.9
e
-
3
1.7
e
-
4
(15, 20)
0.23
1.44
0.0
53
0.17
0.22
0.21
-
1.8
e
-
3
1.1
e
-
4
(30, 5)
0.24
2.28
0.0
29
0
.
13
0
.
12
0.09
-
2.5
e
-
3
3.5
e
-
4
(30, 10)
0.23
1.33
0.0
23
0.13
0.11
0.0
9
-
2.4
e
-
3
1.3
e
-
4
(30, 15)
0.22
0
.
8
6
0.0
21
0.13
0.12
0
.
08
-
1.9
e
-
3
7.6
e
-
5
(30, 20)
0.10
0
.
6
1
0.004
0.12
0.12
0.08
-
1.4
e
-
3
4.6
e
-
5From Table2, it is readily observed that the biases and MSEs become smaller as the number of test units (N) and inspections times (M) increase. The result reveals that the parameters of our model can be accurately estimated with reasonable sample sizes.


## 6. Illustrative Example

Lu and Meeker [14] presented a fatigue-crack-growth dataset. Degradation data from 21 samples are collected. The time unit is in million cycles; all samples are measured every 0.01 million cycles from 0 to 0.09. We choose the first 20 samples from the dataset and divide them into two groups. We assume that samples in the first group represent crack A, while samples in the second group represent crack B. The length growth of the cracks can be seen as degradation processes. Considering a system subject to fatigue cracks in two different positions, the system is defined to be failed if the length of crack A crosses 1.6 inch or the length of crack B crosses 1.3 inch. The fatigue-crack-growth data are listed in Table 3 and depicted in Figure 1. A similar idea is also adopted in other study cases to analyze their models [8]. Now, we further apply our method to the dataset.Table 3
Fatigue crack growth data from Lu and Meeker [14].
Crack
Path
Million cycles
0.00
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
A
1
0.90
0.95
1.00
1.05
1.12
1.19
1.27
1.35
1.48
1.64
2
0.90
0.94
0.98
1.03
1.08
1.14
1.21
1.28
1.37
1.47
3
0.90
0.94
0.98
1.03
1.08
1.13
1.19
1.26
1.35
1.46
4
0.90
0.94
0.98
1.03
1.07
1.12
1.19
1.25
1.34
1.43
5
0.90
0.94
0.98
1.03
1.07
1.12
1.19
1.24
1.34
1.43
6
0.90
0.94
0.98
1.03
1.07
1.12
1.18
1.23
1.33
1.41
7
0.90
0.94
0.98
1.02
1.07
1.11
1.17
1.23
1.32
1.41
8
0.90
0.93
0.97
1.00
1.06
1.11
1.17
1.23
1.30
1.39
9
0.90
0.92
0.97
1.01
1.05
1.09
1.15
1.21
1.28
1.36
10
0.90
0.92
0.96
1.00
1.04
1.08
1.13
1.19
1.26
1.34
B
1
0.90
0.93
0.96
1.00
1.04
1.08
1.13
1.18
1.24
1.31
2
0.90
0.93
0.97
1.00
1.03
1.07
1.10
1.16
1.22
1.29
3
0.90
0.92
0.97
0.99
1.03
1.06
1.10
1.14
1.20
1.26
4
0.90
0.93
0.96
1.00
1.03
1.07
1.12
1.16
1.20
1.26
5
0.90
0.92
0.96
0.99
1.03
1.06
1.10
1.16
1.21
1.27
6
0.90
0.92
0.95
0.97
1.00
1.03
1.07
1.11
1.16
1.22
7
0.90
0.93
0.96
0.97
1.00
1.05
1.08
1.11
1.16
1.20
8
0.90
0.92
0.94
0.97
1.01
1.04
1.07
1.09
1.14
1.19
9
0.90
0.92
0.94
0.97
0.99
1.02
1.05
1.08
1.12
1.16
10
0.90
0.92
0.94
0.97
0.99
1.02
1.05
1.08
1.12
1.16The development of crack length over time: (a) crack A, (b) crack B.
(a)
(b)We first apply the IG process with random drift and time scale transformation to fit the fatigue-crack-growth data. Figure1 shows that degradation path is not linear over time. Lu and Meeker [14] suggested that Power Law is appropriate for crack size modeling, and thus τ(t)=tγ is adopted as the transformed time scale. The estimates of marginal parameters are listed in Table 4. The 90% confidence intervals are estimated using the parametric percentile bootstrap method with 2000 bootstrap replication, as shown in Table 4.Table 4
Estimation and the 90% confidence intervals.
Parameter
Crack A
Crack B
Estimation
CI. lower
CI. upper
Estimation
CI. lower
CI. upper
σ
69226
53583
104220
17778
14183
27403
ω
0.074
0.066
0.079
0.149
0.128
0.160
δ
123.628
90.466
268.429
52.161
38.557
140.941
γ
1.349
1.324
0.390
1.257
1.232
1.312Based on (6), the estimated crack length paths is given by
(26)E[X(t)]=E[μ]τ(t)=E[μ]tγ^.To demonstrate the goodness of fit, the mean crack length paths estimated from the above model and from the sample average are compared. If these two lines overlap, it suggests that the model is correct. The mean crack length paths based on (26) are plotted in conjunction with the sample average, as shown in Figure 2. In addition, the result by bivariate Birnbaum-Saunders method from Pan and Balakrishnan [8] is also depicted in the same figure for comparison. Intuitively, it can be seen that our model tally quite well with sample average, indicating that the model fits the data well. From Figure 2 we can also see that crack A has a higher degradation rate than crack B.Figure 2
Estimated mean paths based on sample average and on our model.The marginal reliability of these two cracks obtained from the estimated IG distribution is displayed in Figure3. Although the degradation rate of crack A is higher than crack B, the reliability curve of crack A tends to be higher than crack B. The reason is that crack B is defined to be failed if crack size exceeds 1.3 inch, which is lower than the threshold of crack A.Figure 3
Comparison of the reliability of two cracks.Given the estimates of marginal parameters in Table4, the copula method is applied to fit the dependent crack size data. To choose the most suitable fitting copula, we test six selected copulas using the AIC as the criteria for goodness of fit. The results are summarized in Table 5.Table 5
Goodness-of-Fit for the four copulas.
Copula
Parameter Estimation
AIC
Kendall’s tau
Spearman’s rho
Ranking
Gumbel
1.5150
-
1114.89
0.3399
0.4836
2
Frank
3.9127
-
1118.13
0.3816
0.5487
1
Clayton
0.3544
-
1095.96
0.1505
0.2241
6
Gaussian
0.4880
-
1112.91
0.3246
0.4708
3
Joe
1.7663
-
1112.25
0.2985
0.4294
4
Amh
0.8235
-
1108.12
0.2436
0.3590
5From Table5, we see that the model with Frank copula has the smallest AIC of −1118.13, followed by the Gumbel copula with AIC of −1114.89. Clayton copula is the worst fitting among all the copulas, with AIC of −1095.96. That is, the AIC criterion favors the Frank copula. The model with Frank copula also has both the largest Kendall’s tau and Spearman’s rho. Therefore, Frank copula may be the best choice to describe the dependence between the two degradation processes.Given the estimated parameters, we can derive the system reliability from (12), as shown in Figure 4. By means of the bivariate Birnbaum-Saunders method proposed by Pan and Balakrishnan [8], their estimated reliability curve is also depicted in Figure 4, which is very similar to the s-dependent reliability curve from our method. By comparing these two curves we can further validate the effectiveness of our method. As we can see from this figure, s-dependent reliability curve and s-independent reliability curve are overlapped at early time. Because at the early stage both of the two degradation processes are in a good status, and unlikely to impact on each other. However, the reliability estimated by Pan and Balakrishnan [8] is a little bit high at early time. This evidence tends to suggest that the bivariate Birnbaum-Saunders method may not be the best choice for the data. After a period of time, s-independent reliability curve tends to be slightly lower than s-dependent reliability curve because degradation processes are subject to the same stress condition; small degradation amount of one degradation process tends to occur with small degradation amount of the other degradation process. That is, there exists some sort of dependence between the degradation processes.Figure 4
Comparison of the system reliability under both s-dependent and s-independent assumptions.


## 7. Conclusions

This study has investigated reliability modeling method for systems with s-dependent degradation processes. We first introduced an inverse Gaussian (IG) process with random drift and time scale transformation to describe marginal degradation process. This is an extended family of IG process models that contains several existing models as its limiting cases. Based on copulas, we then developed a system reliability model subject to s-dependent degradation processes. The statistical inference method was developed using EM algorithm with two-stage procedure. In the crack length growth problem, we applied the model to fit the crack size data. The estimation accuracy was compared in a numerical example. The results revealed that our model is suitable when the degradation path is nonlinear. The system reliability tends to be higher than s-independent situation after early stage. This is expected because our model considers the dependence among degradation processes. If we ignore the dependent relationship, significant information loss is expected. Therefore, our model is potentially very useful for multivariate degradation data analysis, because a common way to model for multivariate data has limitation of the same marginal distribution and can only describe linear correlation.



---

*Source: 829597-2014-06-18.xml*