Reliability inference and remaining useful life prediction for the doubly accelerated degradation model based on Wiener process

1.   Introduction

Accelerated life test (ALT) technology has been widely used in products' reliability analysis. However, even under ALT, little or no failure data can be acquired in a reasonable short life testing time for some high-reliability products. So, it is difficult to evaluate the reliability of those high-reliability products. In such a case, if there exists a quality characteristic (QC) related to reliability which degrades over time, an alternative approach is to collect the degradation data at the higher stress levels and then extrapolate the lifetime information and reliability metrics of products at the normal use stress level. Such an experiment is called an accelerated degradation test (ADT) ^[1,2]. For high-reliability products, the traditional ALT methods can not meet the requirement of product reliability evaluation, and the ADT technique provides another way to solve it ^[3].

As ADTs can greatly shorten the testing time, we often use them to quickly obtain more degradation data and make reliability analyses for products. In the past decades, the research of ADT models has become more and more popular ^[4,5,6]. In the literature, ADTs are broadly classified into constant-stress accelerated degradation test (CSADT), step-stress accelerated degradation test (SSADT) and progressive-step accelerated degradation test (PSADT), according to different stress loading methods. Among them, CSADT is the most popular ADT in real applications. In a CSADT, the testing units are divided into several groups, and each group of units is exposed to a certain severe stress condition to test and collect the degradation data.

In practice, the degradation path of products' QC over time is often described as stochastic $\{X(t); t\geq 0\}$ to account for inherent randomness. Based on the assumption of additive accumulation of degradation, three kinds of degradation processes have been well exploited, i.e., the Gamma process ^[7,8,9,10], the Wiener process ^[6,11,12,13] and the inverse Gaussian process ^[3,14,15,16]. In most cases, the degradation paths of test units are monotonic, so the Gamma process and inverse Gaussian process are often used to model degradation data. However, in a few cases, the degradation path is not monotonic. A distinct feature of the Wiener process is that its sample path is not necessarily monotone, which might be meaningful in some degradation applications. In this situation, the Wiener degradation model may be a good choice.

In view of this, the Wiener process, as an important stochastic process, has been widely studied in degradation data analysis. Pan and Balakrishnan ^[17] discussed the multiple-step SSADT models based on the Wiener and Gamma processes, and they used the Bayesian Markov Chain Monte Carlo method to obtain the maximum likelihood estimates (MLEs) for such analytically intractable models and presented some computational results obtained from their implementation. Motivated by the observation that a unit with a higher degradation rate has a more volatile degradation path, Ye and Chen ^[18] proposed a new class of random effects models for the Wiener process and discussed the statistical inference of the model. Wang et al. ^[13] mainly studied the accurate reliability inference for the Wiener degradation model with random drift parameter and developed an exact test method to test whether there exists population heterogeneity. Other research about the Wiener degradation model with random effects can be found in ^[11,19,20]. Guan et al. ^[21] used the Bayesian method to study the Wiener ADT model. Wang et al. ^[22] discussed the reliability analysis for accelerated degradation data based on the Wiener process with random effects. By constructing pivotal quantities (PQs), they developed generalized confidence intervals (GCIs) for model parameters and some quantities, and the generalized prediction interval (GPI) for future degradation amount at designed stress level was also developed. Hong et al. ^[23] developed interval estimation procedures for the Wiener degradation process with fixed-effects and mixed-effects using the generalized pivotal quantity (GPQ) method. Pan et al. ^[24] studied a reliability estimation approach via Wiener degradation model with measurement errors. Jiang et al. ^[12] proposed a Wiener CSADT model and obtained the exact confidence intervals (ECIs) of model parameters. In addition, they proposed a new optimization criterion from the perspective of degradation prediction, and provided an optimal experimental design scheme.

To construct an ADT model, first is to reasonably determine the relations between model parameters and accelerated stress. For the Wiener ADT model, there are three possible parameter-stress relations. The first is that the drift parameter of the Wiener process depends on the stress level while the diffusion parameter is independent of the stress level ^{[12,20,21,22,25]}. The second is that the diffusion parameter depends on the stress level while the drift parameter is independent of it. In the second case, however, the stress level does not have any effect on the degradation speed while the degradation volatility increases with the stress level. This case is not common in reality. The third is that both the drift and diffusion parameters are affected by the stress level ^[26]. Clearly, the first two relations are special cases of the third one. In addition, due to the complexity of the Wiener ADT model with the third parameter-stress relation, few scholars studied such models, especially for the aspect of interval estimation for reliability metrics. In this paper, we consider a Wiener doubly ADT model in which the parameter-stress relation follows the third one. The relations between parameters and accelerated stress are given by $\mu = a+b\xi, \sigma^2 = \exp(c+d\xi)$, and $\xi$ is the accelerated stress. We mainly focus on the interval estimations of model parameters and some reliability metrics based on the Wiener doubly ADT model.

The remainder of the paper is organized as follows. In Section 2, the general framework of the Wiener doubly ADT model is outlined. In Section 3, the point estimates of model parameters are derived. In Section 4, we mainly study the interval estimation of model parameters and some reliability metrics. In Section 5, a simulation study is conducted to evaluate the performance of the proposed GCIs/GPIs. In Section 6, an example is provided to illustrate the proposed model and GPQ method. Finally, we provide some final conclusions in Section 7.

2.   Model assumptions and data description

Suppose that one accelerated stress $\xi$ is used in a CSADT, and under the stress $\xi$ the degradation path $\{X(t), t\geq0\}$ of testing unit follows a Wiener process given by

where $B(t)$ is a standard Brownian motion, $\mu$ is the drift parameter, and $\sigma^2$ is the diffusion parameter. The unit's lifetime $T$ under the stress $\xi$ is defined as the first-passage-time of $X(t)$ to a pre-specified threshold $L$. It is well known that $T$ follows $IG(L/\mu, L^2/\sigma^2)$ distribution with cumulative distribution function (CDF)

where $\Phi(\cdot)$ is the CDF of a standard normal distribution.

For this study, statistical inference of the Wiener doubly accelerated degradation test model is usually based on the following assumptions:

(A1) The CSADT is conducted using a single stress, which has $K$ levels: $\xi_1 < \xi_2 < \cdots < \xi_K$. $\xi_0$ and $\xi_K$ are the normal use stress level and the highest stress level used in the ADT.

(A2) For each stress level $\xi_i$, the degradation path of a test unit can be described as a Wiener process $X_{i}(t)$,

(A3) Both the drift and the diffusion parameters are affected by the stress levels through the parameter-stress relationships

where, $a, b, c$ and $d$ are unknown parameters. The degradation rate and the diffusion parameter under normal use condition can be obtained by $\mu_0 = a+b\xi_0, \sigma_0^2 = \exp(c+d \xi_0)$, respectively.

Suppose that $n_i$ units are tested under the stress level $\xi_i$, and $r_{i, j}$ is the number of measurements for the $j$th test unit under the stress level $\xi_i$. The degradation characteristics of the $j$th test unit are measured at the times ${\bf t}_{i, j} = \{t_{i, j, k}; k = 0, 1, \ldots, r_{i, j}\}$ under the stress level $\xi_i$. Moreover, let ${\bf X}_{i, j} = \{X_{i, j}(t_{i, j, 0}), X_{i, j}(t_{i, j, 1}), \ldots, X_{i, j}(t_{i, j, r_{i, j}})\}$ denote the observed degradation characteristics of the $j$th test unit under the stress level $\xi_i$. Define $\Delta X_{i, j, k} = X_{i, j}(t_{i, j, k})-X_{i, j}(t_{i, j, k-1})$, and $\Delta t_{i, j, k} = t_{i, j, k}-t_{i, j, k-1}$, where $t_{i, j, 0} = 0$. The data collected from stress level $\xi_i$ is $\mathbb{D}_i = \{({\bf t}_{i, j}, {\bf X}_{i, j}); j = 1, 2, \ldots, n_i\}$, and the data from the whole double ADT is $\mathbb{D} = \bigcup_{i = 1}^{K}\mathbb{D}_i$. Let $N = \sum_{i = 1}^{K}n_i$ be the total number of test units and $M_i = \Sigma_{j = 1}^{n_i} r_{i, j}$ be the total number of measurements under stress level $\xi_i$. Further, define $M = \sum_{i = 1}^K{M_i}$ as the total number of measurements in the whole ADT and $\mathcal{T} = \sum_{i = 1}^{K}\sum_{j = 1}^{n_i}{t_{i, j, r_{i, j}}}$ as the total test time.

3.   Point estimations for model parameters

In this study, we mainly focus on the interval estimation of the proposed model. In order to develop the interval estimation procedures, we need to obtain the point estimations of model parameters first. That is because the point estimates are the basis for constructing PQs of interval estimations. Suppose that $n_i$ units are tested under the stress level $\xi_i$, and the degradation data $\mathbb{D}_i$ is provided. Based on $\mathbb{D}_i$, the log-likelihood function is given by

Let $X_i = \Sigma_{j = 1}^{n_i}\Sigma_{k = 1}^{r_{i, j}}\Delta X_{i, j, k}$, $\mathcal{T}_i = \Sigma_{j = 1}^{n_i}\Sigma_{k = 1}^{r_{i, j}}\Delta t_{i, j, k}$, and then the MLEs of $\mu_i$ and $\sigma_{i}^2$ are given as

Notice that

and we have the following factorization:

According to Cochran's theorem in ^[27], we can easily prove the following facts:

(ⅰ) $\widehat{\mu}_i\sim N(a+b\xi_i, \sigma_{i}^2/\mathcal{T}_i)$;

(ⅱ) $(M_i-1)S_{i}^2/\sigma_{i}^2\sim\chi^2(M_i-1)$;

(ⅲ) $\widehat{\mu}_i$ and $S_{i}^2$ are mutually independent.

3.1.   Point estimations for $a$ and $b$

Let $Y_i = X_i/ \mathcal{T}_i, i = 1, 2, \ldots, K$, and then $Y_i\sim N(\mu_i, \sigma_{i}^2/{\mathcal{T}_i})$. So, the mean and variance of $Y_i$ are given by $E(Y_i) = a+b\xi_i$ and $Var(Y_i) = \sigma_{i}^2/{\mathcal{T}_i}$, respectively. In order to get the estimates of parameters $a$ and $b$, the following linear regression model is considered.

Based on the linear regression model (3.2), the estimates of parameters $a$ and $b$ are provided in the following Theorem 3.1.

Theorem 3.1. Under the linear regression model (3.2), given the degradation data $\mathbb{D}$,

(1) The estimates of parameters $a$ and $b$ are given as

where

(2) The estimates $\widehat{a}$ and $\widehat{b}$ follow the normal distributions, that is,

(3) The covariance of the estimates $\widehat{a}$ and $\widehat{b}$ is given by

In addition, the degradation rate $\mu_0$ at normal use stress level $\xi_0$ can be estimated by $\widehat{\mu}_0 = \widehat{a}+\widehat{b} \xi_0$, and the estimate $\widehat{\mu}_0$ is also unbiased and has the variance $Var(\widehat{\mu}_0) = (G-2I\xi_0+F \xi_0^2)/(FG-I^2)$.

3.2.   Point estimations for $c$ and $d$

Notice that ${(M_i-1)S_i^2}/{\sigma_i^2}\sim \chi^2(M_i-1)$, and let $\Omega_i\hat{ = }\log[{(M_i-1)S_i^2}/{\sigma_i^2}]$. By calculating, the moment generating function (MGF) of $\Omega_i$ is derived.

Using the MGF $M_{\Omega_i}(t)$, the mean and variance of $\Omega_i$ are given by

respectively, where

Let $U_i\hat{ = }\log[{(M_i-1)S_i^2}/{2}]-\psi(\frac{M_i-1}{2})$, and by calculation we find that

To get the estimates of parameters $c$ and $d$, the following linear regression model is constructed.

Some properties about the estimates of $c$ and $d$ are given in the following Theorem 3.2.

Theorem 3.2. Under the linear regression model (3.5), given the degradation data $\mathbb{D}$.

(1) The estimates of parameters $c$ and $d$ are given as

where

(2) The estimates $\widetilde{c}$ and $\widetilde{d}$ are unbiased, that is, $E(\widetilde{c}) = c, E(\widetilde{d}) = d$.

(3) The variance and covariance of the estimates $\widetilde{c}$ and $\widetilde{d}$ are given by

Based on estimates $\widetilde{c}$ and $\widetilde{d}$, the diffusion parameter $\sigma_0^2$ at the normal use stress level $\xi_0$ can be estimated by $\widetilde{\sigma}_0^2 = \exp(\widetilde{c}+\widetilde{d}\xi_0)$. However, the estimate $\widetilde{\sigma}_0^2$ is biased. The following Theorem 3.3 gives an unbiased estimate of $\sigma_0^2$.

Theorem 3.3. Let $\widetilde{c}$ and $\widetilde{d}$ be the estimates of $c$ and $d$ defined in (3.6), and

Then,

(1) If $(M_i-1)/2+D_i > 0 (i = 1, 2, \ldots, K)$, an unbiased estimate of $\sigma_0^2$ can be given by

(2) If $(M_i-1)/2+2D_i > 0 (i = 1, 2, \ldots, K)$, the variance of $\widetilde{\sigma}_{0u}^2$ is given by

(3) If $(M_i-1)/2+2D_i > 0 (i = 1, 2, \ldots, K)$, the estimate $\widetilde{\sigma}_{0u}^2$ has a smaller mean square error than that of $\widetilde{\sigma}_0^2$.

4.   Interval estimations of model parameters and quantities

4.1.   GCIs for parameters $a$ and $b$

In this subsection, we try to derive the GCIs of parameters $a$ and $b$. In order to get the GPQs of $a, b$, we first develop the GPQs of parameters $\sigma_i^2 (i = 1, 2, \ldots, K)$. As $(M_i-1)S_{i}^2/\sigma_{i}^2\sim\chi^2(M_i-1)$, generating a copy $Q_{0, i}$ from the $\chi^2(M_i-1)$ distribution, the GPQ of $\sigma_i^2$ can be obtained by

It is worth emphasizing that $Q_{0, i}$ is treated as a known quantity in generalized inference ^[28].

Based on the model (3.2), substituting $S_i^2$ for the unknown parameter $\sigma_i^2$, the following weighted sum of squares is considered.

By minimizing (4.2), the estimates of $a, b$ are given as

Let $V_1 = \widetilde{a}-a, V_2 = \widetilde{b}-b$, and then $V_1, V_2$ can be presented as

where

It is obvious that in (4.3) and (4.4) the distributions of $V_1$ and $V_2$ only depend on the unknown parameter $\sigma_i^2$. So, generating a series of copies $Z_i^*$ from the standard normal distribution $N(0, 1)$, replace the unknown parameter $\sigma_i$ by its GPQ $\sqrt{\mathcal{W}_i}$ in (4.3) and (4.4) to get the quantities $V_1^*$ and $V_2^*$.

According to the substitute method given in ^[28,29], the GPQs of $a$ and $b$ are obtained by

Based on $\mathcal{G}_1$ and $\mathcal{G}_2$, the confidence intervals of $a$ and $b$ can be constructed. Because the distributions of $\mathcal{G}_1$ and $\mathcal{G}_2$ are very complicated, a simulation procedure can be used. Let $\mathcal{G}_{i, \alpha}$ be the $\alpha$ percentile of $\mathcal{G}_i$, and then $[\mathcal{G}_{1, \alpha/2}, \mathcal{G}_{1, 1-\alpha/2}]$ and $[\mathcal{G}_{2, \alpha/2}, \mathcal{G}_{2, 1-\alpha/2}]$ are the $1-\alpha$ level GCIs of $a$ and $b$, respectively.

4.2.   GCIs for parameters $c$ and $d$

In this subsection, we will derive the GCIs of model parameters $c$ and $d$. Based on (3.6) in Theorem 3.2, let $V_3 = \widetilde{c}-c$, $V_4 = \widetilde{d}-d$, and then $V_3$ and $V_4$ can be represented as

where $Q_i\hat{ = }(M_i-1)S_i^2/\sigma_i^2 \sim \chi^2(M_i-1)$. From (4.8) and (4.9) we know that the distribution of $V_3$ and $V_4$ only depend on the unknown parameters $\sigma_i^2$ through $Q_i$. So, we can generating a series of $Q^{*}_i$ from the distribution $\chi^2(M_i-1)$ to replace $Q_i$ and get $V^*_3$ and $V^*_4$.

Hence, the GPQs of parameters $c$ and $d$ can be given by

Let $\mathcal{G}_{i, \alpha}$ denote the $\alpha$ percentile of $\mathcal{G}_i$, and then $[\mathcal{G}_{3, \alpha/2}, \mathcal{G}_{3, 1-\alpha/2}]$ and $[\mathcal{G}_{4, \alpha/2}, \mathcal{G}_{4, 1-\alpha/2}]$ are the $1-\alpha$ level GCIs of $c$ and $d$, respectively.

4.3.   GCIs for quantities $\mu_0$, $\sigma_0^2$, $R(t_0)$ and $E(T)$

In practical applications, some important quantities for the Wiener double ADT model at the normal use stress level $\xi_0$, such as $\mu_0$, $\sigma_0^2$ and the reliability function of lifetime $T$, may be of more interest than the model parameters $(a, b, c, d)$. However, because these quantities involve two or more parameters, their interval estimations tend to be difficult. Similar to the cases of parameters $(a, b, c, d)$, we can develop the GCIs for these quantities.

Notice that $\mu_0, \sigma_0^2$ and the reliability function of lifetime $T$ are given by $\mu_0 = a+b\xi_0, \sigma_0^2 = \exp(c+d\xi_0)$ and $R(t_0) = 1-F_{T}\left(t_0| \mu_0, \sigma_0^2\right)$, respectively. According to the substitution method given in ^[28], the GPQs for $\mu_0, \sigma_0^2$ and $R(t_0)$ are given by

respectively.

Let $\mathcal{G}_{i, \alpha}$ denote the $\alpha$ percentile of $\mathcal{G}_i$. Then, $[\mathcal{G}_{i, \alpha/2}, \mathcal{G}_{i, 1-\alpha/2}], i = 5, 6, 7$, are the $1-\alpha$ level GCIs of $\mu_0, \sigma_0^2$, and $R(t_0)$, respectively. The percentiles of $\mathcal{G}_i, i = 1, 2, \cdots, 7$ can be obtained by the following simulation, Algorithm 1.

Algorithm 1: GCIs for $a, b, c, d$ and quantities $\mu_0, \sigma_0^2$ and $R(t_0)$.

(1) Given data set $\{(\Delta X_{i, j, k}, \Delta t_{i, j, k}, \xi_i), i = 1, \cdots, K; j = 1, \cdots, n_i; k = 1, \cdots, r_{i, j}\}$, compute $X_{i}, S_i^2$ and $U_i$.

(2) Generate a series of $\{Q_{0, i}\}_{i = 1}^K$ from $\chi^2(M_i-1)$, and then compute a series of $\{\mathcal{W}_{i}\}_{i = 1}^K$ through Eq (4.1).

(3) Generate a series of $\{Z_{i}^*\}_{i = 1}^K$ from $N(0, 1)$, based on $\{\mathcal{W}_{i}\}_{i = 1}^K$, and through Eq (4.7) compute $\mathcal{G}_1$ and $\mathcal{G}_2$.

(4) Generate a series of $\{Q^*_i\}_{i = 1}^K$ from $\chi^2(M_i-1)$, and through Eq (4.12) compute $\mathcal{G}_3$ and $\mathcal{G}_4$.

(5) Based on $\mathcal{G}_1, \mathcal{G}_2, \mathcal{G}_3$ and $\mathcal{G}_4$, through Eqs (4.13)–(4.15) compute $\mathcal{G}_5, \mathcal{G}_6$ and $\mathcal{G}_7$.

(6) Repeat (2)–(5) $B$ times, and then $B$ values of $\mathcal{G}_i, i = 1, 2, \cdots, 7$ are obtained, respectively.

(7) Arrange all $\mathcal{G}_i$ values in ascending order: $\mathcal{G}_{i, (1)} < \mathcal{G}_{i, (2)} < \cdots < \mathcal{G}_{i, (B)}, i = 1, 2, \cdots, 7$. Then, the $\alpha$ percentile of $\mathcal{G}_i$ is estimated by $\mathcal{G}_{i, (\alpha B)}$.

Remark 1. As is known to all, at the normal use stress level $\xi_0$, the unit's mean lifetime $E(T) = L/\mu_0$ is the monotonic function of $\mu_0$. Therefore, the GCI of $E(T)$ can be derived from the GCI of $\mu_0$.

4.4.   GPI for $X_0(t)$ and RUL$(\tau)$

In practical applications, the prediction intervals of the degradation characteristic, the lifetime and the remaining useful lifetime (RUL) of unit under normal use condition may be more practical and interesting for the product designer and user. So, it is important and meaningful to discuss the prediction interval of degradation characteristic $X_0(t)$ and RUL$(\tau)$. Unfortunately, for the proposed Wiener double ADT model, it is hard to obtain the exact prediction interval of these quantities, so we develop the GPI for them.

At the normal use stress level $\xi_0$, the degradation characteristic $X_0(t)$ can be presented as

Notice that given the failure threshold $L$, the product's lifetime $T$ under the normal operating condition follows an inverse Gaussian distribution IG$(L/\mu_0, L^2/\sigma_0^2)$. Given a fixed time $\tau$, if the degradation characteristic $X_0(\tau) = x_\tau$ is known, the remaining useful life RUL$(\tau)$ is defined as

Note that $X_0(t)$ has stationary independent increments, so $RUL(\tau)\sim IG\left((L-x_{\tau})/\mu_0, (L-x_{\tau})^2/\sigma_0^2\right)$. Based on GPQs $\mathcal{G}_5$ and $\mathcal{G}_6$, using the substitution method given in ^[33], the GPQs of $X_0(t)$ and RUL$(\tau)$ are obtained by

Let $\mathcal{G}_{i, \alpha}$ be the $\alpha$ percentile of $\mathcal{G}_i$, and then $[\mathcal{G}_{8, \alpha/2}, \mathcal{G}_{8, 1-\alpha/2}]$ and $[\mathcal{G}_{9, \alpha/2}, \mathcal{G}_{9, 1-\alpha/2}]$ are the $1-\alpha$ level GPIs of degradation characteristic $X_0(t)$ and RUL. The percentiles of $\mathcal{G}_9$ and $\mathcal{G}_{10}$ can be obtained by the following simulation, Algorithm 2.

Algorithm 2: GPIs for $X_0(t)$ and RUL$(\tau)$.

(1) Given data set $\{(\Delta X_{i, j, k}, \Delta t_{i, j, k}, \xi_i), i = 1, \cdots, K; j = 1, \cdots, n_i; k = 1, \cdots, r_{i, j}\}$, compute $X_{i}, S_i^2$ and $U_i$.

(2) Generate a series of $\{Q_{0, i}\}_{i = 1}^K$ from $\chi^2(M_i-1)$, and then compute a series of $\{\mathcal{W}_{i}\}_{i = 1}^K$ through Eq (4.1).

(3) Generate a series of $\{Z_{i}^*\}_{i = 1}^K$ from $N(0, 1)$, based on $\{\mathcal{W}_{i}\}_{i = 1}^K$, and through Eq (4.7) compute the $\mathcal{G}_1$ and $\mathcal{G}_2$.

(4) Generate a series of $\{Q^*_i\}_{i = 1}^K$ from $\chi^2(M_i-1)$, and through Eq (4.12) compute $\mathcal{G}_3$ and $\mathcal{G}_4$.

(5) Based on $\mathcal{G}_1, \mathcal{G}_2, \mathcal{G}_3$ and $\mathcal{G}_4$, through Eqs (4.13) and (4.14) compute $\mathcal{G}_5$ and $\mathcal{G}_6$.

(6) Based on $\mathcal{G}_5, \mathcal{G}_6$, through Eqs (4.16) and (4.17) compute $\mathcal{G}_8$ and $\mathcal{G}_{9}$.

(7) Repeat (2)–(6) $B$ times, and $B$ values of $\mathcal{G}_8, \mathcal{G}_9$ are obtained, respectively.

(8) Arrange all $\mathcal{G}_i$ values in ascending order: $\mathcal{G}_{i, (1)} < \mathcal{G}_{i, (2)} < \cdots < \mathcal{G}_{i, (B)}, i = 8, 9$. Then, the $\alpha$ percentile of $\mathcal{G}_i$ is estimated by $\mathcal{G}_{i, (\alpha B)}$.

Remark 2. Note that at the normal use stress level $\xi_0$, the unit's lifetime $T$ follows the $IG(L/\mu_0, L^2/\sigma_0^2)$ distribution. So, the GCI of lifetime $T$ can be obtained by the GPQ $\mathcal{G}_{10}$, where $\mathcal{G}_{10}\sim IG(L/\mathcal{G}_5, L^2/\mathcal{G}_6)$.

5.   Simulation study and comparative analysis

In this section, a Monte Carlo simulation study is implemented to evaluate the proposed GCIs/GPIs of model parameters and some quantities in terms of coverage percentage (CP) and average interval length (AL). Without loss of generality, we consider an ADT with three stress levels $\xi_1 = 1, \xi_2 = 2$ and $\xi_3 = 3$, and the normal use stress level is $\xi_0 = 0.5$. Four parameter settings are selected and given in Table 1. For convenience, the values of $n_i, r_{i, j}, \Delta t_{i, j, k}$ are chosen to be $n_1 = \cdots = n_K\hat{ = }n = 5, 8, 10; r_{i, j}\hat{ = }r = 4, 6, 8$; and $\Delta t_{i, j, k} = 1$. Three combinations of $(n, r, \Delta t_{i, j, k})$ are considered in the simulation: $(5, 4, 1), (8, 6, 1), (10, 8, 1)$. A total of 12 combinations of $(a, b, c, d, L)$ and $(n, r, \Delta t_{i, j, k})$ are examined in the simulation study. We take $B = 5000$ in the simulation study, and all the results are based on 5000 replications. The simulation results are provided in Tables 2–4.

Tables 2–4 summarize the CPs and ALs of the two-sided equal-tailed 90% and 95% GCIs/GPIs for model parameters and some quantities under the cases of Ⅱ, Ⅲ and Ⅳ. It is observed from the simulation results that the CPs of the proposed GCIs/GPIs are quite close to the nominal levels, even for small sample sizes. Based on the normal approximation to the binomial distribution, CPs between 94% and 96% are considered appropriate for the 95% confidence intervals. For fixed parameter settings, when the sample size $n$ and the number of measurements $r$ increase, the ALs of GCIs/GPIs decrease as expected. These findings show that the proposed confidence interval procedures work well, and the performances of the proposed GCIs/GPIs are satisfactory with respect to the CPs.

As is known to all, the parametric bootstrap method is a classic approach to obtain confidence intervals for model parameters. In order to fully evaluate the performances of the GCIs/GPIs, we also consider the bootstrap CIs for the Wiener double ADT model. For comparison, the confidence limits (CLs), such as lower confidence limit (LCL) and upper confidence limit (UCL), for model parameters and some quantities are also examined. We performed a comparative analysis of the CIs, LCLs and UCLs obtained by the GPQ method and the bootstrap-$p$ method. For saving space, we only give the simulation results under the case Ⅰ, and they are provided in Tables 5–7. The bootstrap-$p$ procedure is also based on 5,000 bootstrap samples.

Tables 5–7 show that the CPs of the CIs, LCLs and UCLs obtained by the GPQ method are all very close to the nominal levels, even for small sample sizes. However, the CPs of the CIs, LCLs and UCLs obtained by the bootstrap-$p$ method are not close to the nominal levels for some parameters and quantities. In particular, from Table 5 we find that the bootstrap-$p$ CIs of RUL$(4)$ are not close to the nominal levels. In addition, from Tables 6 and 7 we also find that the LCLs and UCLs obtained in the bootstrap-$p$ method perform badly. For example, the CPs of model parameters $c, \sigma_0^2$, reliability function $R(5)$, mean lifetime $E(T)$, degradation amount $X_0(10)$ and RUL$(4)$ deviate from the nominal levels.

As the sample size $n$ increases, the CPs of the bootstrap-$p$ CIs/PIs approach the nominal levels. Tables 5–7 also indicate that, for fixed parameter settings, as the sample size $n$ increases, the ALs become shorter, the LCLs become larger, and the UCLs become smaller for both GPQ method and bootstrap-$p$ method as expected. These findings indicate that the CIs, LCLs and UCLs obtained in the GPQ method perform better than the corresponding bootstrap-$p$ ones in terms of CP. Therefore, we recommend the proposed CIs, LCLs and UCLs in the GPQ method for the proposed Wiener double ADT model, especially in small sample cases.

6.   An illustrative example

In this section, we provide the ADT data of commercial white LEDs to illustrate our proposed methods. Degradation in lumen maintenance is the main failure mechanism for LEDs ^[30]. An LED is defined as a failure when the lumen maintenance decreases by 30% of its initial level.

In the ADT, 16 retrofit LED tubes based on a low-power LED are assigned to be tested at stress levels $s_1 = 25\; ^\circ$C and $s_2 = 55\; ^\circ$C, respectively. In general, the normal operating temperature of the LED is $s_0 = 25\; ^\circ$C. The tubes are placed in the climate chamber to ensure the stability of the ambient operating temperature, and the lumen outputs of the tubes are measured regularly. Note that there exists a very short period of increase in the lumen output due to incomplete burn-in ^[31], so this period is identified from the original data and is discarded in the analysis. For each test unit, the decrease in lumen output is normalized by the initial output. In addition, there is an unexpected catastrophic failure, and only 15 tubes under $55\; ^\circ$C are recorded. Figure 1 shows the degradation paths of the LEDs. Similar to Hong et al. ^[23], to protect proprietary information, a power-law time-scale transformation with exponent 0.70 is used to linearize the degradation paths.

The degradation data are not monotone, so the Gamma process and inverse Gaussian process are not suitable to deal with them. Thus, we choose the Wiener process to model the degradation data. In order to assess the goodness-of-fit, a standard normal Q-Q plot for the average standardized degradation increments is given in Figure 2. The points scatter around the line nicely except one point at the left tail, probably due to measurement errors for some observations.

Since temperature acts as the accelerated stress in this example, the Arrhenius model is used to describe the relationship between parameters and the accelerated stress. Hence, we take

Similar to Lim and Yum ^[32], the standardized stress levels are given by

Since $s_0 = s_1$, the temperature levels are standardized as $\zeta_0 = \zeta_1 = 0, \zeta_2 = 1$.

Let $\sigma_i^2$ be replaced by its estimate $S_i^2$ in (3.3), and the point estimates of the parameters $a$ and $b$ are given by $\tilde{a} = 0.2026, \tilde{b} = 0.0619$. Based on Eq (3.6), the point estimates of the parameters $c$ and $d$ are obtained by $\tilde{c} = 3.9533, \tilde{d} = 0.2960$. Notice that $\mu_i = a+b\zeta_i, \sigma_i^2 = \exp(c+d)\zeta_i, i = 1, 2$. The estimates of $\mu_i$ and $\sigma_i^2$ are given by $\tilde{\mu}_1 = 0.2044, \tilde{\mu}_2 = 0.2645, \tilde{\sigma}_1^2 = 52.5497, \tilde{\sigma}_2^2 = 70.0571$, respectively.

Notice that the mean degradation path at the stress level $s_i$ is $E(X_i(t)) = \mu_i t$, and the estimate of $E(X_i(t))$ is then given by $\widetilde{E(X_i(t))} = \tilde{\mu}_i t$ (considering the transformed time scale, here the time $t$ should be replaced by $t^{0.7}$). The degradation paths and the estimates of these mean degradation paths under transformed time scale are given in Figure 1. It is obvious that these degradation data are fitted well by the proposed Wiener doubly accelerated degradation model.

To illustrate the GPQ method described in section 4, we use the degradation data in Figure 1 to derive the GCIs/GPIs of model parameters and some reliability metrics. Based on the $\mathcal{G}_i, i = 1, 2, \cdots, 6$, the GCIs of model parameters $a, b, c, d$ and quantities $\mu_0, \sigma_0^2$ are obtained, for nominal levels 90% and 95%. According to $\mathcal{G}_7$, the GCIs of reliability function $R(100)$ are also obtained. In addition, using $\mathcal{G}_8$ and $\mathcal{G}_9$, the GPIs of degradation amount $X_0(150)$ and remaining useful life RUL(126) are developed. All the results are based on 10,000 replications and presented in Table 8.

7.   Conclusions

In this paper we proposed a doubly accelerated degradation test model of Wiener process, in which both the degradation rate and the diffusion parameter are affected by the stress level. The point estimates of model parameters are obtained by constructing a regression model. Furthermore, based on the point estimates of model parameters, the GCIs of model parameters are developed by constructing GPQs. Utilizing the substitute method, the GCIs of some quantities and reliability metrics are derived, and the GPIs of degradation amount $X_0(t)$ and RUL for units under normal use stress level are also developed.

Extensive simulations are carried out to evaluate the performances of the proposed procedures. The simulation results reveal that the proposed confidence intervals performed well in terms of CPs, and compared with the traditional bootstrap method. Finally, an illustrative example is provided to demonstrate the proposed procedures. As is known to all, the Gamma process and inverse Gaussian process are also widely used in analysis of accelerated degradation data. In the future, we shall mainly focus on studying the small sample inferential methods for doubly accelerated degradation test models based on the Gamma and inverse Gaussian process.

Acknowledgments

The authors thank the Editor and the reviewers for their detailed comments and suggestions, which considerably helped improve the manuscript. The work is supported by the Pre-research Project of the National Science Foundation of Anhui Polytechnic University (Xjky08201903), the Talent Cultivation and Research Start-up Foundation of Anhui Polytechnic University (S022022014), the Support program for outstanding young talents in colleges and universities of Anhui Province (gxyqZD2022046) and the National Social Science Foundation of China (18BTJ034, 20BTJ048).

Conflict of interest

The authors declare that they have no conflicts of interest.

Supplementary

A.   Proof of Theorem 3.1

Let $V = diag(\sigma_1^2/T_1, \cdots, \sigma_K^2/T_K), Y = (Y_1, \cdots, Y_K)^T$, and

Then, the estimates $(\widehat{a}, \widehat{b})$ are given by

So, we have the expectation of $\widehat{a}$ and $\widehat{b}$

Furthermore, the covariance matrix of the estimates $(\widehat{a}, \widehat{b})$ is given by

So, the variance and covariance of the estimates $\widehat{a}$ and $\widehat{b}$ are obtained by

Notice that the $X_i$'s are normal distributions and independent of each other and $\widehat{a}$ and $\widehat{b}$ are linear combinations of them, so the estimates $\widehat{a}$ and $\widehat{b}$ are also normal distributions. That is,

The proof of Theorem 3.1 is completed.

B.   Proof of Theorem 3.2

Let $V = diag(\psi^{'}(\frac{M_1-1}{2}), \cdots, \psi^{'}(\frac{M_K-1}{2})), Y = (U_1, \cdots, U_K)^T$. Similar to the proof of Theorem 3.1, Theorem 3.2 can be easily proved, and here we neglect the detailed proof.

C.   Proof of Theorem 3.3

Let $\zeta_i = \log\left(\frac{(M_i-1)S_i^2}{2\sigma_i^2}\right)-\psi(\frac{M_i-1}{2})$, and then $E(\zeta_i) = 0, Var(\zeta_i) = \psi^{'}(\frac{M_i-1}{2})$. Notice that

So, we have

So,

Because $\frac{(M_i-1)S_i^2}{\sigma_i^2}\sim \chi^2(M_i-1)$, and

we have

From the above formulas, the unbiased estimate of $\sigma_0^2$ and its variance $Var(\widetilde{\sigma}_{0u}^2)$ can be given as

Similar to the proof of Theorem $6$ in Wang and Yu ^[33], we can prove that $\widetilde{\sigma}_{0u}^2$ has a smaller mean square error than that of $\widetilde{\sigma}_0^2$, where we neglect the detailed proof. The proof of Theorem 3.3 is completed.