# Bivariate Nonlinear Diffusion Degradation Process Modeling via Copula and MCMC

Authors: Huibing Hao; Chun Su
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

A novel reliability assessment method for degradation product with two dependent performance characteristics (PCs) is proposed, which is different from existing work that only utilized one dimensional degradation data. In this model, the dependence of two PCs is described by the Frank copula function, and each PC is governed by a random effected nonlinear diffusion process where random effects capture the unit to unit differences. Considering that the model is so complicated and analytically intractable, Markov Chain Monte Carlo (MCMC) method is used to estimate the unknown parameters. A numerical example about LED lamp is given to demonstrate the usefulness and validity of the proposed model and method. Numerical results show that the random effected nonlinear diffusion model is very useful by checking the goodness of fit of the real data, and ignoring the dependence between PCs may result in different reliability conclusion.


---

## Body

## 1. Introduction

Due to the advances in material science and manufacturing technology, the lifetime of modern product has become longer and longer; usually it is difficult to obtain enough failure data during research and development period. This has posed great challenge for reliability assessment of these products with traditional reliability analysis approach [1]. In this situation, degradation data can be used as an alternate resource for reliability analysis. Compared with failure data, they are easier to be obtained [2]. In the last decades, degradation data have played a more important role in reliability assessment than ever before.Degradation (e.g., wear, erosion, and fatigue) is very common phenomenon for many electrical and electromechanical products, and it can be described by a continuous performance process in terms of time [3]. In [3], three kinds of methods for degradation data analysis are proposed, that is, linear regression method, degradation path method, and stochastic process method. Considering that stochastic process can provide a flexible way to describe the failure mechanism and characteristics of operating environment, it has been widely used to model the degradation path, including Gamma process [4, 5], Wiener process [6, 7], and Markov process [8].Among those stochastic process models, Wiener process has become very popular for degradation modeling in recent years, such as Tseng et al. [9] using Wiener process to describe the lifetime for LED lamps and Peng and Tseng [10] studying the reliability assessment for GaAs lasers based on Wiener process. A fundamental problem related to Wiener process is that it can only describe a linearly drifted diffusion process. However, nonlinearity exists extensively in practice, and the linear model cannot trace the dynamics of the nonlinearity degradation process. Therefore, for a model to be realistic, it should incorporate some nonlinear structures.Note that some nonlinear processes can be approximated to be nearly linear by the log-transformation [11, 12] or time-scale transformation on the degradation data [13–15]. But many nonlinear degradation processes cannot be transformed in these ways. In addition, considering that each product possibly experiences different sources of variations during its operation, a random effect model is more appropriate to incorporate unit to unit variability in the degradation process. Wang [13] considered a random effect Wiener process model to fit bridge beam data and corresponding goodness-of-fit tests were carried out, but his article only considered the nonlinear process which could be linear by time-scale transformation. Therefore, a random effect nonlinear diffusion model without the use of the transformations mentioned above is maybe more appropriate for describing the degradation path.Moreover, the above studies consider only one performance characteristic (PC) or component failure mechanism. However, modern products usually have complex structure and more functions, which mean that multiple degradation mechanisms may be involved. For example, with the different purpose of lighting (i.e., illumination and color), the LED system may generate multiple degradation mechanism that leads to failure [16]. There is little literature to consider multivariate PCs. Related work can be found in [17–20]. However, these works use either independence assumption of the multiple PCs or multivariate normal distribution. From a practical point of view, these assumptions may not be realistic. Recently, to relax these assumptions, a copula function is used to describe the dependency of PCs. In [21, 22], Sari et al. answered the question as to how one could quantify the reliability of a system which has two or more degraded PCs, and these PCs can be described through a copula function.But in Sari’s work, they modeled the degradation data with generalized linear regression model with population average approach. Compared with stochastic process model, the regression model ignores the temporal uncertainty of a degradation process, which results in limiting its applications. Hence, based on Sari’s work, Pan et al. [23, 24] used the Wiener process to model the degradation failure mechanism and used the copula function to describe the dependence of multiple degradation failure mechanism. Similar to the Sari’s work, they also used the population average approach and did not take into account the random effect in their model.From the above analysis, it is clear that the reliability assessment of the products with two PCs has not been studied thoroughly. In this paper, we assume that a product has two PCs and the degradation path of each PC can be governed by a nonlinear diffusion process with random effects and that the dependency of the PCs can be described by the Franck copula function. The model in such a situation is very complicated and analytically intractable and becomes cumbersome from a computational viewpoint. For this reason, the Bayesian Markov Chain Monte Carlo (MCMC) method is used to obtain the unknown parameters [25].The MCMC method is a simulation technique in which the analytical posterior distribution is difficult to be computed. The MCMC techniques, including the Metropolis-Hastings (M-H) algorithm [26] and the Gibbs sampler [27], have become very common as methods for generating a sample from a complicated model. Recently, application of the Gibbs sampler to the estimation of parameters or some other vital properties about statistical models has become very popular.The rest of the paper is organized as follows. In Section2, some assumptions and copula basic are described. Then, the bivariate degradation model based on random effects nonlinear diffusion process is introduced in Section 3. In Section 4, the estimation of unknown parameters based on the MCMC algorithm is presented. A numerical example about an actual LED lamp data is given in Section 5. Finally, some conclusions are made in Section 6.


## 2. Some Basic Assumptions and Copula Brief Introduction

### 2.1. Basic Assumptions

To analyze the problem statistically, some assumptions are used for the reliability modeling in this paper. The details of each assumption are explained in the corresponding sections.(1)
The samples are independent, and no catastrophic failures occur during the degradation process.
(2)
The marginal degradation processes can be modeled as a nonlinear diffusion process with random effect.
(3)
For a specific product, the degradation measurements on the two PCs are observed at the same time (balanced data).
(4)
Suppose that the degradation path of a product is a decreasing function; then, a product is announced to be failed if any PC is less than its corresponding failure threshold for the first time. Here, the failure threshold vector is denoted byξ = (ξ
1, ξ
2).
(5)
The two PCs are dependent on each other, and the dependency can be characterized by a Frank copula function.


### 2.2. Copulas and Their Properties

#### 2.2.1. Definition and Base Properties

Copulas provide a very convenient way to model and measure the dependence among multiple PCs since they give the dependence structure which relates the known marginal distributions of each PC to their multivariate joint distribution. In order to see this, we first provide a brief introduction on copulas.A two dimensional copulaC
(
u
,
v
)
(
u
∈
[
0,1
]
,
v
∈
[
0,1
]
) is conventionally defined as a bivariate cumulative distribution function with uniform margins. A probabilistic way to define the copula is provided by the Sklar theorem [28].Theorem A (Sklar, 1959).
LetX and Y be random variables with continuous distributions F
(
x
) and G
(
y
), respectively, and let H
(
x
,
y
) be the two dimensional cumulative distribution function. Then, there exists a two dimensional copula C
(
·
,
·
) such that for all x
,
y
∈
(
-
∞
,
+
∞
),
(1)
H
(
x
,
y
)
=
C
(
F
(
x
)
,
G
(
y
)
)
.
IfF
(
·
) and G
(
·
) are continuous, then C
(
·
,
·
) is unique; conversely, if C
(
·
,
·
) is a copula and F
(
·
) and G
(
·
) are distribution functions, then the function H
(
·
,
·
) defined by (1) is a joint distribution function with margins F
(
·
) and G
(
·
).


#### 2.2.2. Measures of Association

We will consider here the standard dependence measures, Kendall’sτ and Spearman’s ρ. These measures are related to the copula since the latter is an expression of the stochastic relationship between u and v within the entire range of values the variables can take. It is not difficult to show that
(2)
ρ
=
12
∬
[
0,1
]
2
C
(
u
,
v
)
d
u
d
v
-
3
,
τ
=
3
∬
[
0,1
]
2
C
(
u
,
v
)
d
C
(
u
,
v
)
-
1
,
where the range of τ and ρ  can be shown to be [−1, 1]. For other properties of τ and ρ, see in [28].In this paper, the Frank copula is used to depict the dependence among multiple PCs that is given by(3)
C
(
u
,
v
)
=
-
1
α
ln
⁡
{
1
+
[
exp
⁡
⁡
(
-
α
u
)
-
1
]
[
exp
⁡
⁡
(
-
α
v
)
-
1
]
exp
⁡
⁡
(
-
α
)
-
1
}
,
where α is the Frank copula parameter and α
∈
(
-
∞
,
0
)
∪
(
0
,
+
∞
). And the relationship between Kendall’s τ and the Frank copula parameter α is given by
(4)
τ
=
1
+
4
D
1
(
α
)
-
1
α
,
where D
1
(
α
)
=
(
1
/
α
)
∫
0
α
(
t
/
(
e
t
-
1
)
)
d
t is a Debye function of the first kind (see in Nelson [28]).


## 2.1. Basic Assumptions

To analyze the problem statistically, some assumptions are used for the reliability modeling in this paper. The details of each assumption are explained in the corresponding sections.(1)
The samples are independent, and no catastrophic failures occur during the degradation process.
(2)
The marginal degradation processes can be modeled as a nonlinear diffusion process with random effect.
(3)
For a specific product, the degradation measurements on the two PCs are observed at the same time (balanced data).
(4)
Suppose that the degradation path of a product is a decreasing function; then, a product is announced to be failed if any PC is less than its corresponding failure threshold for the first time. Here, the failure threshold vector is denoted byξ = (ξ
1, ξ
2).
(5)
The two PCs are dependent on each other, and the dependency can be characterized by a Frank copula function.


## 2.2. Copulas and Their Properties

### 2.2.1. Definition and Base Properties

Copulas provide a very convenient way to model and measure the dependence among multiple PCs since they give the dependence structure which relates the known marginal distributions of each PC to their multivariate joint distribution. In order to see this, we first provide a brief introduction on copulas.A two dimensional copulaC
(
u
,
v
)
(
u
∈
[
0,1
]
,
v
∈
[
0,1
]
) is conventionally defined as a bivariate cumulative distribution function with uniform margins. A probabilistic way to define the copula is provided by the Sklar theorem [28].Theorem A (Sklar, 1959).
LetX and Y be random variables with continuous distributions F
(
x
) and G
(
y
), respectively, and let H
(
x
,
y
) be the two dimensional cumulative distribution function. Then, there exists a two dimensional copula C
(
·
,
·
) such that for all x
,
y
∈
(
-
∞
,
+
∞
),
(1)
H
(
x
,
y
)
=
C
(
F
(
x
)
,
G
(
y
)
)
.
IfF
(
·
) and G
(
·
) are continuous, then C
(
·
,
·
) is unique; conversely, if C
(
·
,
·
) is a copula and F
(
·
) and G
(
·
) are distribution functions, then the function H
(
·
,
·
) defined by (1) is a joint distribution function with margins F
(
·
) and G
(
·
).


### 2.2.2. Measures of Association

We will consider here the standard dependence measures, Kendall’sτ and Spearman’s ρ. These measures are related to the copula since the latter is an expression of the stochastic relationship between u and v within the entire range of values the variables can take. It is not difficult to show that
(2)
ρ
=
12
∬
[
0,1
]
2
C
(
u
,
v
)
d
u
d
v
-
3
,
τ
=
3
∬
[
0,1
]
2
C
(
u
,
v
)
d
C
(
u
,
v
)
-
1
,
where the range of τ and ρ  can be shown to be [−1, 1]. For other properties of τ and ρ, see in [28].In this paper, the Frank copula is used to depict the dependence among multiple PCs that is given by(3)
C
(
u
,
v
)
=
-
1
α
ln
⁡
{
1
+
[
exp
⁡
⁡
(
-
α
u
)
-
1
]
[
exp
⁡
⁡
(
-
α
v
)
-
1
]
exp
⁡
⁡
(
-
α
)
-
1
}
,
where α is the Frank copula parameter and α
∈
(
-
∞
,
0
)
∪
(
0
,
+
∞
). And the relationship between Kendall’s τ and the Frank copula parameter α is given by
(4)
τ
=
1
+
4
D
1
(
α
)
-
1
α
,
where D
1
(
α
)
=
(
1
/
α
)
∫
0
α
(
t
/
(
e
t
-
1
)
)
d
t is a Debye function of the first kind (see in Nelson [28]).


## 2.2.1. Definition and Base Properties

Copulas provide a very convenient way to model and measure the dependence among multiple PCs since they give the dependence structure which relates the known marginal distributions of each PC to their multivariate joint distribution. In order to see this, we first provide a brief introduction on copulas.A two dimensional copulaC
(
u
,
v
)
(
u
∈
[
0,1
]
,
v
∈
[
0,1
]
) is conventionally defined as a bivariate cumulative distribution function with uniform margins. A probabilistic way to define the copula is provided by the Sklar theorem [28].Theorem A (Sklar, 1959).
LetX and Y be random variables with continuous distributions F
(
x
) and G
(
y
), respectively, and let H
(
x
,
y
) be the two dimensional cumulative distribution function. Then, there exists a two dimensional copula C
(
·
,
·
) such that for all x
,
y
∈
(
-
∞
,
+
∞
),
(1)
H
(
x
,
y
)
=
C
(
F
(
x
)
,
G
(
y
)
)
.
IfF
(
·
) and G
(
·
) are continuous, then C
(
·
,
·
) is unique; conversely, if C
(
·
,
·
) is a copula and F
(
·
) and G
(
·
) are distribution functions, then the function H
(
·
,
·
) defined by (1) is a joint distribution function with margins F
(
·
) and G
(
·
).


## 2.2.2. Measures of Association

We will consider here the standard dependence measures, Kendall’sτ and Spearman’s ρ. These measures are related to the copula since the latter is an expression of the stochastic relationship between u and v within the entire range of values the variables can take. It is not difficult to show that
(2)
ρ
=
12
∬
[
0,1
]
2
C
(
u
,
v
)
d
u
d
v
-
3
,
τ
=
3
∬
[
0,1
]
2
C
(
u
,
v
)
d
C
(
u
,
v
)
-
1
,
where the range of τ and ρ  can be shown to be [−1, 1]. For other properties of τ and ρ, see in [28].In this paper, the Frank copula is used to depict the dependence among multiple PCs that is given by(3)
C
(
u
,
v
)
=
-
1
α
ln
⁡
{
1
+
[
exp
⁡
⁡
(
-
α
u
)
-
1
]
[
exp
⁡
⁡
(
-
α
v
)
-
1
]
exp
⁡
⁡
(
-
α
)
-
1
}
,
where α is the Frank copula parameter and α
∈
(
-
∞
,
0
)
∪
(
0
,
+
∞
). And the relationship between Kendall’s τ and the Frank copula parameter α is given by
(4)
τ
=
1
+
4
D
1
(
α
)
-
1
α
,
where D
1
(
α
)
=
(
1
/
α
)
∫
0
α
(
t
/
(
e
t
-
1
)
)
d
t is a Debye function of the first kind (see in Nelson [28]).


## 3. Reliability Model Based on Nonlinear Diffusion Process

### 3.1. Marginal Reliability Model Based on Nonlinear Diffusion Process

A well-adopted form for the Wiener process{
X
(
t
)
,
t
≥
0
} can be expressed as M
1:
(5)
X
(
t
)
=
μ
Λ
(
t
)
+
σ
B
B
(
Λ
(
t
)
)
,
where  μ is the drift degradation rate, σ
B is the diffusion coefficient, Λ
(
t
) is the transformed time scale, and B
(
·
) is the standard Brownian motion representing a time correlated structure.When theΛ
(
t
)
=
t, the degradation model reduces to the conventional Wiener model with linear drifted as M
2:
(6)
X
(
t
)
=
μ
t
+
σ
B
B
(
t
)
.As we known, with a time-scale transformation ofs
=
Λ
(
t
), the model M
1 reduces to the model M
2. That is,
(7)
Y
(
s
)
=
μ
s
+
σ
B
B
(
s
)
,
where Y
(
s
)
=
X
(
t
) and s
=
Λ
(
t
).From (5)–(7), the degradation models M
1 and M
2 are only considered the linear degradation process or the nonlinear process which can be linear by time-scale transformation. However, nonlinearity exists extensively in practice, and many nonlinear degradation processes cannot be properly linearized. Hence, a nonlinear degradation directly without the use of the transformations may be more appropriate for describing the degradation path, and this creates a unique and challenging problem to analyze the system reliability.In this paper, a diffusion process{
X
(
t
)
,
t
≥
0
} with a nonlinear drift can be written as M
3:
(8)
X
(
t
)
=
x
0
+
∫
0
t
ξ
(
t
;
η
)
d
t
+
σ
B
B
(
t
)
,
where x
0 denotes the initial state and is assumed to be zero, which is frequently used in literature, and ∫
0
t
ξ
(
t
;
η
)
d
t and σ
B are the drift part and diffusion coefficient, respectively. From the degradation model in (8), we can find that if ∫
0
t
ξ
(
t
;
η
)
d
t is a nonlinear function with an unknown parameter vector η, the model can capture the nonlinear characteristics of the degradation process. Moreover, it can contain various nonlinear processes by selecting the different forms of ξ
(
t
;
η
). Clearly, if ∫
0
t
ξ
(
t
;
η
)
d
t
=
μ
t, M
3 will reduce to the traditional Wiener process with the linear drift model M
2. In this paper, it is supposed that ∫
0
t
ξ
(
t
;
η
)
d
t
=
μ
Λ
(
t
,
b
), and the degradation model becomes M
4:
(9)
X
(
t
)
=
μ
Λ
(
t
)
+
σ
B
B
(
t
)
,
where Λ
(
t
)
=
Λ
(
t
,
b
) is a nonlinear function of time t.Suppose that a product has two PCs and each PC is governed by a nonlinear diffusion process. Then, thekth PC is defined by (9) as
(10)
X
k
(
t
)
=
μ
k
Λ
k
(
t
)
+
σ
B
k
B
(
t
)
,
where k
=
1,2 and Λ
k
(
t
)
=
Λ
(
t
,
b
k
).Letξ
k be the threshold value of thekth PC; the lifetime T
k of thekth PC is defined as
(11)
T
k
=
inf
⁡
{
t
∣
X
k
(
t
)
≤
ξ
k
}
.Theorem 1.
For the degradation process{
X
(
t
)
,
t
≥
0
}, if Λ
k
(
t
)
=
Λ
(
t
,
b
k
) is a continuously differentiable function of time t, then the probability density function (PDF) of the lifetime T
k can be obtained with an explicit form as
(12)
f
T
k
(
t
∣
μ
k
,
σ
B
k
2
)
=
ξ
k
2
π
σ
B
k
2
t
3
exp
⁡
(
-
ξ
k
-
μ
k
Λ
k
(
t
)
2
σ
B
k
2
t
)
+
exp
⁡
(
2
ξ
k
μ
k
Λ
k
(
t
)
σ
B
k
2
t
)
(
2
ξ
k
μ
k
σ
B
k
2
(
Λ
k
(
t
)
t
)
T
)
Φ
×
(
-
ξ
k
-
μ
k
Λ
k
(
t
)
σ
B
k
2
t
)
and the corresponding cumulative distribution function (CDF) is
(13)
F
T
k
(
t
∣
μ
k
,
σ
B
k
2
)
=
1
-
Φ
(
ξ
k
-
μ
k
Λ
k
(
t
)
σ
B
k
2
t
)
+
exp
⁡
(
2
ξ
k
μ
k
Λ
k
(
t
)
σ
B
k
2
t
)
Φ
(
-
ξ
k
-
μ
k
Λ
k
(
t
)
σ
B
k
2
t
)
,
where Φ
(
·
) is the CDF of the standard normal distribution and “
T
” denotes first derivative.The proof is given in AppendixB.Then, the marginal reliability at timet can be expressed as
(14)
R
k
(
t
)
=
1
-
F
T
k
(
t
∣
μ
k
,
σ
B
k
2
)
.In most degradation applications, considering that each item possibly experiences different sources of variations during its operation, for a degradation model to be realistic, it is more appropriate to incorporate unit to unit variability in the degradation process. The degradations of such products can be described by random effect models, where random effects capture the unit to unit differences. In this paper, we treat the drift parameterμ  as a random effect representing between item variation, and the diffusion σ
B is a fixed effect that is common to all items. For simplicity, we assume that μ and B
(
t
) are independent and that μ follows N
(
μ
β
,
σ
β
2
). The ideas of random effects and Gaussian assumptions are widely used in degradation modeling [10, 12, 13]. Then, we can get the random effect nonlinear diffusion process model as follows, M
5:
(15)
X
(
t
)
=
μ
Λ
(
t
)
+
σ
B
B
(
t
)
μ
~
N
(
μ
β
,
σ
β
2
)
.This paper is focused on the nonlinear random effect diffusion degradation process modelM
5. The advantages of the degradation model in (15) allow us to takeunit to unit variation,nonlinear degradation characteristics, andtime correlated structure into consideration simultaneously.Similarly, suppose that a product has two PCs and each PC is governed by a random effects nonlinear diffusion process. Thekth PC is defined by (15) as
(16)
X
k
(
t
)
=
μ
k
Λ
k
(
t
)
+
σ
B
k
B
(
t
)
μ
k
~
N
(
μ
β
k
,
σ
β
k
2
)
.Theorem 2.
For the degradation process as defined by (16), if Λ
k
(
t
)
=
Λ
(
t
,
b
k
) is a continuously differentiable function of time t, then the CDF of the lifetime T
k can be obtained with an explicit form as
(17)
F
T
k
(
t
∣
σ
B
k
2
)
=
1
-
Φ
(
ξ
k
-
μ
β
k
Λ
k
(
t
)
σ
B
k
2
t
+
σ
β
k
2
Λ
k
2
(
t
)
)
+
exp
⁡
(
2
ξ
k
μ
k
Λ
k
(
t
)
σ
B
k
2
t
+
2
ξ
k
2
σ
β
k
2
Λ
k
2
(
t
)
σ
B
k
4
t
2
)
×
Φ
(
-
(
ξ
k
+
μ
β
k
Λ
k
(
t
)
)
σ
B
k
2
t
+
2
ξ
k
Λ
k
2
(
t
)
σ
β
k
2
σ
B
k
2
t
σ
B
k
2
t
+
Λ
k
2
(
t
)
σ
β
k
2
)
,
where Φ
(
·
) is the CDF of the standard normal distribution.The proof is given in AppendixD.


### 3.2. Degradation Model Based on Bivariate Degradation Data

Suppose that product has two PCsX
1(t) and X
2(t) at observation time point t and given the failure threshold ξ
k of thekth PC. From the basic assumptions, the marginal reliability at time t can be expressed as
(18)
R
k
(
t
)
=
Pr
(
X
k
(
t
)
>
ξ
k
)
=
1
-
F
T
k
(
t
∣
σ
B
k
2
)
.If the two PCs are dependent, the copula method is utilized to establish the dependent structure among the two PCs. LetR
k
(
t
) be the marginal reliability at time t and C
(
·
,
·
) be the joint copula of the marginal degradation distributions. Similarly to the method in [22, 24], the system reliability at time t can be expressed as
(19)
R
(
t
)
=
C
(
R
1
(
t
)
,
R
2
(
t
)
)
.If the product has two PCs linked by bivariate Frank copula given in (3), then we can obtain the system reliability function as
(20)
R
(
t
)
=
-
1
α
ln
⁡
{
1
+
[
exp
⁡
(
-
α
R
1
(
t
)
)
-
1
]
[
exp
⁡
(
-
α
R
2
(
t
)
)
-
1
]
exp
⁡
⁡
(
-
α
)
-
1
}
.


## 3.1. Marginal Reliability Model Based on Nonlinear Diffusion Process

A well-adopted form for the Wiener process{
X
(
t
)
,
t
≥
0
} can be expressed as M
1:
(5)
X
(
t
)
=
μ
Λ
(
t
)
+
σ
B
B
(
Λ
(
t
)
)
,
where  μ is the drift degradation rate, σ
B is the diffusion coefficient, Λ
(
t
) is the transformed time scale, and B
(
·
) is the standard Brownian motion representing a time correlated structure.When theΛ
(
t
)
=
t, the degradation model reduces to the conventional Wiener model with linear drifted as M
2:
(6)
X
(
t
)
=
μ
t
+
σ
B
B
(
t
)
.As we known, with a time-scale transformation ofs
=
Λ
(
t
), the model M
1 reduces to the model M
2. That is,
(7)
Y
(
s
)
=
μ
s
+
σ
B
B
(
s
)
,
where Y
(
s
)
=
X
(
t
) and s
=
Λ
(
t
).From (5)–(7), the degradation models M
1 and M
2 are only considered the linear degradation process or the nonlinear process which can be linear by time-scale transformation. However, nonlinearity exists extensively in practice, and many nonlinear degradation processes cannot be properly linearized. Hence, a nonlinear degradation directly without the use of the transformations may be more appropriate for describing the degradation path, and this creates a unique and challenging problem to analyze the system reliability.In this paper, a diffusion process{
X
(
t
)
,
t
≥
0
} with a nonlinear drift can be written as M
3:
(8)
X
(
t
)
=
x
0
+
∫
0
t
ξ
(
t
;
η
)
d
t
+
σ
B
B
(
t
)
,
where x
0 denotes the initial state and is assumed to be zero, which is frequently used in literature, and ∫
0
t
ξ
(
t
;
η
)
d
t and σ
B are the drift part and diffusion coefficient, respectively. From the degradation model in (8), we can find that if ∫
0
t
ξ
(
t
;
η
)
d
t is a nonlinear function with an unknown parameter vector η, the model can capture the nonlinear characteristics of the degradation process. Moreover, it can contain various nonlinear processes by selecting the different forms of ξ
(
t
;
η
). Clearly, if ∫
0
t
ξ
(
t
;
η
)
d
t
=
μ
t, M
3 will reduce to the traditional Wiener process with the linear drift model M
2. In this paper, it is supposed that ∫
0
t
ξ
(
t
;
η
)
d
t
=
μ
Λ
(
t
,
b
), and the degradation model becomes M
4:
(9)
X
(
t
)
=
μ
Λ
(
t
)
+
σ
B
B
(
t
)
,
where Λ
(
t
)
=
Λ
(
t
,
b
) is a nonlinear function of time t.Suppose that a product has two PCs and each PC is governed by a nonlinear diffusion process. Then, thekth PC is defined by (9) as
(10)
X
k
(
t
)
=
μ
k
Λ
k
(
t
)
+
σ
B
k
B
(
t
)
,
where k
=
1,2 and Λ
k
(
t
)
=
Λ
(
t
,
b
k
).Letξ
k be the threshold value of thekth PC; the lifetime T
k of thekth PC is defined as
(11)
T
k
=
inf
⁡
{
t
∣
X
k
(
t
)
≤
ξ
k
}
.Theorem 1.
For the degradation process{
X
(
t
)
,
t
≥
0
}, if Λ
k
(
t
)
=
Λ
(
t
,
b
k
) is a continuously differentiable function of time t, then the probability density function (PDF) of the lifetime T
k can be obtained with an explicit form as
(12)
f
T
k
(
t
∣
μ
k
,
σ
B
k
2
)
=
ξ
k
2
π
σ
B
k
2
t
3
exp
⁡
(
-
ξ
k
-
μ
k
Λ
k
(
t
)
2
σ
B
k
2
t
)
+
exp
⁡
(
2
ξ
k
μ
k
Λ
k
(
t
)
σ
B
k
2
t
)
(
2
ξ
k
μ
k
σ
B
k
2
(
Λ
k
(
t
)
t
)
T
)
Φ
×
(
-
ξ
k
-
μ
k
Λ
k
(
t
)
σ
B
k
2
t
)
and the corresponding cumulative distribution function (CDF) is
(13)
F
T
k
(
t
∣
μ
k
,
σ
B
k
2
)
=
1
-
Φ
(
ξ
k
-
μ
k
Λ
k
(
t
)
σ
B
k
2
t
)
+
exp
⁡
(
2
ξ
k
μ
k
Λ
k
(
t
)
σ
B
k
2
t
)
Φ
(
-
ξ
k
-
μ
k
Λ
k
(
t
)
σ
B
k
2
t
)
,
where Φ
(
·
) is the CDF of the standard normal distribution and “
T
” denotes first derivative.The proof is given in AppendixB.Then, the marginal reliability at timet can be expressed as
(14)
R
k
(
t
)
=
1
-
F
T
k
(
t
∣
μ
k
,
σ
B
k
2
)
.In most degradation applications, considering that each item possibly experiences different sources of variations during its operation, for a degradation model to be realistic, it is more appropriate to incorporate unit to unit variability in the degradation process. The degradations of such products can be described by random effect models, where random effects capture the unit to unit differences. In this paper, we treat the drift parameterμ  as a random effect representing between item variation, and the diffusion σ
B is a fixed effect that is common to all items. For simplicity, we assume that μ and B
(
t
) are independent and that μ follows N
(
μ
β
,
σ
β
2
). The ideas of random effects and Gaussian assumptions are widely used in degradation modeling [10, 12, 13]. Then, we can get the random effect nonlinear diffusion process model as follows, M
5:
(15)
X
(
t
)
=
μ
Λ
(
t
)
+
σ
B
B
(
t
)
μ
~
N
(
μ
β
,
σ
β
2
)
.This paper is focused on the nonlinear random effect diffusion degradation process modelM
5. The advantages of the degradation model in (15) allow us to takeunit to unit variation,nonlinear degradation characteristics, andtime correlated structure into consideration simultaneously.Similarly, suppose that a product has two PCs and each PC is governed by a random effects nonlinear diffusion process. Thekth PC is defined by (15) as
(16)
X
k
(
t
)
=
μ
k
Λ
k
(
t
)
+
σ
B
k
B
(
t
)
μ
k
~
N
(
μ
β
k
,
σ
β
k
2
)
.Theorem 2.
For the degradation process as defined by (16), if Λ
k
(
t
)
=
Λ
(
t
,
b
k
) is a continuously differentiable function of time t, then the CDF of the lifetime T
k can be obtained with an explicit form as
(17)
F
T
k
(
t
∣
σ
B
k
2
)
=
1
-
Φ
(
ξ
k
-
μ
β
k
Λ
k
(
t
)
σ
B
k
2
t
+
σ
β
k
2
Λ
k
2
(
t
)
)
+
exp
⁡
(
2
ξ
k
μ
k
Λ
k
(
t
)
σ
B
k
2
t
+
2
ξ
k
2
σ
β
k
2
Λ
k
2
(
t
)
σ
B
k
4
t
2
)
×
Φ
(
-
(
ξ
k
+
μ
β
k
Λ
k
(
t
)
)
σ
B
k
2
t
+
2
ξ
k
Λ
k
2
(
t
)
σ
β
k
2
σ
B
k
2
t
σ
B
k
2
t
+
Λ
k
2
(
t
)
σ
β
k
2
)
,
where Φ
(
·
) is the CDF of the standard normal distribution.The proof is given in AppendixD.


## 3.2. Degradation Model Based on Bivariate Degradation Data

Suppose that product has two PCsX
1(t) and X
2(t) at observation time point t and given the failure threshold ξ
k of thekth PC. From the basic assumptions, the marginal reliability at time t can be expressed as
(18)
R
k
(
t
)
=
Pr
(
X
k
(
t
)
>
ξ
k
)
=
1
-
F
T
k
(
t
∣
σ
B
k
2
)
.If the two PCs are dependent, the copula method is utilized to establish the dependent structure among the two PCs. LetR
k
(
t
) be the marginal reliability at time t and C
(
·
,
·
) be the joint copula of the marginal degradation distributions. Similarly to the method in [22, 24], the system reliability at time t can be expressed as
(19)
R
(
t
)
=
C
(
R
1
(
t
)
,
R
2
(
t
)
)
.If the product has two PCs linked by bivariate Frank copula given in (3), then we can obtain the system reliability function as
(20)
R
(
t
)
=
-
1
α
ln
⁡
{
1
+
[
exp
⁡
(
-
α
R
1
(
t
)
)
-
1
]
[
exp
⁡
(
-
α
R
2
(
t
)
)
-
1
]
exp
⁡
⁡
(
-
α
)
-
1
}
.


## 4. Statistical Inferential Methods for Unknown Parameters

Here, we discuss the estimation of parameters required to implement the reliability function in (20). The unknown parameters are
(21)
θ
=
(
μ
β
1
,
σ
β
1
,
σ
B
1
,
b
1
,
μ
β
2
,
σ
β
2
,
σ
B
2
,
b
2
,
α
)
.From (17), (18), and (20), we know that the model not only has nine parameters but also is very complicated from a computational viewpoint. For this reason, the MCMC with the Gibbs sampling techniques is employed in this study to estimate model parameters. Let π
(
θ
j
∣
θ
-
j
,
Y
) denote the full conditional posterior distribution of θ
j, where θ
-
j
=
(
θ
1
,
…
,
θ
j
-
1
,
θ
j
+
1
,
…
,
θ
n
) and Y is the observed data.Then, the Gibbs sampling algorithm can be summarized as follows.Step 1.
Initializeθ
(
0
)
=
(
θ
1
(
0
)
,
θ
2
(
0
)
,
…
,
θ
n
(
0
)
).Step 2.
Setθ
=
θ
(
t
-
1
), and generate θ
(
t
).Generateθ
1
(
t
) from π
1
*
(
θ
1
∣
θ
2
(
t
-
1
)
,
θ
3
(
t
-
1
)
,
…
,
θ
n
(
t
-
1
)
,
Y
)
,
…
Generateθ
j
(
t
) from π
j
*
(
θ
j
∣
θ
1
(
t
-
1
)
,
…
,
θ
j
-
1
(
t
-
1
)
,
θ
j
+
1
(
t
-
1
)
,
…
,
θ
n
(
t
-
1
)
,
Y
)
,
…
Generateθ
n
(
t
) from π
n
*
(
θ
n
∣
θ
1
(
t
-
1
)
,
θ
2
(
t
-
1
)
,
…
,
θ
n
-
1
(
t
-
1
)
,
Y
).Step 3.
Sett
=
t
+
1, and repeat Step 2, t
=
1,2
,
…
,
N
1.Step 4.
Estimate desired features based on the simulate sampleθ
(
m
)
,
θ
(
m
+
1
)
,
…
,
θ
(
N
1
), where m denotes the burn-in number. For example, we can use the mean of θ
j
(
m
)
,
θ
j
(
m
+
1
)
,
…
,
θ
j
(
N
1
) to be the estimator of θ
j.
We can use the Bayesian software package WinBUGS (see in [25]) to carry out the Gibbs sampling, and then we can obtain the estimator of the model parameters.


## 5. Numerical Example: Application to a LED Lighting System

In recent years, LED has attracted increasing interest in the field of lighting systems owing to its high efficiency, low energy, and long lifetime. Due to its different functions (i.e., illumination and color), there may be multiple failure mechanisms about the LED lighting system. The actual LED lamps experiment data are taken from Chaluvadi’s Ph.D. thesis [29], where the LED lamps were put into test under an accelerated current of 40 mA. In the original data, 12 samples are tested for LED data and the measurements are taken at the same measurement times. The measured frequency of its light intensity is 50 hours. For demonstrating the bivariate degradation model, similarly to [22, 24], we choose 12 samples and the data will be treated as if half of it is the first PC and the other half represents the second PC. That is to say, in the LED lighting system, the performance demand of the first group LED lamps is for lighting, and the performance demand of the second group LED lamps is for chromatic change. The used data are the data measured only until 250 hours. The LED lighting system is considered failed if one of the two LED datasets is less than 50. Table 1 lists the LED data. Figure 1 shows the cumulative degradation of the LED data.Table 1
Light intensity degradation data of 12 LEDs.
Unit
Inspection time (hour)
0
50
100
150
200
250
PC1
1
100
86.6
78.7
76.0
71.6
68.0
2
100
82.1
71.4
65.4
61.7
58.0
3
100
82.7
70.3
64.0
61.3
59.3
4
100
79.8
68.3
62.3
60.0
59.0
5
100
75.1
66.7
62.8
59.0
54.0
6
100
83.7
74.0
67.4
63.0
61.3
PC2
1
100
73.0
65.0
60.7
58.3
58.0
2
100
86.2
67.6
62.7
60.0
59.7
3
100
81.2
65.0
60.6
59.3
57.3
4
100
66.8
63.3
59.3
57.3
56.5
5
100
66.1
64.2
59.4
58.0
55.3
6
100
76.5
61.7
61.3
59.7
56.0The degradation paths of the LEDs under an accelerated current of 40 mA.
(a)
(b)


### 5.1. Estimation of Unknown Parameters

Now, we use the LED data to illustrate the proposed model and method in this paper. Firstly, we use the random effect nonlinear diffusion process model to fit the degradation data of PC1 and PC2, respectively. Tseng et al. [9] and Noortwijk [30] suggested that Λ
(
t
)
=
Λ
(
t
,
b
)
=
t
b is appropriate for the LED degradation modeling, and this form is also adopted in this paper. Considering the dependency between the two PCs, by using the MCMC method in Section 4, we generate 50,000 samples. A burn-in of 10,000 samples is used, with an additional 40,000 Gibbs samples that are used to estimate parameters. Table 2 tabulates posterior summaries including parameters posterior mean, standard error, Monte Carlo error, and 95% HPD interval.Table 2
Parameter estimation results considering the dependency.
Parameter
Mean
Standard error
MC error
95% HPD interval
μ
β
1
−2.850
0.4472
0.0389
(−3.810, −2.072)
σ
β
1
0.2788
0.2298
0.0134
(0.0299, 0.8273)
σ
B
1
0.3557
0.0550
0.0016
(0.2660, 0.4792)
b
1
0.4850
0.0305
0.0027
(0.4222, 0.5413)
μ
β
2
−7.357
1.7390
0.1635
(−11.100, −4.282)
σ
β
2
0.5801
0.7412
0.0494
(0.0309, 2.444)
σ
B
2
0.7259
0.1030
0.0024
(0.5546, 0.9584)
b
2
0.3291
0.0535
0.0049
(0.2274, 0.4367)
α
3.5550
2.0620
0.1152
(0.0000, 7.4160)Meanwhile, the estimates of the unknown parameters without considering the dependency between the two PCs are presented in Table3.Table 3
Parameter estimation results without considering the dependency.
Parameter
Mean
Standard errors
MC error
95% HPD interval
μ
β
1
−3.005
0.5444
0.0214
(−4.166, −2.082)
σ
β
1
0.3082
0.2619
0.0085
(0.0317, 0.9514)
σ
B
1
0.3500
0.0547
0.0010
(0.2611, 0.4750)
b
1
0.4719
0.0355
0.0014
(0.4004, 0.5389)
μ
β
2
−7.894
2.2860
0.1033
(−12.65, −4.198)
σ
β
2
0.7973
0.9306
0.0356
(0.0331, 3.228)
σ
B
2
0.7228
0.1056
0.0015
(0.5486, 0.9623)
b
2
0.3088
0.0628
0.0029
(0.1941, 0.4275)To demonstrate the goodness of fit, the mean degradation paths estimated from the model and from the sample average are compared. Based on (16), it is easy to see that the estimated degradation path based on the random effect nonlinear diffusion process is given by
(22)
E
[
X
k
(
t
)
]
=
μ
^
β
k
Λ
^
k
(
t
)
=
μ
^
β
k
t
b
^
k
.Intuitively, if the nonlinear diffusion process model is correct, then the estimated mean degradation paths based on (22) should agree with the sample average. Therefore, we plot the estimated mean degradation paths based on (22) in conjunction with the sample average in Figure 2. The two curves tally quite well, implying the goodness of fit of the nonlinear diffusion process model as shown in Figure 2.Estimated mean paths based on the sample average and on the fit model.
(a)
(b)


### 5.2. Reliability Assessment

Based on the estimated results of the unknown parameters, the marginal reliability curves of PC1 and PC2 considering the dependency and independency are presented in Figures3, 4, 5, and 6.Figure 3
Marginal reliability curves of PC1.Figure 4
Marginal reliability curves of PC2.Figure 5
Marginal reliability curves of system under independent cases.Figure 6
Marginal reliability curves of system under dependent cases.It can be concluded from the Figures3-4 that either PC1 or PC2 has a relatively higher reliability under without considering the dependency case than considering the dependency case. In other words, the incorrect independent assumption may overrate the reliability. In addition, take the dependent case as an example in Figure 6, the PC1 has a relatively higher reliability than PC2 before about 400 hours; therefore, the system reliability may tend more to the marginal reliability of PC2. Similarly, after about 400 hours, the system reliability may tend more to the marginal reliability of PC1.Furthermore, the system reliability curves under the independent and dependent (with Frank copula assumption) case are plotted in Figure7. Comparing Figure 7 with Figures 3–6, we can find that there has been larger difference between the marginal reliability and the system reliability. Furthermore, from Figure 7, we can also find that there are some differences between the dependent and independent cases. That is to say that ignoring the dependence between PCs may result in different reliability conclusion. Therefore, it is necessary to analyze the possibility of the failure mechanisms dependency and perform the dependent reliability analysis.Figure 7
Reliability curves of system based on dependent and independent cases.


## 5.1. Estimation of Unknown Parameters

Now, we use the LED data to illustrate the proposed model and method in this paper. Firstly, we use the random effect nonlinear diffusion process model to fit the degradation data of PC1 and PC2, respectively. Tseng et al. [9] and Noortwijk [30] suggested that Λ
(
t
)
=
Λ
(
t
,
b
)
=
t
b is appropriate for the LED degradation modeling, and this form is also adopted in this paper. Considering the dependency between the two PCs, by using the MCMC method in Section 4, we generate 50,000 samples. A burn-in of 10,000 samples is used, with an additional 40,000 Gibbs samples that are used to estimate parameters. Table 2 tabulates posterior summaries including parameters posterior mean, standard error, Monte Carlo error, and 95% HPD interval.Table 2
Parameter estimation results considering the dependency.
Parameter
Mean
Standard error
MC error
95% HPD interval
μ
β
1
−2.850
0.4472
0.0389
(−3.810, −2.072)
σ
β
1
0.2788
0.2298
0.0134
(0.0299, 0.8273)
σ
B
1
0.3557
0.0550
0.0016
(0.2660, 0.4792)
b
1
0.4850
0.0305
0.0027
(0.4222, 0.5413)
μ
β
2
−7.357
1.7390
0.1635
(−11.100, −4.282)
σ
β
2
0.5801
0.7412
0.0494
(0.0309, 2.444)
σ
B
2
0.7259
0.1030
0.0024
(0.5546, 0.9584)
b
2
0.3291
0.0535
0.0049
(0.2274, 0.4367)
α
3.5550
2.0620
0.1152
(0.0000, 7.4160)Meanwhile, the estimates of the unknown parameters without considering the dependency between the two PCs are presented in Table3.Table 3
Parameter estimation results without considering the dependency.
Parameter
Mean
Standard errors
MC error
95% HPD interval
μ
β
1
−3.005
0.5444
0.0214
(−4.166, −2.082)
σ
β
1
0.3082
0.2619
0.0085
(0.0317, 0.9514)
σ
B
1
0.3500
0.0547
0.0010
(0.2611, 0.4750)
b
1
0.4719
0.0355
0.0014
(0.4004, 0.5389)
μ
β
2
−7.894
2.2860
0.1033
(−12.65, −4.198)
σ
β
2
0.7973
0.9306
0.0356
(0.0331, 3.228)
σ
B
2
0.7228
0.1056
0.0015
(0.5486, 0.9623)
b
2
0.3088
0.0628
0.0029
(0.1941, 0.4275)To demonstrate the goodness of fit, the mean degradation paths estimated from the model and from the sample average are compared. Based on (16), it is easy to see that the estimated degradation path based on the random effect nonlinear diffusion process is given by
(22)
E
[
X
k
(
t
)
]
=
μ
^
β
k
Λ
^
k
(
t
)
=
μ
^
β
k
t
b
^
k
.Intuitively, if the nonlinear diffusion process model is correct, then the estimated mean degradation paths based on (22) should agree with the sample average. Therefore, we plot the estimated mean degradation paths based on (22) in conjunction with the sample average in Figure 2. The two curves tally quite well, implying the goodness of fit of the nonlinear diffusion process model as shown in Figure 2.Estimated mean paths based on the sample average and on the fit model.
(a)
(b)


## 5.2. Reliability Assessment

Based on the estimated results of the unknown parameters, the marginal reliability curves of PC1 and PC2 considering the dependency and independency are presented in Figures3, 4, 5, and 6.Figure 3
Marginal reliability curves of PC1.Figure 4
Marginal reliability curves of PC2.Figure 5
Marginal reliability curves of system under independent cases.Figure 6
Marginal reliability curves of system under dependent cases.It can be concluded from the Figures3-4 that either PC1 or PC2 has a relatively higher reliability under without considering the dependency case than considering the dependency case. In other words, the incorrect independent assumption may overrate the reliability. In addition, take the dependent case as an example in Figure 6, the PC1 has a relatively higher reliability than PC2 before about 400 hours; therefore, the system reliability may tend more to the marginal reliability of PC2. Similarly, after about 400 hours, the system reliability may tend more to the marginal reliability of PC1.Furthermore, the system reliability curves under the independent and dependent (with Frank copula assumption) case are plotted in Figure7. Comparing Figure 7 with Figures 3–6, we can find that there has been larger difference between the marginal reliability and the system reliability. Furthermore, from Figure 7, we can also find that there are some differences between the dependent and independent cases. That is to say that ignoring the dependence between PCs may result in different reliability conclusion. Therefore, it is necessary to analyze the possibility of the failure mechanisms dependency and perform the dependent reliability analysis.Figure 7
Reliability curves of system based on dependent and independent cases.


## 6. Conclusion

In this paper, we establish a reliability assessment model with two dependent PCs, and each PC is governed by a random effect nonlinear diffusion process. The nonlinear diffusion model is very useful by checking the goodness of fit of the real data. We suppose that the two PCs are dependent and the dependency is described by a Frank copula function. The Bayesian MCMC method is used to obtain the unknown parameters. From the numerical example of the LED lighting system, we know that ignoring the dependence between PCs may result in different reliability conclusion.In this work, we have only considered the case that products have two marginal distributions with a nonlinear diffusion process. In the future research, it can be extended to consider multiple PCs by using multidimensional copula function. It also can be extended to the Gamma process, Geometric Brownian motion, and so on. In addition, we can use the different copula function (such as Gaussian copula and Gumbel copula) to describe the relationship. Moreover, from the practice point of view, how to make effective maintenance decisions for the products with PCs based on the proposed estimation results is necessary to be studied.



---

*Source: 510929-2014-03-30.xml*