# Misspecification analysis of gamma- and inverse Gaussian-based perturbed degradation processes

Nicola Esposito, Agostino Mele, Bruno BC Castanier, Massimiliano Giorgio

► To cite this version:

Nicola Esposito, Agostino Mele, Bruno BC Castanier, Massimiliano Giorgio. Misspecification analysis of gamma- and inverse Gaussian-based perturbed degradation processes. Applied Stochastic Models in Business and Industry, 2023, 237 (3), pp.109320. 10.1002/asmb.2824. hal-04631570

HAL Id: hal-04631570

https://hal.science/hal-04631570v1

Submitted on 2 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Misspecification analysis of gamma- and inverse Gaussian-based perturbed degradation processes

Nicola Esposito<fnref n="1" /> ID | Agostino Mele<fnref n="2" />,<fnref n="3" /> ID | Bruno Castanier<fnref n="4" /> ID | Massimiliano Giorgio<fnref n="1" /> ID

## Correspondence

Massimiliano Giorgio, Università di Napoli Federico II, Napoli, Italy.
Email: massimiliano.giorgio@unina.it

## Funding information

Université Franco Italien within the frame of the chapitre 2 of Vinci project, Grant/Award Number: subvention N° C2-221; Università di Napoli Federico II in the frame of the international agreement between Dipartimento di Ingegneria Industriale and Polytech Angers, Grant/Award Numbers: codice identificativo 000011--ALTRI-2021-M-GIORGIO_001_001, prot. 40162 del 20/4/2021; WISE project of the Région Pays de la Loire (France)

## Abstract

Albeit not equivalent, in many applications the gamma and the inverse Gaussian processes are treated as if they were. This circumstance makes the misspecification problem of these models interesting and important, especially when data are affected by measurement errors, since noisy/perturbed data do not allow to verify whether the selected model is actually able to adequately fit the real (hidden) degradation process. Motivated by the above considerations, in this paper we conduct a large Monte Carlo study to evaluate whether and how the presence of measurement errors affects this misspecification issue. The study is performed considering as reference models a perturbed gamma process recently proposed in the literature and a new perturbed inverse Gaussian process that share the same non-Gaussian distributed error term. As an alternative option, we also analyze the more classical case where the error term is Gaussian distributed. We consider both the situation where the true model is the perturbed gamma and the one where it is the perturbed inverse Gaussian. Model parameters are estimated from perturbed data using the maximum likelihood method. Estimates are retrieved by using a new sequential Monte Carlo EM algorithm, which use allows to hugely mitigate the severe numerical issues posed by the direct maximization of the likelihood. The risk of incurring in a misspecification is evaluated as percentage of times the Akaike information criterion leads to select the wrong model. The severity of a misspecification is evaluated in terms of its impact on maximum likelihood estimate of the mean remaining useful life.

## KEYWORDS

expectation-maximization algorithm, gamma process, inverse Gaussian process, measurement errors, model misspecification, particle filter

# 1 | INTRODUCTION

The gamma and the inverse Gaussian processes are widely applied in engineering and reliability.<fnref n="1" />⁻¹² Indeed, both are considered natural choices for modeling monotonic degradation phenomena. In particular, even though they are not fully equivalent, in many applications these models are treated as if they were. This situation makes the model misspecification

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2023 The Authors. Applied Stochastic Models in Business and Industry published by John Wiley & Sons Ltd.

<footnote n="1">¹Università di Napoli Federico II, Napoli, Italy</footnote>

<footnote n="2">²Dipartimento di Ingegneria, Università degli studi della Campania “Luigi Vanvitelli”, Aversa, Italy</footnote>

<footnote n="3">³Kineton, Napoli, Italy</footnote>

<footnote n="4">⁴Laris, Polytech Angers/Université d'Angers, Angers, France</footnote>
issue of gamma and inverse Gaussian processes interesting and important. However, discriminating between these models is not a simple task, especially when available experimental data are contaminated by measurement errors, a situation that is often encountered in practical applications, mainly when data are collected through in-service and/or non-destructive inspection methods.¹³ In fact, noisy/perturbed data do not allow to directly evaluate the ability of the model to describe the hidden degradation process, which is the primary aim of using a degradation model selection criterion in reliability and maintenance applications. Despite its importance, the misspecification issue of gamma and inverse Gaussian processes has been rarely addressed in the literature. In fact, more specifically, while the basic problem is addressed in Zhang and Revie¹⁴ and Tseng and Yao,¹⁵ to the best of our knowledge, in the literature there are no papers dealing with this misspecification issue that focus on situations where available data are contaminated by measurement errors.

Motivated by these considerations, in this article, we conduct a large Monte Carlo study aimed at investigating whether and how the presence of measurement errors affects the misspecification issue of gamma and inverse Gaussian processes.

The study is carried out considering as competing models a perturbed gamma process (PGP) recently proposed by Giorgio et al.¹⁶ and a new perturbed inverse Gaussian process (PIGP) that share the same error term. In addition, to facilitate the comparative analysis, the (hidden) inverse Gaussian process (IGP) is formulated by using a new parameterization that allows the considered competing models to share the same parameters and the same functional forms of the mean and variance functions.

As in Giorgio et al.,¹⁶ the error term is supposed to depend (in stochastic sense) on the hidden degradation level and, conditionally to the degradation level, is modeled as a 3-parameter inverse gamma random variable. This modeling solution distinguishes the considered models from other perturbed degradation models suggested in the literature, where the error is modeled by using a Gaussian distribution and is assumed to be stochastically independent of the hidden degradation process.¹³,¹⁷⁻²⁰ Moreover, it ensures that the perturbed measurement is nonnegative, a result that is not guaranteed in the case where the error term is described by using the Gaussian model, especially when the magnitude of the standard deviation of the error term is comparable to that of the measured degradation level.

Nonetheless, for the sake of generality, as an alternative modeling solution, we have also examined the case where the error term is modeled by using a Gaussian distribution. Also in this second case, by following Pulcini,²¹ we have assumed that the measurement error depends in stochastic sense on the measured degradation level. Yet, as a special case of this latter model, we have also considered the classical assumption where the error is independent of the hidden degradation process.

Both the misspecification of a PGP with a PIGP and the symmetric case of the misspecification of a PIGP with a PGP are considered. Model parameters are estimated from perturbed data by using the maximum likelihood (ML) method.

The fitting ability of the considered competing perturbed models is evaluated by using the Akaike information criterion (AIC).²² The risk of incurring in a misspecification is evaluated as percentage of times the AIC leads to select the wrong model. The severity of a misspecification is evaluated in terms of its impact on ML estimate of the mean remaining useful life (MRUL). The impact of the presence of measurement errors on risk and consequences of incurring in a misspecification is evaluated by comparing the obtained results with those obtained by carrying out the same misspecification analysis in the absence of measurement errors.

Unfortunately, as it also occurs in the case of gamma- and inverse Gaussian-based perturbed degradation models, computing the likelihood functions, which are not available in closed form, requires intensive numerical methods that, at the same time, increase the computational burden and exacerbate convergence issues of numerical optimization algorithms used to retrieve the ML estimates. Indeed, mainly due to numerical problems and long computational times, the performance of ML estimators of parameters of the gamma- and inverse Gaussian-based perturbed degradation processes, and/or functions thereof, are typically investigated by using a relatively small number (e.g., 100 or 200) of synthetic datasets (i.e., number of replications), a situation that does not allow to obtain results characterized by an adequate degree of accuracy.

In order to overcome this limitation, as an additional contribution of this paper, we suggest and adopt a new sequential Monte Carlo expectation-maximization (EM) algorithm, which allows to drastically simplify the estimation task and (consequently) to find a better tradeoff between precision and computational affordability of the Monte Carlo study. In fact, the use of this algorithm allows us to perform a Monte Carlo study where each index of interest is evaluated on the basis of 2000 simulated datasets.

However, even by adopting the suggested algorithm, the considered misspecification study remains a very time consuming task. This is mainly because a correct application of the Akaike information criterion requires a very accurate
evaluation of ML estimates of the 5-parameter log-likelihood functions of the competing models, which in about 6% of cases differ only in the fourth significant figure and more rarely (in less than 1% of cases) only in the fifth significant figure. In the case of the suggested sequential Monte Carlo EM algorithm, the needed accuracy can be obtained only by using a very high number of particles (say, 500,000 or more), a circumstance that significantly affects the computational burden.

The rest of the paper is structured as it follows. Section 2 introduces the considered competing perturbed degradation models. Section 3 deals with the formulation of the cumulative distribution function of the remaining useful life. Section 4 is devoted to the formulation of the likelihood function. Section 5 describes the performed misspecification study in detail, illustrates the Monte Carlo approach adopted to develop it, and reports obtained results and related comments. Section 6 provides some conclusions. Appendices A and B illustrate the EM algorithm and the particle filter method used to compute the ML estimates of model parameters and functions thereof.

## 2 | PERTURBED GAMMA AND PERTURBED INVERSE GAUSSIAN DEGRADATION PROCESSES

A perturbed stochastic process {$Z(t); t \ge 0$} is customarily formulated as in (1):

$$ \{Z(t) = W(t) + \varepsilon(t); t \ge 0\}, \quad (1) $$

where {$W(t); t \ge 0$} is the actual (hidden) degradation process and $\varepsilon(t)$ is a perturbing term, here intended as a measurement error. To provide a full stochastic description of the process {$Z(t); t \ge 0$} it is necessary to define the hidden process {$W(t); t \ge 0$}, the error term $\varepsilon(t)$, and their mutual stochastic relationship.

In this paper, it is assumed that:

i. the hidden process might be either a gamma process (GP) or an IGP. Given that both processes have independent increments, both are fully defined by the distribution of their generic increment $\Delta W(t, t + \tau) = W(t + \tau) - W(t)$ and by an initial condition, here $W(0) = 0$.

The probability density function (pdf) of the increment $\Delta W(t, t + \tau)$ of the GP is expressed by using the classical parametrization (2):

$$ f_{\Delta W(t,t+\tau)}(\delta) = \frac{\delta^{\Delta\eta(t,t+\tau)-1}}{\theta^{\Delta\eta(t,t+\tau)} \cdot \Gamma(\Delta\eta(t,t+\tau))} \cdot e^{-\delta/\theta}, \quad \delta > 0, \quad (2) $$

where $\Gamma(\cdot)$ is the complete gamma function, $\theta(\theta > 0)$ is the scale parameter, $\eta(t)$ is a non-negative monotonic increasing function, here referred to as the age function, and $\Delta\eta(t, t + \tau) = \eta(t + \tau) - \eta(t)$.

Differently, in order to facilitate the misspecification study, the pdf of the increment $\Delta W(t, t + \tau)$ of the IGP is expressed by using the following special form:

$$ f_{\Delta W(t,t+\tau)}(\delta) = \frac{\Delta\eta(t,t+\tau)}{\sqrt{2 \cdot \pi \cdot \theta^{-1} \cdot \delta^3}} \cdot e^{-\frac{\delta-\theta \cdot \Delta\eta(t,t+\tau)^2}{2 \cdot \theta \cdot \delta}}, \quad \delta \ge 0, \quad (3) $$

which can be easily obtained from the classical functional form of the pdf of the increment of the IGP:⁸

$$ f_{\Delta W(t,t+\tau)}(\delta) = \sqrt{\frac{\lambda \cdot [\Delta\Lambda(t,t+\tau)]^2}{2 \cdot \pi \cdot \delta^3}} \cdot e^{-\frac{\lambda \cdot [\delta - \mu \cdot \Delta\Lambda(t,t+\tau)]^2}{2 \cdot \mu^2 \cdot \delta}}, \quad \delta \ge 0, \quad (4) $$

where $\Lambda(t)$ is a nonnegative monotonic increasing function, $\Delta\Lambda(t, t + \tau) = \Lambda(t + \tau) - \Lambda(t)$, and $\mu$ and $\lambda$ are positive valued parameters, by setting $\mu^2/\lambda = \theta$ and $\mu \cdot \Delta\Lambda(t, t + \tau) = \Delta\eta(t, t + \tau)$.

The main advantage of using the parametrization (3) is that it allows the considered competing processes to share the same parameters and the same functional forms of the mean $E\{W(t)\}$ and variance $V\{W(t)\}$ functions, which can be
expressed as in (5) and (6), respectively:

$$E\{W(t)\} = \theta \cdot \eta(t), \qquad (5)$$

$$V\{W(t)\} = \theta^2 \cdot \eta(t). \qquad (6)$$

In fact, under the classical parameterization (4) the mean and variance functions of the IGP are equal to $E\{W(t)\} = \mu \cdot \Lambda(t)$ and $V\{W(t)\} = [\mu \cdot \Lambda(t)]^3 / \lambda$, respectively.

From (2), the cumulative distribution function (cdf) of the increment $\Delta W(t, t + \tau)$ of the GP results in:

$$F_{\Delta W(t,t+\tau)}(\delta) = \frac{\gamma(\Delta\eta(t, t+\tau), \delta/\theta)}{\Gamma(\Delta\eta(t, t+\tau))}, \delta \ge 0, \qquad (7)$$

where $\gamma(\cdot)$ is the lower incomplete gamma function. Likewise, from (3), under the IGP the same cdf has the following (special) expression:

$$F_{\Delta W(t,t+\tau)}(\delta) = \Phi\left(\frac{\delta - \theta \cdot \Delta\eta(t, t+\tau)}{\sqrt{\theta \cdot \delta}}\right) + e^{\theta \cdot \Delta\eta(t,t+\tau)} \cdot \Phi\left(-\frac{\delta + \theta \cdot \Delta\eta(t, t+\tau)}{\sqrt{\theta \cdot \delta}}\right), \delta \ge 0, \qquad (8)$$

where $\Phi(\cdot)$ is the standard normal cdf.

In this article, both under the gamma and the inverse Gaussian models, the age function $\eta(t)$ is modeled by using the very flexible and largely adopted power law functional form $\eta(t) = (t/a)^b$.

ii. To model the perturbing term $\varepsilon(t)$ we consider two different options:

(1) by following Giorgio et al.,¹⁶ we assume that $\varepsilon(t)$ depends in a stochastic sense on the hidden degradation level $W(t)$ and that, given $W(t)$, is conditionally distributed as a 3-parameter inverse gamma random variable. The conditional pdf of the error term $\varepsilon(t)$, given $W(t) = w$, is expressed as follows:

$$f_{\varepsilon(t)|W(t)}(\varepsilon|w) = \frac{(\alpha(w))^{\beta(w)} \cdot (\varepsilon + w)^{-\beta(w)-1}}{\Gamma(\beta(w))} \cdot e^{-\frac{\alpha(w)}{\varepsilon+w}}, \varepsilon \ge -w, \qquad (9)$$

where $\beta(w) = \varphi \cdot w^{2-\nu} + 2$, $\varphi > 0$, $-\infty < \nu < \infty$, and $\alpha(w) = (\beta(w) - 1) \cdot w$. As remarked in Giorgio et al.,¹⁶ this modeling solution ensures that the perturbed measurement $Z(t) = W(t) + \varepsilon(t)$ is non-negative and enables to easily model cases where the measurement error depends (in stochastic sense) on the measured degradation level, two situations that are often encountered in practical applications.

Under this setting, the conditional mean and variance of $\varepsilon(t)$, given $W(t) = w$, can be expressed as in (10) and (11), respectively:

$$E\{\varepsilon(t)|W(t) = w\} = \frac{\alpha(w)}{\beta(w) - 1} - w = 0, \qquad (10)$$

$$V\{\varepsilon(t)|W(t) = w\} = \frac{(\alpha(w))^2}{(\beta(w) - 1)^2 \cdot (\beta(w) - 2)} = \frac{w^{\nu}}{\varphi}. \qquad (11)$$

Thus, from (1), the perturbed measurement $Z(t)$, given $W(t) = w$, has conditional pdf:

$$f_{Z(t)|W(t)}(z|w) = \frac{(\alpha(w))^{\beta(w)} \cdot z^{-\beta(w)-1}}{\Gamma(\beta(w))} \cdot e^{-\frac{\alpha(w)}{z}}, z \ge 0, \qquad (12)$$

conditional mean:

$$E\{Z(t)|W(t) = w\} = E\{\varepsilon(t)|W(t) = w\} + w = w, \qquad (13)$$
and conditional variance:

$$V\{Z(t)|W(t) = w\} = V\{\varepsilon(t)|W(t) = w\} = \frac{w^{\nu}}{\varphi}. \qquad (14)$$

(2) inspired by Pulcini,²¹ as an alternative modeling solution, we assume that $\varepsilon(t)$ is Gaussian distributed with zero mean and variance that depends on $W(t)$. Specifically, in this second case, the conditional pdf of the error term $\varepsilon(t)$, given $W(t) = w$, is formulated as:

$$f_{\varepsilon(t)|W(t)}(\varepsilon|w) = \frac{1}{\sqrt{2 \cdot \pi \cdot \sigma^2(w)}} \cdot e^{-\frac{1}{2} \cdot \frac{\varepsilon^2}{\sigma^2(w)}}, \quad -\infty < \varepsilon < \infty, \qquad (15)$$

where, $\sigma^2(w) = w^\nu/\varphi$, $\varphi > 0$, and $-\infty < \nu < \infty$.

Hence, from (1), the perturbed measurement $Z(t)$, given $W(t) = w$, has conditional pdf:

$$f_{Z(t)|W(t)}(z|w) = \frac{1}{\sqrt{2 \cdot \pi \cdot \sigma^2(w)}} \cdot e^{-\frac{1}{2} \cdot \frac{(z-w)^2}{\sigma^2(w)}}, \quad -\infty < z < \infty, \qquad (16)$$

conditional mean $E\{Z(t)|W(t) = w\} = w$, and conditional variance $V\{Z(t)|W(t) = w\} = w^\nu/\varphi$.

Indeed, this alternative setting has been specifically chosen because it allows $\varepsilon(t)$ and $Z(t)$, given $W(t) = w$, to have the same conditional mean and variance as they have under the error model (1). Pulcini²¹ used this modeling solution to formulate a PGP. Successively, Hao et al.²³ and Sun et al.²⁴ used it to formulate a PIGP.

iii. for any $n > 1$, the measurement error $\varepsilon(t_j)$ given $W(t_j)$ is conditionally independent both of $\varepsilon(t_k)$ and $W(t_k) \forall k \neq j$ ($j, k = 1, \dots, n$). Thus, equivalently, the perturbed observation $Z(t_j)$ given $W(t_j)$ is conditionally independent both of $Z(t_k)$ and $W(t_k) \forall k \neq j$, ($j, k = 1, \dots, n$).

Under these assumptions, from (5), (10), and (13), by using the law of total mean, the (marginal) means of $\varepsilon(t)$ and $Z(t)$, under all the considered perturbed processes (i.e., under all the possible combinations of hidden processes and error terms) result in:

$$E\{\varepsilon(t)\} = E\{E\{\varepsilon(t)|W(t)\}\} = E\{0\} = 0, \qquad (17)$$

and

$$E\{Z(t)\} = E\{E\{Z(t)|W(t)\}\} = E\{W(t)\} = \theta \cdot \eta(t). \qquad (18)$$

Similarly, from (6), (11), (13), and (14), by exploiting the law of total variance, the marginal variances of $\varepsilon(t)$ and $Z(t)$ are equal to:

$$V\{\varepsilon(t)\} = V\{E\{\varepsilon(t)|W(t)\}\} + E\{V\{\varepsilon(t)|W(t)\}\} = V\{0\} + E\{V\{\varepsilon(t)|W(t)\}\} = \frac{E\{(W(t))^\nu\}}{\varphi}, \qquad (19)$$

and

$$V\{Z(t)\} = V\{E\{Z(t)|W(t)\}\} + E\{V\{Z(t)|W(t)\}\} = V\{W(t)\} + E\left\{\frac{(W(t))^\nu}{\varphi}\right\} = \theta^2 \cdot \eta(t) + \frac{E\{(W(t))^\nu\}}{\varphi}, \qquad (20)$$

where the fractal moment $E\{(W(t))^\nu\}$ (in general) depends on the hidden process. In fact, in the case of the PGP $E\{(W(t))^\nu\}$ is equal to:

$$E\{(W(t))^\nu\} = \theta^\nu \cdot \frac{\Gamma(\eta(t) + \nu)}{\Gamma(\eta(t))},$$
while under the PIGP it is given by the following integral:

$$E\{(W(t))^\nu\} = \frac{\eta(t)}{\sqrt{2 \cdot \pi \cdot \theta^{-1}}} \cdot \int_0^\infty w^{\nu - \frac{3}{2}} \cdot e^{-\frac{[w - \theta \cdot \eta(t)]^2}{2 \cdot \theta \cdot w}} \cdot dw$$

that (in general) is not available in closed form.

However, if $\nu = 0, 1$, or $2$, also $V\{\varepsilon(t)\}$ and $V\{Z(t)\}$ do not depend on the hidden model. In fact, if $\nu = 0$ it is

$$V\{\varepsilon(t)\} = \frac{1}{\varphi}$$

and

$$V\{Z(t)\} = \theta^2 \cdot \eta(t) + \frac{1}{\varphi}.$$

Likewise, if $\nu = 1$ it is

$$V\{\varepsilon(t)\} = \frac{\theta \cdot \eta(t)}{\varphi}$$

and

$$V\{Z(t)\} = \theta^2 \cdot \eta(t) + \frac{\theta \cdot \eta(t)}{\varphi}.$$

Notably, in this case the ratio between $V\{Z(t)\}$ and $V\{W(t)\}$ does not depend on $t$:

$$\frac{V\{Z(t)\}}{V\{W(t)\}} = \frac{\theta^2 \cdot \eta(t) + \frac{\theta \cdot \eta(t)}{\varphi}}{\theta^2 \cdot \eta(t)} = 1 + \frac{1}{\theta \cdot \varphi} \qquad (21)$$

Again, by using similar arguments, if $\nu = 2$, given that under the considered parameterization the two hidden processes have the same mean and variance (and hence the same second moment), under both processes it is:

$$E\{(W(t))^2\} = \theta^2 \cdot \eta(t) \cdot (\eta(t) + 1)$$

and thus, it results:

$$V\{\varepsilon(t)\} = \theta^2 \cdot \frac{\eta(t) \cdot (\eta(t) + 1)}{\varphi}$$

and

$$V\{Z(t)\} = \theta^2 \cdot \left\{ \eta(t) + \frac{\eta(t) \cdot (\eta(t) + 1)}{\varphi} \right\}.$$

Finally, it is worth to remark that:

* both when the hidden process is the gamma or the inverse Gaussian one, the perturbed process $\{Z(t), t \ge 0\}$ is non-Markovian.

* under the error model (1), when $\nu = 0$, although the conditional mean and variance of the error term do not depend on $W(t)$, $\varepsilon(t)$ still stochastically depends on $W(t)$ because its support depends on it. In fact, from (9), under this setting it is $\varepsilon \ge -w$.

* under the error model (2), when $\nu = 0$, $\varepsilon(t)$ is stochastically independent of $W(t)$. Hence, in this latter case, the error model (2) reduces to the classical model adopted in the majority of existing literature on perturbed degradation models.¹³,¹⁷⁻²⁰
## 3 | CDF OF RUL

A unit is assumed to fail when its true degradation level passes a preassigned failure threshold $D$. Hence, the lifetime $X$ of the unit is defined as the first passage time of the hidden degradation process to the threshold $D$:

$$X = \inf\{x : W(x) > D\}.$$

Accordingly, its remaining useful life $RUL(t)$ at time $t$ is defined as follows:

$$RUL(t) = \max\{0, X - t\}.$$

It is assumed that failures are not self-announcing$^{25,26}$ and that a failed unit may continue to operate, albeit with (possibly) reduced performance. Thus, due to the presence of measurement errors, it is also supposed that a perturbed measurement of the degradation state of the unit does not permit to assess with certainty whether a unit is already failed or not.

Under these assumptions, by following Esposito et al.,$^{27}$ considered that both the gamma and inverse Gaussian processes are monotonic increasing and that both have independent increments, the conditional cdf $F_{RUL(t)|Z(t)}(\tau|\mathbf{z}(t))$ of the $RUL(t)$ given $Z(t) = \mathbf{z}(t)$, is formulated as follows:

$$ 
\begin{aligned}
F_{RUL(t)|Z(t)}(\tau|\mathbf{z}(t)) &= P[RUL(t) \le \tau | Z(t) = \mathbf{z}(t)] = P[W(t + \tau) > D | Z(t) = \mathbf{z}(t)] \\
&= 1 - F_{W(t+\tau)|Z(t)}(D|\mathbf{z}(t)) = 1 - \int_{0}^{D} F_{\Delta W(t,t+\tau)}(D - w) \cdot f_{W(t)|Z(t)}(w|\mathbf{z}(t)) \cdot dw,
\end{aligned} 
\quad (22) $$

where $Z(t) = \{Z(t_j); j \ge 1, t_j \le t\}$ denotes the set of measurements collected up to and included the time $t$ and $\mathbf{z}(t) = \{\mathbf{z}(t_j); j \ge 1, t_j \le t\}$ denotes its realization.

Note that, due to the presence of measurement errors, given that failures are not self-announcing, in general it results $F_{RUL(t)}(0|\mathbf{z}(t)) > 0$. Moreover, for $t < t_1$ (that is, before the first measurement time) $Z(t)$ should be intended as the empty set so that, in this case, $F_{W(t+\tau)|Z(t)}(D|\mathbf{z}(t))$ coincides with $F_{W(t+\tau)}(D)$.

From (22), given $Z(t) = \mathbf{z}(t)$, the conditional mean of the $RUL(t)$, $MRUL(t|Z(t) = \mathbf{z}(t))$, can be computed as:

$$ MRUL(t|Z(t) = \mathbf{z}(t)) = \int_{0}^{\infty} (1 - F_{RUL(t)}(\tau|\mathbf{z}(t))) \cdot d\tau = \int_{0}^{\infty} F_{W(t+\tau)|Z(t)}(D|\mathbf{z}(t)) \cdot d\tau. \quad (23) $$

Unfortunately, under both the considered models, the cdf $F_{W(t+\tau)|Z(t)}(D|\mathbf{z}(t))$ is not available in closed form. Hence, the cdf of $RUL(t)$ in (22) and the related conditional mean $MRUL(t|Z(t) = \mathbf{z}(t))$ in (23) are computed by using the particle filter algorithm illustrated in Appendix B.

## 4 | LIKELIHOOD FUNCTION

Let us consider $m$ identical units, operating under homogeneous conditions, that are subjected to periodic inspections, aimed at measuring their degradation level, and let us suppose that measurements are affected by errors. Then, let us denote by $t_{i,j}$ ($j = 1, \dots, n_i; n_i \ge 1$) the age of the $i$th unit ($i = 1, \dots, m$) at the epoch of the $j$th inspection, by $Z_{i,j} = Z(t_{i,j})$ the perturbed degradation level of the unit $i$ at $t_{i,j}$, and by $\mathbf{z}_{i,j}$ its realization.

With these notations, under both the considered models, the likelihood function $L(\xi; \mathbf{z})$ of the perturbed data can be formulated as:

$$ L(\xi; \mathbf{z}) = \prod_{i=1}^{m} \prod_{j=1}^{n_i} f_{Z_{i,j}|Z_{i,j-1}}(z_{i,j} | z_{i,j-1}), \quad (24) $$

where $\xi = (a, b, \theta, \phi, \nu)$ is the vector of model parameters, $\mathbf{z} = \{\mathbf{z}_{1,n_1}, \dots, \mathbf{z}_{m,n_m}\}$ is the realization of the whole set of available noisy measurements $\mathbf{Z} = \{\mathbf{Z}_{1,n_1}, \dots, \mathbf{Z}_{m,n_m}\}$, $\mathbf{Z}_{i,j} = \{Z_{i,1}, \dots, Z_{i,j}\}$ is the set of perturbed measurements of the
degradation state of the unit *i* collected up to time $t_{i,j}$, $\mathbf{z}_{i,j} = \{z_{i,1}, \dots, z_{i,j}\}$ is the realization of $\mathbf{Z}_{i,j}$, $t_{i,0} = 0$, $\mathbf{Z}_{i,0}$ and $\mathbf{z}_{i,0}$ are the empty set, and $f_{\mathbf{Z}_{i,1}|\mathbf{Z}_{i,0}}(z_{i,1}|\mathbf{Z}_{i,0}) = f_{\mathbf{Z}_{i,1}}(z_{i,1})$.

Unfortunately, due to the fact that {$Z(t); t \ge 0$} is not a Markov process, the pdfs in the likelihood function (24) are not available in closed form and their computation is extremely demanding. However, they can be efficiently computed, for any $i = 1, \dots, m$ and $j = 1, \dots, n_i$, by using the recursive Equations (25)–(27):

$$f_{W_{i,j}|\mathbf{Z}_{i,j-1}}(w_{i,j}|\mathbf{Z}_{i,j-1}) = \int_{0}^{w_{i,j}} f_{\Delta W_{i,j}|W_{i,j-1}}(w_{i,j} - w_{i,j-1}|w_{i,j-1}) \cdot f_{W_{i,j-1}|\mathbf{Z}_{i,j-1}}(w_{i,j-1}|\mathbf{Z}_{i,j-1}) \cdot dw_{i,j-1} \quad (25)$$

$$f_{\mathbf{Z}_{i,j}|\mathbf{Z}_{i,j-1}}(z_{i,j}|\mathbf{Z}_{i,j-1}) = \int_{0}^{\infty} f_{\mathbf{Z}_{i,j}|W_{i,j}}(z_{i,j}|w_{i,j}) \cdot f_{W_{i,j}|\mathbf{Z}_{i,j-1}}(w_{i,j}|\mathbf{Z}_{i,j-1}) \cdot dw_{i,j}, \quad (26)$$

$$f_{W_{i,j}|\mathbf{Z}_{i,j}}(w_{i,j}|\mathbf{Z}_{i,j}) = \frac{f_{\mathbf{Z}_{i,j}|W_{i,j}}(z_{i,j}|w_{i,j}) \cdot f_{W_{i,j}|\mathbf{Z}_{i,j-1}}(w_{i,j}|\mathbf{Z}_{i,j-1})}{f_{\mathbf{Z}_{i,j}|\mathbf{Z}_{i,j-1}}(z_{i,j}|\mathbf{Z}_{i,j-1})}, \quad (27)$$

where $W_{i,j} = W(t_{i,j})$ denotes the hidden (true) degradation level of the unit *i* at time $t_{i,j}$, $\Delta W_{i,j} = \Delta W(t_{i,j-1}, t_{i,j}) = W_{i,j} - W_{i,j-1}$ is the hidden (true) degradation increment of the unit *i* in the interval $(t_{i,j-1}, t_{i,j})$, $w_{i,j}$ is the realization of $W_{i,j}$, and $\Delta w_{i,j} = w_{i,j} - w_{i,j-1}$ is the realization of $\Delta W_{i,j}$. Details about the derivation of (25)–(27) are given in Esposito et al.²⁷

In this paper the likelihood function (24) is computed (numerically) via the particle filter described in Appendix B, which takes advantage of Equations (25)–(27). The ML estimate $\hat{\xi}$ of $\xi$ is the value of $\xi$ that maximizes (over the parameter space) the likelihood function. However, the direct maximization of this likelihood function still poses serious convergence issues and carries a heavy computational burden. For this reason, the ML estimates of model parameters are retrieved by using a new expectation maximization particle filter algorithm (described in Appendix A) that significantly reduces these issues. In fact, the adoption of this maximization strategy is crucial for carrying out the misspecification analysis performed in this paper which is based on a Monte Carlo study where the estimation procedure is repeated thousands of times.

## 5 | MISSPECIFICATION ANALYSIS

Two misspecification issues have been addressed: namely, the misspecification of a PGP with a PIGP and the symmetric case of a misspecification of a PIGP with a PGP.

To this aim, we have developed a large Monte Carlo study where three realistic experimental scenarios are simulated by using the setups described in Table 1. As already mentioned in Section 2, the age function is modeled by using the widely adopted power law function $\eta(t) = (t/a)^b$.

In particular, this specific choice of the parameters $\varphi$ and $\nu$ allows (see (21)) to calibrate the error term so that the ratio between the variance of perturbed and hidden processes (i.e., $V\{Z(t)\}/V\{W(t)\}$) is time independent and equal to $1 + 1/(\theta \cdot \varphi) = 1.10$, under all the considered setups. This means that the variance of the error term depends on the size of the measured degradation level and that, due to the presence of the measurement errors, $(\forall t > 0)$ the variance of $Z(t)$ is 10% higher than the variance of $W(t)$.

We have firstly investigated the case where the error term is modeled by using the option (1). Under each setup, we have used the true model (i.e., either the PGP or the PIGP, depending on the misspecification issue of concern) to generate $N_t = 2000$ synthetic datasets. Each dataset consists of $m = 6$ degradation paths, which simulate the evolution of the perturbed degradation levels of as many degrading units. Each path consists of $n_i = 6$ perturbed measurements, taken at equally spaced inspection times $t_1 = 1, t_2 = 2, t_3 = 3, t_4 = 4, t_5 = 5$, and $t_6 = 6$, expressed in time units, that are the same for all the units (that is: $t_{i,j} = t_j = j, \forall i, j, i = 1, \dots, m$ and $j = 1, \dots, 6$). Together with each perturbed measurement

<table><caption>TABLE 1 Setups $\mathcal{A}$, $\mathcal{B}$, and $\mathcal{C}$ used to generate the datasets.</caption><thead><tr><th>Setup</th><th>a</th><th>b</th><th>θ</th><th>φ</th><th>ν</th></tr></thead><tbody><tr><td><i>A</i></td><td>1.25</td><td>1</td><td>1.25</td><td>8</td><td>1</td></tr><tr><td><i>B</i></td><td>1</td><td>1</td><td>1</td><td>10</td><td>1</td></tr><tr><td><i>C</i></td><td>0.5</td><td>1</td><td>0.5</td><td>20</td><td>1</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>FIGURE 1 Pdf of the increments $\Delta W(t_{j-1}, t_j)$ of the considered hidden competing processes under the adopted setups.</figcaption></figure>

we also kept the value of the measured (hidden) degradation level, which is generated (as an intermediate result) by the algorithm that we have adopted to generate the perturbed data. The true values of the hidden (i.e., measured) degradation levels are used to perform comparative analyses.

Note that, under all the considered setups the hidden processes are homogeneous (i.e., the parameter $b$ is always set to 1). Hence, given that they also have independent increments, and that measurement times are equally spaced, the hidden increments $W(t_1), \Delta W(t_1, t_2), \dots, \Delta W(t_5, t_6)$ (where $W(t_1) \equiv \Delta W(t_0, t_1)$) of any considered degradation path are always independent and identically distributed. Moreover, as remarked in Section 2, given the setup, the mean and variance of these increments do not depend on the considered (i.e., either gamma or inverse Gaussian) hidden model.

These convenient setups allow to create, without loss of generality, a direct and easily interpretable link between the risk of misspecification and the parameter $a$ of the age function.

In fact, considered that, both in the case of PGP and PIGP, as $\Delta\eta(t_{j-1}, t_j)$ increases, due to the central limit theorem, the increment $\Delta W(t_{j-1}, t_j)$ tends in distribution to a Gaussian random variable, it is reasonable to expect that, under the considered setups, the similarity between the distributions of the mentioned increments, and consequently the risk of incurring in a misspecification, increases with $\Delta\eta(t_{j-1}, t_j)$ and thus, being $\Delta\eta(t_{j-1}, t_j) = a^{-1}$, decreases with $a$. Indeed, $\Delta\eta(t_{j-1}, t_j)$ is equal to 0.8 under the setup $\mathcal{A}$, to 1 under the setup $\mathcal{B}$, and to 2 under the setup $\mathcal{C}$.

This situation is clearly highlighted in Figure 1, where the solid lines (in red) are the pdfs of the mentioned increments under the PGP and the dashed lines (in blue) are the corresponding pdfs under the PIGP. From Figure 1 it is evident that the similarity between the two increments increases moving from the setup $\mathcal{A}$ to the setup $\mathcal{C}$.

In fact, the setups in Table 1 have been specially designed to simulate scenarios with increasing (i.e., from $\mathcal{A}$ to $\mathcal{C}$) risk of incurring in a model misspecification. The values of $\theta$ have been calibrated so that the mean function of the PGP and PIGP (as well as the mean function of the corresponding hidden models) is the same under all considered setups. The variance of $Z(t)$ and $W(t)$, for any given $t > 0$, decreases moving from the setup $\mathcal{A}$ to the setup $\mathcal{C}$.

Each dataset has been used to perform ML estimates of the parameters $(a, b, \theta, \varphi, \nu)$ of both the PGP and PIGP. Then, dataset by dataset, the AIC has been used to select the model that provides the best fit for that dataset. A misspecification is assumed to occur when the AIC leads to select the wrong model (that is, the PIGP in the case of a dataset generated by using the PGP and the PGP in the case of a dataset generated by using the PIGP). Under each considered setup, results
obtained by performing these analyses have been used to evaluate the percent risk of misspecification as:

$$r_m(\%) = \left( \frac{N_m}{N_t} \right) \cdot 100, \qquad (28)$$

where $N_m$ indicates the number of datasets where a misspecification occurred and $N_t$ is the total number of datasets (i.e., 2000) used to conduct the analyses. For the index $r_m(\%)$ we have computed six values, three values from the datasets generated under the PGP process and three from those generated under the PIGP, which provide (under each setup) the risk of incurring in a misspecification of a PGP with a PIGP and the risk of incurring in the misspecification of a PIGP with a PGP, respectively.

The consequences of a misspecification have been evaluated in terms of its impact on the estimate of the mean remaining useful life, measured by the root mean square error (RMSE) of the ML estimators of the MRUL(6) as:

$$RMSE_{M,d} = \sqrt{\sum_{k=1}^{N_d} \sum_{i=1}^{6} \frac{(\widehat{MRUL}(6)_{M,k,i} - MRUL(6)_{k,i})}{6 \cdot N_d}}, \qquad (29)$$

where:

1. The subscript $d$ of $RMSE_{M,d}$ specifies which are the $N_d$ datasets used to compute the index. In particular, $d = t$ indicates that $RMSE_{M,d}$ is calculated by using all the $N_t$ datasets, $d = \bar{m}$ indicates that it is calculated by using only the $N_{\bar{m}}$ datasets that did not lead to a misspecification, and $d = m$ indicates that it is calculated by considering only the $N_m$ datasets that led to a misspecification.

2. The subscript $M$ of $RMSE_{M,d}$ specifies which is the model used to estimate the MRUL(6). In particular, $M = \text{PGP}$ indicates the estimates are obtained by using the PGP while $M = \text{PIGP}$ indicates that the estimates are obtained by using the PIGP.

3. MRUL(6)$_{k,i}$ is the true value of the MRUL of the unit whose degradation path, up to $t_6 = 6$, is described by the $i$th path ($i = 1, \dots, m$) of the $k$th dataset ($k = 1, \dots, N_d$). In fact, MRUL(6)$_{k,i}$ is computed, path by path, under the true model (i.e., the model used to simulate the data) as follows:

$$\begin{aligned} MRUL(6)_{k,i} &= MRUL(t|W(6)_{k,i} = w(6)_{k,i}) = \int_{0}^{\infty} F_{W(t+\tau)|W(6)_{k,i}}(D|w(6)_{k,i}) \cdot d\tau \\ &= \int_{0}^{\infty} F_{\Delta W(t,t+\tau)}(D - w(6)_{k,i}) \cdot d\tau \end{aligned} \qquad (30)$$

where $F_{\Delta W(t,t+\tau)}(\cdot)$ is either the CDF (7) or (8), depending on the true model (i.e., PGP or PIGP, respectively), with parameters set to the values reported in Table 1 (according to the considered setup), and $w(6)_{k,i}$ is the true (hidden) value of the degradation level of $i$th unit of the $k$th dataset at $t_6 = 6$ (i.e., the true value of $W(6)_{k,i}$), which in this simulation study is known. From each dataset six values of MRUL(6)$_{i,k}$ are obtained.

4. $\widehat{MRUL}(6)_{M,k,i}$ is the ML estimate of the MRUL of the unit whose observed degradation history, up to $t_6 = 6$, is described by the $i$th perturbed path of the $k$th dataset, here denoted as $Z(6)_{k,i}$. $\widehat{MRUL}(6)_{M,k,i}$ is computed as in (23), by using the model $M$ (i.e., either $M = \text{PGP}$ or the $M = \text{PIGP}$) with parameters set to the corresponding ML estimates obtained, under the same model $M$, from the $k$th dataset. In this case the distribution of the RUL used to compute the MRUL is conditional to $Z(6)_{k,i} = Z(6)_{k,i}$. This estimate depends on the dataset, on the model $M$, and on the whole perturbed degradation path of the considered unit (note that both the PGP and PIGP are non-Markovian). From each dataset are obtained twelve values of $\widehat{MRUL}(6)_{M,k,i}$, six under the PGP and six under the PIGP. Hence, $\widehat{MRUL}(6)_{M,k,i}$ differs from the true value MRUL(6)$_{k,i}$ because it is conditional to the perturbed measurements instead than on the true (hidden) degradation level and because model parameters are estimated. Moreover, it is worth to underline that, when $\widehat{MRUL}(6)_{M,k,i}$ is computed by using the wrong model the ML estimates of model parameters are obtained by using the wrong model;

Next, we have repeated the same analyses by modeling the error term according to the option (2).
<table><caption>TABLE 2 Setups A, B, and C used to generate the datasets when $v = 0$.</caption><thead><tr><th>Setup</th><th>a</th><th>b</th><th>θ</th><th>φ</th><th>v</th></tr></thead><tbody><tr><td>A</td><td>1.25</td><td>1</td><td>1.25</td><td>8/3</td><td>0</td></tr><tr><td>B</td><td>1</td><td>1</td><td>1</td><td>10/3</td><td>0</td></tr><tr><td>C</td><td>0.5</td><td>1</td><td>0.5</td><td>20/3</td><td>0</td></tr></tbody></table>

Furthermore, under both the error models (1) and (2), we have also considered the case where $v = 0$. In this latter case, the simulated data used to conduct the Monte Carlo study have been generated by using the setups described in Table 2, which differ from those given in Table 1 for the parameters of the error term only.

In particular, differently than under the setups reported in Table 1, in this case the ratio between the variance of perturbed and hidden processes (i.e., $V\{Z(t)\}/V\{W(t)\}$) is not time independent. In fact, under each setup (i.e., A, B, and C, respectively); here, the value of the parameter $\varphi$ has been selected to set the variance of the error term (that is equal to $V\{\varepsilon(t)\} = 1/\varphi$) to the value that $V\{\varepsilon(t)\}$assumes at $t = 3$, when $\varphi$ and $v$ are those given in Table 1. Note that under the setups given in Table 1 $V\{\varepsilon(t)\}$ increases linearly from $t = 0$, where it is null, to $t = 6$, where it is equal to $6/\varphi$. Thus, adopting the setups given in Table 2 allows to set $V\{\varepsilon(t)\}$ to the value that the same variance assumes in mean over the time interval (0, 6) under the setups given in Table 1. Accordingly, under the setups given in Table 2, given that $V\{\varepsilon(t)\}$ does not depend on $t$, the ratio $V\{Z(t)\}/V\{W(t)\}$ decreases as $t$ increases.

It is worth to remark again that, under the option (2) when $v = 0$, the adopted modeling solution reduces to the classical case where the error term is Gaussian distributed and independent of the measured degradation level.

Here, for the convenience of the readers, we report a scheme of the step-by-step procedure we have used to evaluate the risk and consequences of model misspecification under a given setup and a given true model:

1. select the true model (i.e., either PGP or the PIGP with a given error model) and a setup (i.e., either A, B, or C) from Table 1 or Table 2, depending on $v$;

2. repeat steps 3–6 $N_t = 2000$ times;

3. under the true model, simulate a synthetic dataset. At this step the true hidden degradation values, generated at an intermediate step of the simulation, are also kept;

4. use the EM algorithm and the particle filter described in Appendices A and B to compute the ML estimates of the parameters of both the PGP and the PIGP. This step defines the “estimated processes” (i.e., the PGP and PIGP calibrated by using the ML estimates of model parameters);

5. use the AIC to select the best model. In case of incorrect diagnosis (i.e., if the selected model is not the one used to generate the dataset) a misspecification is assumed to have occurred;

6. path by path (i.e., for any $k = 1, 2, \dots, 2000$ and $i = 1, 2, \dots, 6$) compute $MRUL(6)_{k,i}$ under the model used to generate the data and $MRUL(6)_{M,k,i}$ both under the (estimated) PGP and the PIGP;

7. use the results obtained by these 2000 iterations to compute the indices (28) and (29).

For the sake of comparison, the same analysis has also been performed in the absence of measurement errors (i.e., by assuming that measurements provide exact values of measured degradation levels). In this latter case the competing models are the GP and the IGP and the datasets used to perform the Monte Carlo study have been generated under the hidden models with parameters (i.e., a, b, and $\theta$) set to the values given in Table 1 (which, as mentioned above, coincide with those given in Table 2).

The results of the mentioned analyses are reported in Sections 5.1 and 5.2.

## 5.1 | Risk of incurring in a misspecification

Table 3 reports the values of the percent risk of incurring in a misspecification both when the true model is the PGP and when the true model is the PIGP (i.e., in the symmetric case). The index $r_m(\%)$ is evaluated as in (28). This table refers to the case where the error term is modeled by using the option (1). The setups used to generate the data are those given in Table 1. As mentioned above, $N_t$, $N_{\bar{m}}$, and $N_m$ indicate the total number of datasets used to perform the analysis, the number of datasets that did not lead to a misspecification, and the number of datasets that led to a misspecification, respectively.
<table><caption>TABLE 3 Risk of misspecification when the error term is modeled by using the option (1).</caption><thead><tr><th>True process</th><th>Setup</th><th>N<sub>t</sub></th><th>N<sub>m</sub></th><th>N<sub>m</sub></th><th>r<sub>m</sub>(%)</th></tr></thead><tbody><tr><td rowspan="3">PGP</td><td>A</td><td>2000</td><td>1362</td><td>638</td><td>31.9</td></tr><tr><td>B</td><td>2000</td><td>1252</td><td>748</td><td>37.4</td></tr><tr><td>C</td><td>2000</td><td>1065</td><td>935</td><td>46.8</td></tr><tr><td rowspan="3">PIGP</td><td>A</td><td>2000</td><td>1484</td><td>516</td><td>25.8</td></tr><tr><td>B</td><td>2000</td><td>1470</td><td>530</td><td>26.5</td></tr><tr><td>C</td><td>2000</td><td>1368</td><td>632</td><td>31.6</td></tr></tbody></table>

<table><caption>TABLE 4 Risk of misspecification in the absence of measurement errors.</caption><thead><tr><th>True process</th><th>Setup</th><th>N<sub>t</sub></th><th>N<sub>m</sub></th><th>N<sub>m</sub></th><th>r<sub>m</sub>(%)</th></tr></thead><tbody><tr><td rowspan="3">GP</td><td>A</td><td>2000</td><td>1888</td><td>112</td><td>5.60</td></tr><tr><td>B</td><td>2000</td><td>1806</td><td>194</td><td>9.70</td></tr><tr><td>C</td><td>2000</td><td>1541</td><td>459</td><td>22.9</td></tr><tr><td rowspan="3">IGP</td><td>A</td><td>2000</td><td>1843</td><td>157</td><td>7.85</td></tr><tr><td>B</td><td>2000</td><td>1814</td><td>186</td><td>9.30</td></tr><tr><td>C</td><td>2000</td><td>1665</td><td>335</td><td>16.7</td></tr></tbody></table>

Table 4 reports, by using the same notations, the results of the same analyses in the absence of measurement errors.

Table 3 shows that, under all setups, the risk of incurring in a misspecification of a PGP with a PIGP is about 5% – 10% higher than the one of misspecifying a PIGP by a PGP. As expected, in both cases the risk increases moving from the setup A to the setup C. In the worst-case scenario, represented by the setup C, the risk of a wrong diagnosis when the true model is the PGP is close to 50%, whereas when the true model is the PIGP it is close to 30%.

Results reported in Table 3 also show that the risk of misspecification depends on the setup more when the true model is the PGP than when the true model is the PIGP. An intuitive explanation for this is that passing from the setup A to the setup C the shape of the pdf of the increments of the hidden IGP depicted in Figure 1 changes less than the shape of the pdf of the increments of the GP. Moreover, it seems that, under every setup, when the true process is the PIGP, it is often possible to find a PGP that fits the simulated data in an acceptable manner, while the PIGP more rarely allows to adequately fit PGP data generated under the setups A and B, where the increments $W(t_1), \Delta W(t_1, t_2), \dots, \Delta W(t_5, t_6)$ of the hidden process are gamma distributed with shape parameter (i.e., $\Delta\eta(t_1), \Delta\eta(t_1, t_2), \dots, \Delta\eta(t_5, t_6)$) smaller than 1. In fact, from Figure 1, it is apparent that while the pdf of the increments of the GP obtained under the setup C is relatively similar to the pdfs of the increments of the IGP obtained under the setups A, B, and C, none of pdfs obtained under the IGP is similar to the pdf of the increments of the GP obtained under the setups A and B.

The comparison with the results reported in Table 4 shows that, under all the setups, the presence of measurement errors increases the risk of incurring in a misspecification. Obviously, this result was expected, since perturbed data do not allow to directly verify whether the selected model is actually able to adequately fit the true (hidden) degradation data, being only useful to check the ability of the perturbed model to fit the perturbed measurements.

Table 4 also shows that in the absence of measurement errors the risk of misspecifying a GP with an IGP is close to the one of misspecifying an IGP with a GP. Nonetheless, it seems again that the risk of misspecifying an IGP with a GP depends on the setup less than the risk of misspecifying a GP with an IGP.

Table 5 reports the values of the percent risk of incurring in a misspecification in the case the error term is modeled according to the option (2). The index $r_m(%)$ is evaluated as in (28). The setups used to generate the data are again those given in Table 1.

Results reported in Tables 3 and 5 show that the risk of misspecifying a PGP with a PIGP does not significantly depend on the option adopted to model the error term. Moreover, they also show that, when the true model is the PIGP, results obtained by adopting the error term (2) differ from, and are closer to each other than, those obtained by adopting the error term (1). In fact, it seems that modeling the error term according to the option (2) increases the risk of misspecifying a
<table><caption>TABLE 5 Risk of misspecification when the error term is modeled by using the option (2).</caption><thead><tr><th>True process</th><th>Setup</th><th>N<sub>t</sub></th><th>N<sub>m</sub></th><th>N<sub>m</sub></th><th>r<sub>m</sub>(%)</th></tr></thead><tbody><tr><td rowspan="3">PGP</td><td>A</td><td>2000</td><td>1376</td><td>624</td><td>31.2</td></tr><tr><td>B</td><td>2000</td><td>1256</td><td>744</td><td>37.2</td></tr><tr><td>C</td><td>2000</td><td>1067</td><td>933</td><td>46.7</td></tr><tr><td rowspan="3">PIGP</td><td>A</td><td>2000</td><td>1296</td><td>704</td><td>35.2</td></tr><tr><td>B</td><td>2000</td><td>1292</td><td>708</td><td>35.4</td></tr><tr><td>C</td><td>2000</td><td>1334</td><td>666</td><td>33.3</td></tr></tbody></table>

<table><caption>TABLE 6 Risk of misspecification when the error term is modeled by using the option (1) with $v = 0$.</caption><thead><tr><th>True process</th><th>Setup</th><th>N<sub>t</sub></th><th>N<sub>m</sub></th><th>N<sub>m</sub></th><th>r<sub>m</sub>(%)</th></tr></thead><tbody><tr><td rowspan="3">PGP</td><td>A</td><td>2000</td><td>1307</td><td>693</td><td>34.6</td></tr><tr><td>B</td><td>2000</td><td>1278</td><td>722</td><td>36.1</td></tr><tr><td>C</td><td>2000</td><td>1157</td><td>843</td><td>42.2</td></tr><tr><td rowspan="3">PIGP</td><td>A</td><td>2000</td><td>1459</td><td>541</td><td>27.1</td></tr><tr><td>B</td><td>2000</td><td>1439</td><td>561</td><td>28.1</td></tr><tr><td>C</td><td>2000</td><td>1326</td><td>674</td><td>33.7</td></tr></tbody></table>

PIGP with a PGP and reduces the sensitivity of the mentioned risk on the setup. Indeed, somewhat surprisingly, Table 5 also shows that the risk of misspecifying a PIGP with a PGP under the setup $\mathcal{C}$ is very close to (and even slightly smaller than) the corresponding risk computed under the setups $A$ and $B$. We have carefully checked and confirmed these latter results in various ways, yet we do not have an intuitive explanation for them. In fact, we cannot exclude that the value obtained for $r_m(\%)$ under the setup $\mathcal{C}$ is smaller than those obtained under the other setups only for a matter of numerical accuracy. Indeed, the log-likelihood functions computed under the PIGP and PGP (i.e., those used to calculate the AIC index) in about 6% of cases differ only in the fourth significant figure and in less than 1% of cases only in the fifth significant figure, which the (time-demanding) numerical procedure used to compute the log-likelihood does not always allow to calculate in a sufficiently accurate manner.

As we have already mentioned above, the error term (2) allows the perturbed measurement to assume negative values, a result that in many applications can be unrealistic. About this point, we note that, in the case of the datasets used to produce the results reported in Table 5, when the true model is the PGP we have obtained 989 out of 72,000 negative perturbed measurements under the setup $A$, 607 out of 72,000 negative perturbed measurements under the setup $B$, and 69 out of 72,000 negative perturbed measurements under the setup $\mathcal{C}$. While, when the true model is the PIGP we have obtained 443 negative perturbed measurements under the setup $A$, 227 under the setup $B$, and 17 under the setup $\mathcal{C}$.

Finally, Tables 6 and 7 report the values of the percent risk of incurring in a misspecification in the cases where the error term is modeled by using the setups reported in Table 2 (i.e., in the case $v = 0$). In particular, Table 6 refers to the cases where the error term is modeled by using the option (1), while Table 7 reports the results obtained by using the option (2).

Here, in the case of the datasets used to produce the results reported in Table 7, we have obtained 3390 negative perturbed measurements under the setup $A$, 2687 negative perturbed measurements under the setup $B$, and 990 negative perturbed measurements under the setup $\mathcal{C}$ when the true model is the PGP. Similarly, when the true model is the PIGP, we have obtained 2643 negative perturbed measurements under the setup $A$, 2016 under the setup $B$, and 718 under the setup $\mathcal{C}$.

In general, by comparing Tables 3 and 5 with Tables 6 and 7, it seems that, when $v = 0$ the risk of misspecifying a PGP with a PIGP is similar to the one computed when $v = 1$. Yet, it also seems that, when $v = 0$ the mentioned risk depends on the setup less than when $v = 1$, especially when the error term is modeled by using the option (2). The same effect is also observed when the true model is the PIGP. However, from the same tables, it also seems that the risk of misspecifying the PIGP with a PGP under the setups given in Table 2 is higher than under those given in Table 1.
<table><caption>TABLE 7 Risk of misspecification when the error term is modeled by using the option (2) with $v = 0$.</caption><thead><tr><th>True process</th><th>Setup</th><th>N<sub>t</sub></th><th>N<sub>m</sub></th><th>N<sub>m</sub></th><th>r<sub>m</sub>(%)</th></tr></thead><tbody><tr><td rowspan="3">PGP</td><td>A</td><td>2000</td><td>1247</td><td>753</td><td>37.7</td></tr><tr><td>B</td><td>2000</td><td>1242</td><td>758</td><td>37.9</td></tr><tr><td>C</td><td>2000</td><td>1139</td><td>861</td><td>43.1</td></tr><tr><td rowspan="3">PIGP</td><td>A</td><td>2000</td><td>1216</td><td>784</td><td>39.2</td></tr><tr><td>B</td><td>2000</td><td>1195</td><td>805</td><td>40.3</td></tr><tr><td>C</td><td>2000</td><td>1185</td><td>815</td><td>40.7</td></tr></tbody></table>

It should be emphasized that, strictly speaking, the results obtained by setting $v = 0$ and $v = 1$ are not entirely comparable to each other because, as mentioned before, while under the setups reported in Table 1 the variance of the error term, $\epsilon(t)$, is proportional to the measured degradation level $W(t)$, under the setups given in Table 2 the variance of $\epsilon(t)$ does not depend on $W(t)$. However, a possible intuitive explanation of the difference existing between the results obtained when $v = 0$ and when $v = 1$ can be given by focusing on the importance of the measurement obtained at the first measurement epoch. Indeed, as shown in Figure 1, the shapes of the pdfs of the increments of the hidden GP and IGP mainly differ in the left tail. The same figure also shows that the difference diminishes as $\eta(t)$ increases (i.e., moving from the setup A to the setup C). For the same reasons the considered shapes become more and more similar as $t$ increases, because $\eta(t)$ increases with $t$ and both the gamma and the inverse Gaussian random variables, as $\eta(t)$ increases, tend in distribution to a Gaussian random variable (i.e., to the same Gaussian random variable, because the considered competing hidden processes have identical mean and variance functions). Based on this reasoning, it is reasonable to expect that the first perturbed measurements are the most useful to identify the hidden processes. Consequently, given that $V\{\epsilon(t)\}$ at $t = 1$ is larger when $v = 0$ than when $v = 1$, it results that the perturbed measurements performed at $t = 1$ when $v = 0$ are less useful to identify the hidden model than the corresponding measurements performed when $v = 1$.

## 5.2 | Consequences of incurring in a misspecification

Table 8 reports the values of the index $RMSE_{M,d}$ computed when the error term is modeled according to the option (1). To assess the prognostic abilities of the competing models on a short, medium, and long timespan, we used three different values of the threshold: namely, $w_M = 7.5$, $w_M = 9$, and $w_M = 12$. For the sake of simplicity, this piece of information is not included in the notation of the RMSE, but is directly provided in the headings of the tables. The first column of Table 8 specifies which is the true model. The second column indicates the setup under which the RMSE is computed. Details about these setups are given in Table 1. The third column specifies the model used to compute $MRUL(6)_{M,k,i}$, which is also the one used to obtain the MLEs of model parameters. Finally, the fourth column specifies the datasets used to compute the index. With the same notations as Tables 3–7, the subscripts $t$, $m$, and $\bar{m}$ indicate that the RMSE is evaluated by using all the $N_t$ datasets, only the $N_m$ datasets where a misspecification occurred, and only the $N_{\bar{m}}$ datasets where a misspecification did not occur, respectively.

So, for example, the values reported in the fourth row of Table 8 are the RMSEs of the ML estimator of $MRUL(6)$, computed as in (29), for all the considered thresholds, by using a PIGP with model parameters set to their ML estimates obtained from the $N_m$ datasets generated under a PGP (calibrated according to the setup A) that led to a misspecification.

Table 8 shows that the prognostic abilities of the considered competing perturbed models are quite similar both when the true model is the PGP and when the true model is the PIGP. In addition, as expected, obtained results also show that the difference between the RMSEs obtained under the considered perturbed models decreases moving from the setup A to the setup C. In fact, this is a direct consequence of the circumstance that (as shown by Figure 1) passing from the setup A to the setup C the difference between the pdfs of the already mentioned increments of the hidden processes diminishes.

Table 9 reports the results obtained by performing the same analyses in the absence of measurement errors.

The comparison between the results reported in Table 8 and those reported in Table 9 shows how the presence of measurement errors impacts on the estimation of the $MRUL(6)$. In fact, it is evident that if by one side the presence of measurement errors negatively affects the performances of the ML estimator of $MRUL(6)$ constructed under the right perturbed model (i.e., by using the PGP when the data are generated under a PGP or by using a PIGP when the true model is a PIGP), on the other side it also allows to mitigate the effect of using the wrong model. This is especially clear if one
<table><caption>TABLE 8 RMSE<sub>M,d</sub> of the ML estimators of MRUL(6) when the error term is modeled by using the option (1).</caption><thead><tr><th rowspan="2">True model</th><th rowspan="2">Setup</th><th rowspan="2">M</th><th rowspan="2">d</th><th colspan="3">RMSE<sub>M,d</sub></th></tr><tr><th>w<sub>M</sub> = 7.5</th><th>w<sub>M</sub> = 9</th><th>w<sub>M</sub> = 12</th></tr></thead><tbody><tr><td rowspan="12">PGP</td><td rowspan="6">A</td><td>PGP</td><td><span style="text-decoration: overline;">m</span></td><td>1.67</td><td>2.34</td><td>3.99</td></tr><tr><td>PIGP</td><td><span style="text-decoration: overline;">m</span></td><td>1.53</td><td>2.13</td><td>3.60</td></tr><tr><td>PGP</td><td>m</td><td>1.89</td><td>2.64</td><td>4.44</td></tr><tr><td>PIGP</td><td>m</td><td>1.94</td><td>2.69</td><td>4.46</td></tr><tr><td>PGP</td><td>t</td><td>1.74</td><td>2.44</td><td>4.14</td></tr><tr><td>PIGP</td><td>t</td><td>1.67</td><td>2.32</td><td>3.89</td></tr><tr><td rowspan="6">B</td><td>PGP</td><td><span style="text-decoration: overline;">m</span></td><td>1.38</td><td>1.92</td><td>3.21</td></tr><tr><td>PIGP</td><td><span style="text-decoration: overline;">m</span></td><td>1.32</td><td>1.84</td><td>3.09</td></tr><tr><td>PGP</td><td>m</td><td>1.37</td><td>1.91</td><td>3.20</td></tr><tr><td>PIGP</td><td>m</td><td>1.33</td><td>1.85</td><td>3.07</td></tr><tr><td>PGP</td><td>t</td><td>1.38</td><td>1.92</td><td>3.20</td></tr><tr><td>PIGP</td><td>t</td><td>1.32</td><td>1.84</td><td>3.08</td></tr><tr><td rowspan="6">C</td><td>PGP</td><td><span style="text-decoration: overline;">m</span></td><td>0.76</td><td>1.09</td><td>1.89</td></tr><tr><td>PIGP</td><td><span style="text-decoration: overline;">m</span></td><td>0.74</td><td>1.07</td><td>1.85</td></tr><tr><td>PGP</td><td>m</td><td>0.78</td><td>1.12</td><td>1.92</td></tr><tr><td>PIGP</td><td>m</td><td>0.78</td><td>1.11</td><td>1.91</td></tr><tr><td>PGP</td><td>t</td><td>0.77</td><td>1.10</td><td>1.90</td></tr><tr><td>PIGP</td><td>t</td><td>0.76</td><td>1.09</td><td>1.88</td></tr><tr><td rowspan="12">PIGP</td><td rowspan="6">A</td><td>PGP</td><td><span style="text-decoration: overline;">m</span></td><td>1.56</td><td>2.20</td><td>3.72</td></tr><tr><td>PIGP</td><td><span style="text-decoration: overline;">m</span></td><td>1.52</td><td>2.14</td><td>3.60</td></tr><tr><td>PGP</td><td>m</td><td>1.54</td><td>2.17</td><td>3.64</td></tr><tr><td>PIGP</td><td>m</td><td>1.38</td><td>1.91</td><td>3.14</td></tr><tr><td>PGP</td><td>t</td><td>1.56</td><td>2.19</td><td>3.70</td></tr><tr><td>PIGP</td><td>t</td><td>1.49</td><td>2.08</td><td>3.49</td></tr><tr><td rowspan="6">B</td><td>PGP</td><td><span style="text-decoration: overline;">m</span></td><td>1.26</td><td>1.77</td><td>2.98</td></tr><tr><td>PIGP</td><td><span style="text-decoration: overline;">m</span></td><td>1.21</td><td>1.69</td><td>2.83</td></tr><tr><td>PGP</td><td>m</td><td>1.29</td><td>1.84</td><td>3.15</td></tr><tr><td>PIGP</td><td>m</td><td>1.20</td><td>1.69</td><td>2.87</td></tr><tr><td>PGP</td><td>t</td><td>1.27</td><td>1.79</td><td>3.03</td></tr><tr><td>PIGP</td><td>t</td><td>1.21</td><td>1.69</td><td>2.84</td></tr><tr><td rowspan="6">C</td><td>PGP</td><td><span style="text-decoration: overline;">m</span></td><td>0.74</td><td>1.04</td><td>1.77</td></tr><tr><td>PIGP</td><td><span style="text-decoration: overline;">m</span></td><td>0.73</td><td>1.02</td><td>1.74</td></tr><tr><td>PGP</td><td>m</td><td>0.77</td><td>1.10</td><td>1.90</td></tr><tr><td>PIGP</td><td>m</td><td>0.75</td><td>1.06</td><td>1.83</td></tr><tr><td>PGP</td><td>t</td><td>0.75</td><td>1.06</td><td>1.81</td></tr><tr><td>PIGP</td><td>t</td><td>0.73</td><td>1.04</td><td>1.77</td></tr></tbody></table>
<table><caption>TABLE 9 RMSE<sub>M,d</sub> of the ML estimators of MRUL(6) in the absence of measurement errors.</caption><thead><tr><th rowspan="2">True model</th><th rowspan="2">Setup</th><th rowspan="2">M</th><th rowspan="2">d</th><th colspan="3">RMSE<sub>M,d</sub></th></tr><tr><th>w<sub>M</sub> = 7.5</th><th>w<sub>M</sub> = 9</th><th>w<sub>M</sub> = 12</th></tr></thead><tbody><tr><td rowspan="12">GP</td><td rowspan="6">A</td><td>GP</td><td>m̄</td><td>1.24</td><td>1.79</td><td>3.08</td></tr><tr><td>IGP</td><td>m̄</td><td>4.01 · 10<sup>11</sup></td><td>1.70 · 10<sup>12</sup></td><td>1.70 · 10<sup>13</sup></td></tr><tr><td>GP</td><td>m</td><td>0.94</td><td>1.36</td><td>2.34</td></tr><tr><td>IGP</td><td>m</td><td>2.03</td><td>2.75</td><td>4.38</td></tr><tr><td>GP</td><td>t</td><td>1.23</td><td>1.76</td><td>3.04</td></tr><tr><td>IGP</td><td>t</td><td>3.90 · 10<sup>11</sup></td><td>1.66 · 10<sup>12</sup></td><td>1.62 · 10<sup>13</sup></td></tr><tr><td rowspan="6">B</td><td>GP</td><td>m̄</td><td>0.96</td><td>1.41</td><td>2.47</td></tr><tr><td>IGP</td><td>m̄</td><td>2.92 · 10<sup>14</sup></td><td>2.10 · 10<sup>15</sup></td><td>3.79 · 10<sup>16</sup></td></tr><tr><td>GP</td><td>m</td><td>0.75</td><td>1.09</td><td>1.90</td></tr><tr><td>IGP</td><td>m</td><td>1.22</td><td>1.73</td><td>2.92</td></tr><tr><td>GP</td><td>t</td><td>0.94</td><td>1.39</td><td>2.43</td></tr><tr><td>IGP</td><td>t</td><td>2.77 · 10<sup>14</sup></td><td>1.99 · 10<sup>15</sup></td><td>3.60 · 10<sup>16</sup></td></tr><tr><td rowspan="6">C</td><td>GP</td><td>m̄</td><td>0.51</td><td>0.81</td><td>1.53</td></tr><tr><td>IGP</td><td>m̄</td><td>1.63</td><td>2.54</td><td>4.98</td></tr><tr><td>GP</td><td>m</td><td>0.50</td><td>0.80</td><td>1.50</td></tr><tr><td>IGP</td><td>m</td><td>0.60</td><td>0.93</td><td>1.73</td></tr><tr><td>GP</td><td>t</td><td>0.51</td><td>0.81</td><td>1.53</td></tr><tr><td>IGP</td><td>t</td><td>1.46</td><td>2.28</td><td>4.45</td></tr><tr><td rowspan="12">IGP</td><td rowspan="6">A</td><td>GP</td><td>m̄</td><td>0.96</td><td>1.47</td><td>2.67</td></tr><tr><td>IGP</td><td>m̄</td><td>1.07</td><td>1.61</td><td>2.89</td></tr><tr><td>GP</td><td>m</td><td>0.81</td><td>1.21</td><td>2.12</td></tr><tr><td>IGP</td><td>m</td><td>0.85</td><td>1.25</td><td>2.17</td></tr><tr><td>GP</td><td>t</td><td>0.83</td><td>1.23</td><td>2.16</td></tr><tr><td>IGP</td><td>t</td><td>0.87</td><td>1.28</td><td>2.24</td></tr><tr><td rowspan="6">B</td><td>GP</td><td>m̄</td><td>0.71</td><td>1.12</td><td>2.08</td></tr><tr><td>IGP</td><td>m̄</td><td>0.73</td><td>1.13</td><td>2.08</td></tr><tr><td>GP</td><td>m</td><td>0.70</td><td>1.06</td><td>1.91</td></tr><tr><td>IGP</td><td>m</td><td>1.73</td><td>1.09</td><td>1.95</td></tr><tr><td>GP</td><td>t</td><td>0.70</td><td>1.07</td><td>1.93</td></tr><tr><td>IGP</td><td>t</td><td>0.73</td><td>1.10</td><td>1.96</td></tr><tr><td rowspan="6">C</td><td>GP</td><td>m̄</td><td>0.37</td><td>0.62</td><td>1.20</td></tr><tr><td>IGP</td><td>m̄</td><td>0.38</td><td>0.63</td><td>1.22</td></tr><tr><td>GP</td><td>m</td><td>0.43</td><td>0.69</td><td>1.32</td></tr><tr><td>IGP</td><td>m</td><td>0.43</td><td>0.69</td><td>1.30</td></tr><tr><td>GP</td><td>t</td><td>0.42</td><td>0.68</td><td>1.30</td></tr><tr><td>IGP</td><td>t</td><td>0.42</td><td>0.68</td><td>1.29</td></tr></tbody></table>
<table><caption>TABLE 10 RMSE<sub>M,d</sub> of the ML estimators of MRUL(6), when the error term is modeled by using the option (2).</caption><thead><tr><th rowspan="2">True model</th><th rowspan="2">Setup</th><th rowspan="2">M</th><th rowspan="2">d</th><th colspan="3">RMSE<sub>M,d</sub></th></tr><tr><th>w<sub>M</sub> = 7.5</th><th>w<sub>M</sub> = 9</th><th>w<sub>M</sub> = 12</th></tr></thead><tbody><tr><td rowspan="12">PGP</td><td rowspan="6">A</td><td>PGP</td><td>m̄</td><td>1.73</td><td>2.44</td><td>4.15</td></tr><tr><td>PIGP</td><td>m̄</td><td>1.70</td><td>2.36</td><td>3.94</td></tr><tr><td>PGP</td><td>m</td><td>2.04</td><td>2.85</td><td>4.81</td></tr><tr><td>PIGP</td><td>m</td><td>2.02</td><td>2.79</td><td>4.64</td></tr><tr><td>PGP</td><td>t</td><td>1.83</td><td>2.57</td><td>4.37</td></tr><tr><td>PIGP</td><td>t</td><td>1.80</td><td>2.50</td><td>4.17</td></tr><tr><td rowspan="6">B</td><td>PGP</td><td>m̄</td><td>1.33</td><td>1.85</td><td>3.12</td></tr><tr><td>PIGP</td><td>m̄</td><td>1.29</td><td>1.80</td><td>3.02</td></tr><tr><td>PGP</td><td>m</td><td>1.50</td><td>2.14</td><td>3.67</td></tr><tr><td>PIGP</td><td>m</td><td>1.51</td><td>2.14</td><td>3.68</td></tr><tr><td>PGP</td><td>t</td><td>1.39</td><td>1.96</td><td>3.33</td></tr><tr><td>PIGP</td><td>t</td><td>1.37</td><td>1.93</td><td>3.28</td></tr><tr><td rowspan="6">C</td><td>PGP</td><td>m̄</td><td>0.76</td><td>1.08</td><td>1.86</td></tr><tr><td>PIGP</td><td>m̄</td><td>0.75</td><td>1.08</td><td>1.85</td></tr><tr><td>PGP</td><td>m</td><td>0.75</td><td>1.06</td><td>1.82</td></tr><tr><td>PIGP</td><td>m</td><td>0.75</td><td>1.07</td><td>1.84</td></tr><tr><td>PGP</td><td>t</td><td>0.75</td><td>1.07</td><td>1.84</td></tr><tr><td>PIGP</td><td>t</td><td>0.75</td><td>1.07</td><td>1.84</td></tr><tr><td rowspan="12">PIGP</td><td rowspan="6">A</td><td>PGP</td><td>m̄</td><td>1.59</td><td>2.21</td><td>3.70</td></tr><tr><td>PIGP</td><td>m̄</td><td>1.49</td><td>2.07</td><td>3.45</td></tr><tr><td>PGP</td><td>m</td><td>1.51</td><td>2.14</td><td>3.69</td></tr><tr><td>PIGP</td><td>m</td><td>1.46</td><td>2.05</td><td>3.50</td></tr><tr><td>PGP</td><td>t</td><td>1.56</td><td>2.18</td><td>3.70</td></tr><tr><td>PIGP</td><td>t</td><td>1.48</td><td>2.06</td><td>3.46</td></tr><tr><td rowspan="6">B</td><td>PGP</td><td>m̄</td><td>1.27</td><td>1.78</td><td>2.99</td></tr><tr><td>PIGP</td><td>m̄</td><td>1.18</td><td>1.64</td><td>2.74</td></tr><tr><td>PGP</td><td>m</td><td>1.32</td><td>1.88</td><td>3.26</td></tr><tr><td>PIGP</td><td>m</td><td>1.26</td><td>1.78</td><td>3.03</td></tr><tr><td>PGP</td><td>t</td><td>1.29</td><td>1.82</td><td>3.08</td></tr><tr><td>PIGP</td><td>t</td><td>1.21</td><td>1.69</td><td>2.84</td></tr><tr><td rowspan="6">C</td><td>PGP</td><td>m̄</td><td>0.76</td><td>1.08</td><td>1.86</td></tr><tr><td>PIGP</td><td>m̄</td><td>0.73</td><td>1.04</td><td>1.79</td></tr><tr><td>PGP</td><td>m</td><td>0.77</td><td>1.10</td><td>1.92</td></tr><tr><td>PIGP</td><td>m</td><td>0.74</td><td>1.06</td><td>1.83</td></tr><tr><td>PGP</td><td>t</td><td>0.76</td><td>1.09</td><td>1.88</td></tr><tr><td>PIGP</td><td>t</td><td>0.74</td><td>1.05</td><td>1.80</td></tr></tbody></table>
<table><caption>TABLE 11 RMSE<sub>M,d</sub> of the ML estimators of MRUL(6) when the error term is modeled by using the option (1) with v = 0.</caption><thead><tr><th rowspan="2">True model</th><th rowspan="2">Setup</th><th rowspan="2">M</th><th rowspan="2">d</th><th colspan="3">RMSE<sub>M,d</sub></th></tr><tr><th>w<sub>M</sub> = 7.5</th><th>w<sub>M</sub> = 9</th><th>w<sub>M</sub> = 12</th></tr></thead><tbody><tr><td rowspan="12">PGP</td><td rowspan="6">A</td><td>PGP</td><td>m&#773;</td><td>1.42</td><td>1.97</td><td>3.29</td></tr><tr><td>PIGP</td><td>m&#773;</td><td>1.28</td><td>1.77</td><td>2.93</td></tr><tr><td>PGP</td><td>m</td><td>1.64</td><td>2.31</td><td>3.95</td></tr><tr><td>PIGP</td><td>m</td><td>1.71</td><td>2.42</td><td>4.17</td></tr><tr><td>PGP</td><td>t</td><td>1.50</td><td>2.09</td><td>3.53</td></tr><tr><td>PIGP</td><td>t</td><td>1.45</td><td>2.02</td><td>3.41</td></tr><tr><td rowspan="6">B</td><td>PGP</td><td>m&#773;</td><td>1.27</td><td>1.80</td><td>3.07</td></tr><tr><td>PIGP</td><td>m&#773;</td><td>1.16</td><td>1.64</td><td>2.81</td></tr><tr><td>PGP</td><td>m</td><td>1.35</td><td>1.88</td><td>3.16</td></tr><tr><td>PIGP</td><td>m</td><td>1.21</td><td>1.68</td><td>2.81</td></tr><tr><td>PGP</td><td>t</td><td>1.30</td><td>1.83</td><td>3.1</td></tr><tr><td>PIGP</td><td>t</td><td>1.18</td><td>1.66</td><td>2.81</td></tr><tr><td rowspan="6">C</td><td>PGP</td><td>m&#773;</td><td>0.73</td><td>1.07</td><td>1.92</td></tr><tr><td>PIGP</td><td>m&#773;</td><td>0.72</td><td>1.05</td><td>1.89</td></tr><tr><td>PGP</td><td>m</td><td>0.74</td><td>1.07</td><td>1.89</td></tr><tr><td>PIGP</td><td>m</td><td>0.73</td><td>1.07</td><td>1.88</td></tr><tr><td>PGP</td><td>t</td><td>0.73</td><td>1.07</td><td>1.90</td></tr><tr><td>PIGP</td><td>t</td><td>0.72</td><td>1.06</td><td>1.89</td></tr><tr><td rowspan="12">PIGP</td><td rowspan="6">A</td><td>PGP</td><td>m&#773;</td><td>1.36</td><td>1.90</td><td>3.21</td></tr><tr><td>PIGP</td><td>m&#773;</td><td>1.25</td><td>1.74</td><td>2.91</td></tr><tr><td>PGP</td><td>m</td><td>1.83</td><td>2.64</td><td>4.61</td></tr><tr><td>PIGP</td><td>m</td><td>1.61</td><td>2.31</td><td>3.99</td></tr><tr><td>PGP</td><td>t</td><td>1.50</td><td>2.13</td><td>3.65</td></tr><tr><td>PIGP</td><td>t</td><td>1.36</td><td>1.91</td><td>3.24</td></tr><tr><td rowspan="6">B</td><td>PGP</td><td>m&#773;</td><td>1.11</td><td>1.56</td><td>2.66</td></tr><tr><td>PIGP</td><td>m&#773;</td><td>1.04</td><td>1.45</td><td>2.46</td></tr><tr><td>PGP</td><td>m</td><td>1.40</td><td>2.04</td><td>3.64</td></tr><tr><td>PIGP</td><td>m</td><td>1.24</td><td>1.79</td><td>3.13</td></tr><tr><td>PGP</td><td>t</td><td>1.20</td><td>1.71</td><td>2.97</td></tr><tr><td>PIGP</td><td>t</td><td>1.10</td><td>1.56</td><td>2.67</td></tr><tr><td rowspan="6">C</td><td>PGP</td><td>m&#773;</td><td>0.66</td><td>0.96</td><td>1.7</td></tr><tr><td>PIGP</td><td>m&#773;</td><td>0.63</td><td>0.92</td><td>1.62</td></tr><tr><td>PGP</td><td>m</td><td>0.72</td><td>1.07</td><td>1.94</td></tr><tr><td>PIGP</td><td>m</td><td>0.69</td><td>1.02</td><td>1.82</td></tr><tr><td>PGP</td><td>t</td><td>0.68</td><td>1.00</td><td>1.78</td></tr><tr><td>PIGP</td><td>t</td><td>0.65</td><td>0.95</td><td>1.69</td></tr></tbody></table>
<table><caption>TABLE 12 RMSE<sub>M,d</sub> of the ML estimators of MRUL(6) when the error term is modeled by using the option (2) with v = 0.</caption><thead><tr><th rowspan="2">True model</th><th rowspan="2">Setup</th><th rowspan="2">M</th><th rowspan="2">d</th><th colspan="3">RMSE<sub>M,d</sub></th></tr><tr><th>w<sub>M</sub> = 7.5</th><th>w<sub>M</sub> = 9</th><th>w<sub>M</sub> = 12</th></tr></thead><tbody><tr><td rowspan="12">PGP</td><td rowspan="6">A</td><td>PGP</td><td><span style="text-decoration: overline;">m</span></td><td>1.76</td><td>2.48</td><td>4.23</td></tr><tr><td>PIGP</td><td><span style="text-decoration: overline;">m</span></td><td>1.77</td><td>2.46</td><td>4.14</td></tr><tr><td>PGP</td><td>m</td><td>2.03</td><td>2.95</td><td>5.27</td></tr><tr><td>PIGP</td><td>m</td><td>2.13</td><td>3.08</td><td>5.49</td></tr><tr><td>PGP</td><td>t</td><td>1.87</td><td>2.67</td><td>4.65</td></tr><tr><td>PIGP</td><td>t</td><td>1.91</td><td>2.71</td><td>4.70</td></tr><tr><td rowspan="6">B</td><td>PGP</td><td><span style="text-decoration: overline;">m</span></td><td>1.36</td><td>1.93</td><td>3.29</td></tr><tr><td>PIGP</td><td><span style="text-decoration: overline;">m</span></td><td>1.31</td><td>1.84</td><td>3.12</td></tr><tr><td>PGP</td><td>m</td><td>1.39</td><td>1.99</td><td>3.49</td></tr><tr><td>PIGP</td><td>m</td><td>1.38</td><td>1.98</td><td>3.46</td></tr><tr><td>PGP</td><td>t</td><td>1.37</td><td>1.95</td><td>3.36</td></tr><tr><td>PIGP</td><td>t</td><td>1.33</td><td>1.90</td><td>3.25</td></tr><tr><td rowspan="6">C</td><td>PGP</td><td><span style="text-decoration: overline;">m</span></td><td>0.71</td><td>1.03</td><td>1.82</td></tr><tr><td>PIGP</td><td><span style="text-decoration: overline;">m</span></td><td>0.69</td><td>1.02</td><td>1.79</td></tr><tr><td>PGP</td><td>m</td><td>0.65</td><td>0.95</td><td>1.70</td></tr><tr><td>PIGP</td><td>m</td><td>0.64</td><td>0.94</td><td>1.69</td></tr><tr><td>PGP</td><td>t</td><td>0.68</td><td>1.00</td><td>1.77</td></tr><tr><td>PIGP</td><td>t</td><td>0.67</td><td>0.98</td><td>1.75</td></tr><tr><td rowspan="12">PIGP</td><td rowspan="6">A</td><td>PGP</td><td><span style="text-decoration: overline;">m</span></td><td>1.59</td><td>2.24</td><td>3.85</td></tr><tr><td>PIGP</td><td><span style="text-decoration: overline;">m</span></td><td>1.47</td><td>2.08</td><td>3.55</td></tr><tr><td>PGP</td><td>m</td><td>1.61</td><td>2.31</td><td>4.07</td></tr><tr><td>PIGP</td><td>m</td><td>1.54</td><td>2.21</td><td>3.86</td></tr><tr><td>PGP</td><td>t</td><td>1.59</td><td>2.27</td><td>3.94</td></tr><tr><td>PIGP</td><td>t</td><td>1.50</td><td>2.13</td><td>3.67</td></tr><tr><td rowspan="6">B</td><td>PGP</td><td><span style="text-decoration: overline;">m</span></td><td>1.21</td><td>1.71</td><td>2.94</td></tr><tr><td>PIGP</td><td><span style="text-decoration: overline;">m</span></td><td>1.13</td><td>1.60</td><td>2.73</td></tr><tr><td>PGP</td><td>m</td><td>1.33</td><td>1.92</td><td>3.39</td></tr><tr><td>PIGP</td><td>m</td><td>1.26</td><td>1.82</td><td>3.20</td></tr><tr><td>PGP</td><td>t</td><td>1.26</td><td>1.80</td><td>3.13</td></tr><tr><td>PIGP</td><td>t</td><td>1.19</td><td>1.69</td><td>2.93</td></tr><tr><td rowspan="6">C</td><td>PGP</td><td><span style="text-decoration: overline;">m</span></td><td>0.68</td><td>0.99</td><td>1.79</td></tr><tr><td>PIGP</td><td><span style="text-decoration: overline;">m</span></td><td>0.66</td><td>0.96</td><td>1.73</td></tr><tr><td>PGP</td><td>m</td><td>0.71</td><td>1.05</td><td>1.87</td></tr><tr><td>PIGP</td><td>m</td><td>0.68</td><td>1.00</td><td>1.79</td></tr><tr><td>PGP</td><td>t</td><td>0.69</td><td>1.02</td><td>1.82</td></tr><tr><td>PIGP</td><td>t</td><td>0.67</td><td>0.98</td><td>1.75</td></tr></tbody></table>
compares the results reported in the first half of Table 8, where the true model is the PGP, with those reported in the first half of Table 9, where the true model is the GP. In fact, for example, while the RMSEs reported in the fifth row of Table 9 (obtained under the setup $\mathcal{A}$ by using the right GP) differ by several orders of magnitude from the RMSEs reported in the sixth row of the same table (obtained under the setup $\mathcal{A}$, by using the wrong IGP), the RMSEs reported in the fifth and sixth row of Table 8 (obtained under the right PGP and the wrong PIGP, respectively) are very close to each other. A similar situation is evidenced by the results reported in the eleventh and twelfth row of the same tables (i.e., in the case of the setup $\mathcal{B}$).

Indeed, more specifically, the results reported in Table 9 also show that, even in the absence of measurement error, the use of a good model selection criterion (such as the AIC) allows to greatly mitigate the consequences of adopting the wrong model. In fact, the results reported in the third and fourth row of Table 9 show that the ML estimators constructed under the wrong and right models in the case of the datasets, generated under the setup $\mathcal{A}$, that led to a misspecification (that is, those for which the AIC suggests the use of the wrong model) behave in a very similar manner. This conclusion is further strengthened by the results reported in the first and second row of Table 9 (i.e., those obtained by using only the datasets for which the AIC leads to select the right model) which as those reported in the fifth and sixth row of the same table, differ from each other by several orders of magnitude. In other words, from the results obtained under the setup $\mathcal{A}$ when the true model is the GP (i.e., in the absence of measurement error) it appears that the huge difference existing between the RMSEs computed on the basis of the estimates of the MRUL obtained from all the datasets by using the right GP and the wrong IGP (i.e., the difference existing between the results reported in the fifth row and in the sixth row of Table 9) is due to the (poor) estimates of the MRUL obtained under the wrong IGP from the datasets that do not cause a misspecification. The results reported in the seventh, eighth, ninth, and tenth row of Table 9 show that the same situation occurs also in the case of the setup $\mathcal{B}$.

It is also interesting to note that, as shown by the results reported in the second part of Table 9, when the true model is the IGP the performances of the estimators constructed under the GP and the IGP seem to have very similar performances, regardless of whether a misspecification has occurred or not. As we have already remarked by commenting the results in Table 3, even in the absence of measurement error, the differences existing between the results obtained when the true process is the GP and when it is the IGP are probably due to the circumstance that when the true process is the IGP under all the setups it is possible to find a GP that fits the simulated data in an acceptable manner, whereas the IGP more rarely allows to adequately fit GP data generated under the setups $\mathcal{A}$ and $\mathcal{B}$. In fact, as shown in Figure 1, while the gamma pdf can assume shapes that are similar to those of all the inverse Gaussian pdfs depicted in that figure, the inverse Gaussian pdf cannot assume shapes just as similar to the gamma pdfs corresponding to the setups $\mathcal{A}$ and $\mathcal{B}$.

Table 10 reports the values of the index $RMSE_{M,d}$ computed when the error term is modeled according to the option (2). The $RMSE_{M,d}$ computed in these cases are very similar to those reported in Table 9. That is, it seems that the shape of the distribution of the error term does not significantly affect the properties of the ML estimator of the MRUL.

Finally, Tables 11 and 12 report the values of the index $RMSE_{M,d}$ computed when $v = 0$ and the error term is modeled according to the options (1) and (2), respectively. The $RMSE_{M,d}$ computed in these cases are very similar to those reported in Tables 8 and 10. That is, it seems that also this specific feature of the distribution of the error term (i.e., setting $v = 0$ or $v = 1$) does not significantly affect the properties of the ML estimator of MRUL.

# 6 | CONCLUSIONS

In this paper, we have investigated risk and consequences of the misspecification of a PGP with a PIGP and the symmetrical misspecification problem of a PIGP with a PGP. Two different models are used to describe the error term.

To facilitate the misspecification study, the (hidden) IGP is formulated by using a new parameterization that allows the considered competing models to share the same parameters and the same functional forms of mean and variance functions.

To conduct the analyses we have carried out a large Monte Carlo study where six realistic experimental scenarios, characterized by different misspecification risk, are simulated by using as many model setups. A setup defines the true model, which is either a PGP or a PIGP, depending on the misspecification issue of concern. Hence, with six setups, two hidden models, and two error models, we have defined 24 true models (i.e., 12 PGP and 12 PIGPs). Each true model has been used to generate 2000 synthetic datasets, each one consisting of the degradation paths of six units. Together with each perturbed measurement we also kept the value of the measured (hidden) degradation level, which is generated at
an intermediate step of the algorithm that we have adopted to generate the perturbed data. The true values of the hidden degradation level are used in the analysis for comparative purposes.

From each dataset, we have estimated the parameters of both the competing perturbed processes by using the ML method. Hence, path by path we have estimated the mean remaining useful life by using the model estimated from the dataset which contains the considered path.

The ML estimation of the parameters of the competing perturbed models has been performed by using a new expectation-maximization particle filter algorithm. Dataset by dataset the model selection has been performed by using the Akaike information criterion. A misspecification has been assumed to occur when the Akaike information criterion has led to select the wrong model (i.e., for example, if it has led to select the PGP in the case of a dataset that has been generated under the PIGP).

All the estimates have been used to compute the risk of incurring in a misspecification and to evaluate its effect on the ML estimates of the mean remaining useful life. Finally, for the sake of comparison, the same study has also been repeated in the absence of measurement errors (that is, by considering as competing models a GP and an IGP).

Obtained results demonstrate that, when the error term is modeled by adopting the option (1) (i.e., by using the 3-parameter inverse gamma distribution) the risk of incurring in a misspecification (as expected) is significantly influenced by the value of the shape parameter of the degradation increment of the hidden process (be it gamma or inverse Gaussian) defined between successive measurement times. Furthermore, and more specifically, it also seems that the risk of misspecifying a PGP by a PIGP depends on the setup more than the risk of misspecifying a PIGP with a PGP. In fact, while the PGP often allows to fit PIGP data in an acceptable manner independently of the setup, the PIGP more rarely allows to adequately fit PGP data generated under the setups $\mathcal{A}$ and $\mathcal{B}$ (i.e., cases where the increments of the hidden process between successive inspection epochs are gamma distributed with shape parameter smaller than or equal to 1).

The results obtained by modeling the error term according to the option (2) (i.e., by using a Gaussian distribution) indicate that the shape of the error term does not significantly influence the risk of misspecifying a PGP with a PIGP. Contrarily, it seems that modeling the error term according to the option (2) increases the risk of misspecifying a PIGP with a PGP and reduces the sensitivity of the mentioned risk on the setup.

Finally, the results obtained by setting $\nu = 0$ indicate that also this modeling option does not significantly affect the risk of misspecifying a PGP with a PIGP. Yet, it also seems that, when $\nu = 0$ the mentioned risk depends on the setup less than when $\nu = 1$, especially when the error term is modeled by using the option (2). The same effect is also observed when the true model is the PIGP. However, in this case, it also seems that when $\nu = 0$ the risk of misspecifying the PIGP with a PGP is higher than when $\nu = 1$.

Obtained results also clearly show that the presence of measurement errors significantly increases the risk of selecting the wrong model. Obviously, this result was expected, since perturbed data do not allow to directly check whether the selected model is actually able to adequately fit the real (hidden) degradation process, being only useful to check the ability of the perturbed model to fit the perturbed measurements.

With respect to the consequences produced by a misspecification, obtained results show that the maximum likelihood estimators of the (perturbed measurements based) mean remaining useful life constructed under the competing perturbed processes, irrespectively of the model used to describe the error term, have very similar performances, regardless of whether a misspecification has occurred. Moreover, the comparison of the results obtained in the presence of measurement errors with those obtained in its absence shows that if on one side the presence of measurement errors negatively affects the performances of the maximum likelihood estimator of the mean remaining useful life constructed under the right perturbed model, it also mitigates (with respect to the case where measurement errors are absent) the consequences of using the wrong model, especially when the true hidden process is the GP. However, the results obtained in the absence of measurement errors also show that, when the true model is the GP, even in experimental situations where the consequences of the use of the wrong model could be extremely severe, adopting an appropriate model selection procedure (such as the Akaike information criterion), allows to hugely mitigate the effect of a misspecification.

## ACKNOWLEDGMENTS

The authors thank the associate editor and the anonymous reviewers for their precious comments and suggestions that allowed to significantly improve the quality of the paper.
## FUNDING INFORMATION

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research activity was supported by Università di Napoli Federico II in the frame of the international agreement between Dipartimento di Ingegneria Industriale and Polytech Angers (codice identificativo 000011-ALTRI-2021-M-GIORGIO_001_001, prot. 40162 del 20/4/2021), by the Université Franco Italienne within the frame of the chapitre 2 of Vinci project (subvention N° C2-221), and by the WISE project of the Région Pays de la Loire (France).

## CONFLICT OF INTEREST STATEMENT

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

## DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## ORCID

Nicola Esposito https://orcid.org/0000-0001-6392-1294

Agostino Mele https://orcid.org/0000-0001-6796-8508

Bruno Castanier https://orcid.org/0000-0002-3735-3331

Massimiliano Giorgio https://orcid.org/0000-0002-5348-5289

## REFERENCES

1. Abdel-Hameed MA. A gamma wear process. *IEEE Trans Reliab.* 1975;24(2):152-154. 10.1109/TR.1975.5215123

2. van Noortwijk JM. A survey of the application of gamma processes in maintenance. *Reliab Eng Syst Saf.* 2009; 94(1): 2-21. 10.1016/j.ress.2007.03.019

3. Wu Z, Wang Z, Feng Q, Sun B, Qian C, Ren Y, Jiang X. A gamma process-based prognostics method for CCT shift of high-power white LEDs. *IEEE Trans Electron Dev.* 2018;65(7):2909-2916. 10.1109/TED.2018.2835651

4. Wang X, Wang BX, Hong Y, Jiang PH. Degradation data analysis based on gamma process with random effects. *Eur J Oper Res.* 2021;292(3):1200-1208. 10.1016/j.ejor.2020.11.036

5. Wasan MT. On an inverse Gaussian process. *Scandinavian Actuarial Journal*. 1968;1968(1-2):69-96. doi:10.1080/03461238.1968.10413264

6. Wang X, Xu D. An inverse Gaussian process model for degradation data. *Dent Tech.* 2010;52(2):188-197. doi:10.1198/TECH.2009.08197

7. Ye ZS, Chen N. The inverse Gaussian process as a degradation model. *Dent Tech.* 2014;56(3):302-311. doi:10.1080/00401706.2013.830074

8. Peng CY. Inverse Gaussian processes with random effects and explanatory variables for degradation data. *Dent Tech.* 2015;57(1):100-111. doi:10.1080/00401706.2013.879077

9. Rodríguez-Picón LA, Rodríguez-Picón AP, Alvarado-Iniesta A. Degradation modeling of 2 fatigue-crack growth characteristics based on inverse Gaussian processes: A case study. *Appl Stochast Models Bus Indus.* 2019;35(3):504-521. doi:10.1002/asmb.2329

10. Morita LH, Tomazella VL, Ferreira PH, Ramos PL, Balakrishnan N, Louzada F. Optimal burn-in policy based on a set of cutoff points using mixture inverse Gaussian degradation process and copulas. *Appl Stoch Models Bus Indus.* 2021;37(3):612-627. doi:10.1002/asmb.2601

11. Ye ZS, Xie M. Stochastic modelling and analysis of degradation for highly reliable products. *Appl Stoch Models Bus Indus.* 2015;31(1):16-32. doi:10.1002/asmb.2063

12. Kahle W, Mercier S, Paroissin C. Degradation Processes in Reliability. ISTE Ltd and. John Wiley & Sons; 2016.

13. Lu D, Pandey MD, Xie WC. An efficient method for the estimation of parameters of stochastic gamma process from noisy degradation measurements. *Proc Instit Mech Eng, O: J Risk Reliab.* 2013;227(4):425-433. doi:10.1177/1748006X13477008

14. Zhang M, Revie M. Model selection with application to gamma process and inverse Gaussian process. In: Walls L, Revie M, Bedford T, eds. *Risk, Reliability & Safety*. CRC Press, Taylor & Francis Group; 2016:1505-1511. doi:10.1201/9781315374987-226

15. Tseng ST, Yao YC. Misspecification analysis of a gamma with inverse Gaussian degradation processes. In: Chen DG, Lio Y, Ng H, Tsai TR, eds. *Statistical Modeling for Degradation Data ICSA Book Series in Statistics*. Springer; 2017:193-208. doi:10.1007/978-981-10-5194-4_10

16. Giorgio M, Mele A, Pulcini G. perturbed gamma degradation process with degradation dependent non-Gaussian measurement errors. *Appl Stoch Models Bus Indus.* 2019;35(2):198-210. doi:10.1002/asmb.2377

17. Le Son K, Fouladirad M, Barros A. Remaining useful lifetime estimation and noisy gamma deterioration process. *Reliab Eng Syst Saf.* 2016;149:76-87. doi:10.1016/j.ress.2015.12.016

18. Bordes L, Paroissin C, Salami A. Parametric inference in a perturbed gamma degradation process. *Commun Stat – Theory methods.* 2015;45(9):2730-2747. doi:10.1080/03610926.2014.892133

19. Hao S, Yang J, Berenguer C. Degradation analysis based on an extended inverse Gaussian process model with skew-normal random effects and measurement errors. *Reliab Eng Syst Saf.* 2019;189:261-270. doi:10.1016/j.ress.2019.04.031
20. Chen X, Ji G, Sun X, Li Z. Inverse Gaussian-based model with measurement errors for degradation analysis. Proc Instit Mech Eng, O: J Risk Reliab. 2019;233(6):1086-1098. 10.1177/1748006X19860682

21. Pulcini G. A perturbed gamma process with statistically dependent measurement errors. Reliab Eng Syst Saf. 2016;152:296-306. 10.1016/j.ress.2016.03.024

22. Akaike H. A new look at the statistical model identification. IEEE Trans Autom Contr. 1974;19(6):716-723. doi:10.1109/TAC.1974.1100705

23. Hao S, Yang J, Berenguer C. A perturbed inverse Gaussian process model with time varying variance-to-mean ratio. In: Beer M, Zio E, eds. Proceedings of the 29th European Safety and Reliability Conference (ESREL2019). Research Publishing Services; 2019:739-745. doi:10.3850/978-981-11-2724-3_0144-cd

24. Sun B, Li Y, Wang Z, Ren Y, Feng Q, Yang, D. An improved inverse Gaussian process with random effects and measurement errors for RUL prediction of hydraulic piston pump. Measurement. 2021;173:108604. 10.1016/j.measurement.2020.108604

25. Bautista L, Castro IT, Landesa L. Condition-based maintenance for a system subject to multiple degradation processes with stochastic arrival intensity. Eur J Oper Res. 2022;302(2):560-574. 10.1016/j.ejor.2022.01.004

26. Bismut E, Pandey MD, Straub D. Reliability-based inspection and maintenance planning of a nuclear feeder piping system. Reliab Eng Syst Saf. 2022;224:108521. 10.1016/j.ress.2022.108521

27. Esposito N, Mele A, Castanier B, Giorgio M. A new gamma degradation process with random effect and state-dependent measurement error. Proc Instit Mech Eng, O: J Risk Reliab. 2023;237(5):868-885. 10.1177/1748006X211067299

28. Dempster AP, Laird NM, Rubin DB. Maximum likelihood from incomplete data via the EM algorithm. JR Statist Soc Ser B. 1977; 3981:1-22. 10.1111/j.2517-6161.1977.tb01600.x

29. Doucet A, Johansen AM. A tutorial on particle filtering and smoothing: Fifteen years later. In: Crisan D, Rozovskii, B, eds. The Oxford Handbook on Nonlinear Filtering. Oxford University Press; 2011:656-704.

**How to cite this article:** Esposito N, Mele A, Castanier B, Giorgio M. Misspecification analysis of gamma- and inverse Gaussian-based perturbed degradation processes. *Appl Stochastic Models Bus Ind.* 2024;40(3):640-667. doi: 10.1002/asmb.2824

## APPENDIX A. THE EM ALGORITHM

The EM algorithm is a general approach for the iterative computation of ML estimates in the presence of missing values and/or incomplete observations (see, e.g., Dempster et al.²⁸). It consists of a two-step sequence, an expectation step (E-step) and a maximization step (M-step), that is iterated until a predefined convergence condition is reached.

In this paper, we treat as missing data the (unknown) values $\mathbf{w} = \{\mathbf{w}_1, \dots, \mathbf{w}_m\}$ of the (unobservable) true degradation levels $\mathbf{W} = \{\mathbf{W}_1, \dots, \mathbf{W}_m\}$ of the $m$ units at the measurement times, where $\mathbf{w}_i = \{w_{i,1}, \dots, w_{i,n_i}\}$ is the realization of $\mathbf{W}_i = \{W_{i,1}, \dots, W_{i,n_i}\}$. Obviously, the observed data consist in the realizations $\mathbf{z} = \{\mathbf{z}_1, \dots, \mathbf{z}_m\}$ of the related perturbed measurements $\mathbf{Z} = \{\mathbf{Z}_1, \dots, \mathbf{Z}_m\}$, where $\mathbf{z}_i = \{z_{i,1}, \dots, z_{i,n_i}\}$ is the realization of $\mathbf{Z}_i = \{Z_{i,1}, \dots, Z_{i,n_i}\}$.

Under this setup, the complete likelihood (i.e., the likelihood formulated considering all data, available and missing) can be formulated as follows:

$$L(\xi; \mathbf{z}, \mathbf{w}) = \prod_{i=1}^{m} f_{\mathbf{z}_i|\mathbf{w}_i}(\mathbf{z}_i|\mathbf{w}_i) \cdot f_{\mathbf{w}_i}(\mathbf{w}_i) \cdot d\mathbf{w}_i,$$

where

$$f_{\mathbf{z}_i|\mathbf{w}_i}(\mathbf{z}_i|\mathbf{w}_i) = \prod_{j=1}^{n_i} f_{\mathbf{z}_{i,j}|\mathbf{w}_{i,j}}(\mathbf{z}_{i,j}|\mathbf{w}_{i,j}),$$

$$f_{\mathbf{w}_i}(\mathbf{w}_i) = \prod_{j=1}^{n_i} f_{\Delta\mathbf{w}_{i,j}}(\Delta\mathbf{w}_{i,j}),$$

where $f_{\Delta\mathbf{w}_{i,j}}(\Delta\mathbf{w}_{i,j})$ is either the pdf in (2) or in (3), depending on the hidden model, and $f_{\mathbf{z}_{i,j}|\mathbf{w}_{i,j}}(\mathbf{z}_{i,j}|\mathbf{w}_{i,j})$ is either the pdf in (12) or in (16), depending on the model used to describe the error term.
Thus, the complete likelihood function can be expressed as:

$$L(\xi; \mathbf{z}, \mathbf{w}) = \prod_{i=1}^{m} \prod_{j=1}^{n_i} f_{Z_{i,j}|W_{i,j}}(z_{i,j} | w_{i,j}) \cdot \prod_{i=1}^{m} \prod_{j=1}^{n_i} f_{\Delta W_{i,j}}(\Delta w_{i,j}) \quad (A1)$$

and the corresponding log-likelihood function can be written as:

$$l(\xi; \mathbf{z}, \mathbf{w}) = \sum_{i=1}^{m} \sum_{j=1}^{n_i} \ln(f_{Z_{i,j}|W_{i,j}}(z_{i,j} | w_{i,j})) + \sum_{i=1}^{m} \sum_{j=1}^{n_i} \ln(f_{\Delta W_{i,j}}(\Delta w_{i,j})). \quad (A2)$$

At the (h + 1)th iteration, the E-step consists in computing the conditional mean $Q(\xi|\xi^{(h)})$:

$$Q(\xi|\xi^{(h)}) = E\{l(\xi; \mathbf{z}, \mathbf{W})|\mathbf{Z} = \mathbf{z}, \xi^{(h)}\} \quad (A3)$$

of the complete log-likelihood with respect to W by using its conditional distribution, given $\mathbf{Z} = \mathbf{z}$, with parameter vector $\xi$ set to its tentative estimate available after the $h$th iteration, denoted by $\xi^{(h)}$. Under the perturbed models described in Section 2, the conditional distribution of W given $\mathbf{Z} = \mathbf{z}$ does not allow for a closed form expression. Hence, the E-step is performed by using the particle filter-based approach described in Appendix B.

The corresponding M-step consists in finding the new tentative estimate $\xi^{(h+1)}$ that, by definition, is the value of $\xi$ that maximizes the function $Q(\xi|\xi^{(h)})$.

In the case of the considered perturbed models, the advantages of adopting the EM approach mainly lie in this latter step, because the conditional expectation $Q(\xi|\xi^{(h)}) = Q(a, b, \theta, \varphi, v|\xi^{(h)})$ splits into the sum of the functions $Q_E(\varphi, v|\xi^{(h)})$ and $Q_H(a, b, \theta|\xi^{(h)})$, where $Q_E(\varphi, v|\xi^{(h)})$:

$$Q_E(\varphi, v|\xi^{(h)}) = \sum_{i=1}^{m} \sum_{j=1}^{n_i} E\{\ln(f_{Z_{i,j}|W_{i,j}}(z_{i,j} | W_{i,j}))|Z_i = z_i, \xi^{(h)}\}$$

given $\xi^{(h)}$, depends only on the parameters $\varphi$ and $v$ of the perturbing term, and $Q_H(a, b, \theta|\xi^{(h)})$

$$Q_H(a, b, \theta|\xi^{(h)}) = \sum_{i=1}^{m} \sum_{j=1}^{n_i} E\{\ln(f_{\Delta W_{i,j}}(\Delta W_{i,j}))|Z_i = z_i, \xi^{(h)}\}$$

depends only on the parameters of the hidden model.

Obviously, the functional form of $Q_E(\varphi, v|\xi^{(h)})$ and $Q_H(a, b, \theta|\xi^{(h)})$ also depends on the model adopted to describe the error term and hidden process, respectively.

In fact, when the error term is modeled by using the option (1) (i.e., by using the inverse gamma conditional pdf (9)) $Q_E(\varphi, v|\xi^{(h)})$ reduces to:

$$Q_E^1(\varphi, v|\xi^{(h)}) = -\sum_{i=1}^{m} \sum_{j=1}^{n_i} \frac{E\{\alpha(W_{i,j})|Z_i = z_i, \xi^{(h)}\}}{Z_{i,j}} + \sum_{i=1}^{m} \sum_{j=1}^{n_i} E\{\beta(W_{i,j}) \cdot \ln(\alpha(W_{i,j}))|Z_i = z_i, \xi^{(h)}\} \\ - \sum_{i=1}^{m} \sum_{j=1}^{n_i} E\{\beta(W_{i,j})|Z_i = z_i, \xi^{(h)}\} \cdot \ln(z_{i,j}) - \sum_{i=1}^{m} \sum_{j=1}^{n_i} E\{\ln(\Gamma(\beta(W_{i,j})))|Z_i = z_i, \xi^{(h)}\}. \quad (A4)$$

whereas, when the error term is modeled by using the option (2) (i.e., by using the Gaussian conditional pdf (15)), $Q_E(\varphi, v|\xi^{(h)})$ reduces to:

$$Q_E^2(\varphi, v|\xi^{(h)}) = -\frac{n_t}{2} \cdot \ln(2\pi) + \frac{n_t}{2} \cdot \ln(\varphi) - \frac{\nu}{2} \cdot \sum_{i=1}^{m} \sum_{j=1}^{n_i} E\{\ln(W_{i,j})|Z_i = z_i, \xi^{(h)}\} \\ - \frac{\varphi}{2} \cdot \sum_{i=1}^{m} \sum_{j=1}^{n_i} E\left\{\frac{(z_{i,j} - W_{i,j})^2}{W_{i,j}^{\nu}}\middle|Z_i = z_i, \xi^{(h)}\right\}, \quad (A5)$$
where the superscripts 1 and 2 in $Q_E^1(\varphi, v|\xi^{(h)})$ and $Q_E^2(\varphi, v|\xi^{(h)})$ indicate the option adopted to model the error term, in agreement with the symbols introduced in Section 2.

Similarly, from (2), under the GP, the function $Q_H(a, b, \theta|\xi^{(h)})$, takes the form:

$$Q_H^G(a, b, \theta|\xi^{(h)}) = -\frac{\sum_{i=1}^m E\{W_{i,n_i} | \mathbf{Z}_i = \mathbf{z}_i, \xi^{(h)}\}}{\theta} + \sum_{i=1}^m \sum_{j=1}^{n_i} \left[ \left(\frac{t_{i,j}}{a}\right)^b - \left(\frac{t_{i,j-1}}{a}\right)^b \right] \cdot E\{\ln(\Delta W_{i,j}) | \mathbf{Z}_i = \mathbf{z}_i, \xi^{(h)}\} \\ - \ln(\theta) \cdot \sum_{i=1}^m \left(\frac{t_{i,n_i}}{a}\right)^b - \sum_{i=1}^m \sum_{j=1}^{n_i} \ln\left(\Gamma\left(\left(\frac{t_{i,j}}{a}\right)^b - \left(\frac{t_{i,j-1}}{a}\right)^b\right)\right), \quad (A6)$$

whereas, from (3), under the inverse Gaussian process it reduces to

$$Q_H^{IG}(a, b, \theta|\xi^{(h)}) = -\frac{n_t}{2} \cdot \ln(2 \cdot \pi) + \frac{n_t}{2} \cdot \ln(\theta) - n_t \cdot b \cdot \ln(a) - \frac{3}{2} \cdot \sum_{i=1}^{m} \sum_{j=1}^{n_i} E\{\ln(\Delta W_{i,j}) | \mathbf{Z}_i = \mathbf{z}_i, \xi^{(h)}\} \\ - \frac{1}{2 \cdot \theta} \cdot \sum_{i=1}^{m} E\{W_{i,n_i} | \mathbf{Z}_i = \mathbf{z}_i, \xi^{(h)}\} - \frac{\theta}{2 \cdot a^{2 \cdot b}} \cdot \sum_{i=1}^{m} \sum_{j=1}^{n_i} E\left\{\frac{\left(t_{i,j}^b - t_{i,j-1}^b\right)^2}{\Delta W_{i,j}} \bigg| \mathbf{Z}_i = \mathbf{z}_i, \xi^{(h)}\right\} \\ + \frac{1}{a^b} \cdot \sum_{i=1}^{m} t_{i,n_i}^b + \sum_{i=1}^{m} \sum_{j=1}^{n_i} \ln(t_{i,j}^b - t_{i,j-1}^b). \quad (A7)$$

where $n_t = \sum_{i=1}^m n_i$ and the superscripts G and IG in $Q_H^G(a, b, \theta|\xi^{(h)})$ and $Q_H^{IG}(a, b, \theta|\xi^{(h)})$ indicate that the hidden process used to compute $Q_H(a, b, \theta|\xi^{(h)})$ are the gamma and the inverse Gaussian, respectively. Both in (A6) and (A7) the age function is modeled as $\eta(t) = (t/a)^b$.

The expressions in (A5), (A6) and (A7) can be further simplified. In fact, from (A5), by solving with respect to $\varphi$ the equation:

$$\frac{\partial Q_E^2(\varphi, v|\xi^{(h)})}{\partial \varphi} = 0$$

the explicit form

$$\tilde{\varphi}(v|\xi^{(h)}) = \frac{n_t}{\sum_{i=1}^{m} \sum_{j=1}^{n_i} E\left\{\frac{(z_{i,j} - W_{i,j})^2}{W_{i,j}^\gamma} \bigg| \mathbf{Z}_i = \mathbf{z}_i, \xi^{(h)}\right\}} \quad (A8)$$

is obtained for the value $\tilde{\varphi}(v|\xi^{(h)})$ that maximizes (A6) with respect to $\varphi$ when $v$ is set to the values indicated in the parentheses on the left side of the equation. By exploiting this result, $v^{(h+1)}$ can be obtained by numerically maximizing this 1-parameter function:

$$Q_E^2(v|\tilde{\varphi}, \xi^{(h)}) = -\frac{n_t}{2} \cdot \ln(2\pi) + \frac{n_t}{2} \cdot \ln(\tilde{\varphi}) - \frac{v}{2} \cdot \sum_{i=1}^{m} \sum_{j=1}^{n_i} E\{\ln(W_{i,j}) | \mathbf{Z}_i = \mathbf{z}_i, \xi^{(h)}\} - \frac{n_t}{2}.$$

Then, $\varphi^{(h+1)}$ can be obtained from (A8) by setting $v = v^{(h+1)}$.

Similarly, from (A6), by solving with respect to $\theta$ the equation:

$$\frac{\partial Q_H^G(a, b, \theta|\xi^{(h)})}{\partial \theta} = 0,$$

the explicit form

$$\tilde{\theta}(a, b|\xi^{(h)}) = \frac{\sum_{i=1}^{m} E\{W_{i,n_i} | \mathbf{Z}_i = \mathbf{z}_i, \xi^{(h)}\}}{\sum_{i=1}^{m} \left(\frac{t_{i,n_i}}{a}\right)^b} \quad (A9)$$
is obtained for the value $\tilde{\theta}(a, b)$ that maximizes (A6) with respect to $\theta$ when $a$ and $b$ are set to the values indicated in the parentheses on the left side of the equation. By exploiting this result, $a^{(h+1)}$ and $b^{(h+1)}$ can be obtained by numerically maximizing the following 2-parameter function:

$$Q_H^G(a, b|\tilde{\theta}, \xi^{(h)}) = -\sum_{i=1}^m \sum_{j=1}^{n_i} \ln\left(\Gamma\left(\left(\frac{t_{i,j}}{a}\right)^b - \left(\frac{t_{i,j-1}}{a}\right)^b\right)\right) + \sum_{i=1}^m \sum_{j=1}^{n_i} \left[\left(\frac{t_{i,j}}{a}\right)^b - \left(\frac{t_{i,j-1}}{a}\right)^b\right] \cdot E\left\{\ln(\Delta W_{i,j}) | Z_i = z_i, \xi^{(h)}\right\} \\ + \left[1 + \ln\left(\frac{\sum_{i=1}^m E\{W_{i,n_i} | Z_i = z_i, \xi^{(h)}\}}{\sum_{i=1}^m \left(\frac{t_{i,n_i}}{a}\right)^b}\right)\right] \cdot \sum_{i=1}^m \left(\frac{t_{i,n_i}}{a}\right)^b$$

Then, $\theta^{(h+1)}$ can be obtained from (A9) by setting $a = a^{(h+1)}$ and $b = b^{(h+1)}$.

Finally, in the case of the IGP, from (A7), by solving with respect to $\theta$ and $a$ the system of equations:

$$\left\{ \begin{aligned} \frac{\partial Q_H^{IG}(a, b, \theta | \xi^{(h)})}{\partial \theta} &= 0 \\ \frac{\partial Q_H^{IG}(a, b, \theta | \xi^{(h)})}{\partial a} &= 0 \end{aligned} \right.$$

the following explicit forms are obtained:

$$\tilde{\theta}(b|\xi^{(h)}) = \frac{\left[ \frac{\sum_{i=1}^{m} E\{W_{i,n_i} | Z_i = z_i, \xi^{(h)}\}}{\sum_{i=1}^{m} t_{i,n_i}^b} \right]^2 \cdot \sum_{i=1}^{m} \sum_{j=1}^{n_i} E\left\{ \frac{\left( \frac{t_{i,j}^b - t_{i,j-1}^b}{\Delta W_{i,j}} \right)^2}{\Delta W_{i,j}} \middle| Z_i = z_i, \xi^{(h)} \right\} - \sum_{i=1}^{m} E\{W_{i,n_i} | Z_i = z_i, \xi^{(h)}\}}{n_t} \quad (A10)$$

$$\tilde{a}(b|\xi^{(h)}) = \left[ \frac{\frac{\sum_{i=1}^{m} E\{W_{i,n_i} | Z_i = z_i, \xi^{(h)}\}}{\sum_{i=1}^{m} t_{i,n_i}^b} \cdot \sum_{i=1}^{m} \sum_{j=1}^{n_i} E\left\{ \frac{\left( \frac{t_{i,j}^b - t_{i,j-1}^b}{\Delta W_{i,j}} \right)^2}{\Delta W_{i,j}} \middle| Z_i = z_i, \xi^{(h)} \right\} - \sum_{i=1}^{m} t_{i,n_i}^b}{n_t} \right]^{\frac{1}{b}} \quad (A11)$$

for the values $\tilde{\theta}(b)$ and $\tilde{a}(b)$ of $\theta$ and $a$ that (jointly) maximize $Q_H^G(a, b, \theta | \xi^{(h)})$ when the parameter $b$ is set to the value reported in the parentheses. By exploiting these results, $b^{(h+1)}$ can be obtained by numerically maximizing the 1-parameter function that is obtained by setting $\theta$ and $a$ in (A7) to $\tilde{\theta}(b)$ and $\tilde{a}(b)$ (we do not report here the resulting function because it is rather bulky). Then, $\theta^{(h+1)}$ and $a^{(h+1)}$ can be obtained by evaluating (A10) and (A11) at $b = b^{(h+1)}$.

Therefore, when the error term is modeled according to the option (1), the M-step results in maximizing (numerically) two 2-parameter functions in the case of the GP, and one 2-parameter function and one 1-parameter function in the case of the IGP. Whereas, when the error term is modeled by using the option (2), the M-step results in maximizing (numerically) one 2-parameter function and one 1-parameter function in the case of the GP, and two 1-parameter functions in the case of the IGP.

It is clear that maximizing these functions is far easier than directly maximizing (numerically) the 5-parameter likelihood function $L(\xi; z)$ in (24).

The convergence condition is satisfied when the absolute difference:

$$\left| \frac{\ln(L(\xi^{(h)}; z))}{\ln(L(\xi^{(h+1)}; z))} - 1 \right|$$

drops below a predetermined value. In case the convergence condition is met at the $h$th iteration, $\xi^{(h+1)}$ is assumed to be the ML estimate of the parameter vector $\xi$. Otherwise, the procedure continues. The iterative algorithm is initialized by assigning at the first iteration a tentative estimate, say $\xi^{(0)}$.
APPENDIX B. THE PARTICLE FILTER.

In this paper, we use a particle filter algorithm (see, e.g., Doucet and Johansen²⁹) to generate random samples from the (joint) conditional distribution of W given Z = z. The data generated by using this procedure are then used to compute (empirically) the likelihood function, the distribution of the RUL, and the conditional expectations requested in the M-step of the EM algorithm described in Appendix A.

The procedure described below allows to obtain N pseudorandom realizations of $W_i$ given $Z_i = z_i$. Thus, to generate pseudorandom realizations of W given $Z = z$ it is necessary to replicate its use for any $i = 1, ..., m$.

The method consists in iterating $n_i$ times the 2-step sequence described below. To apply the procedure it is necessary to assign a value to the parameter vector $\xi$. The conditional pdf of $\Delta W(t_{i,j-1}, t_{i,j})$ mentioned in the prediction step is either the pdf in (2) or the one in (3), depending on the hidden process model (i.e., gamma or inverse Gaussian, respectively). Analogously, the conditional pdf of $Z_{i,j}|W_{i,j}$ mentioned in the update step is either the pdf in (9) (i.e., inverse gamma) or the one in (15) (i.e., Gaussian), according to the option adopted to model the error term.

* Step 1 (prediction step), jth iteration:

for any $k = 1, ..., N$, compute the term $\sum_{k=1}^{j-1} w_{i,j} = \sum_{k=1}^{j-1} \Delta w_{i,j} + \sum_{k=1}^{j-1} w_{i,j-1}$, where $\sum_{k=1}^{j-1} w_{i,j-1}$ is the $k$th pseudorandom realization of $W(t_{i,j-1})$ given $Z_{i,j-1} = Z_{i,j-1}$ generated at the $(j-1)$th iteration and $\sum_{k=1}^{j-1} \Delta w_{i,j}$ is a pseudorandom realization of $\Delta W(t_{i,j-1}, t_{i,j})$. Then, append it to the particle vector $\sum_{k=1}^{j-1} w_{i,1}, ..., \sum_{k=1}^{j-1} w_{i,j-1}$ defined at the $(j-1)$th iteration. The output of this prediction step is a set of N vectors:

$$ \begin{gathered} \sum_{k=1}^{j-1} w_{i,1}, \dots, \sum_{k=1}^{j-1} w_{i,j-1}, \sum_{k=1}^{j-1} w_{i,j} \\ \vdots \\ \sum_{k=1}^{j-1} w_{i,1}, \dots, \sum_{k=1}^{j-1} w_{i,j-1}, \sum_{k=1}^{j-1} w_{i,j} \end{gathered} $$

that we will refer to as particles.

* Step 2 (update step), jth iteration:

for any $k = 1, ..., N$, compute the importance weight of the $k$th particle as:

$$ q_{i,j} = \frac{f_{Z_{i,j}|W_{i,j}} \left( Z_{i,j} \middle| \sum_{k=1}^{j-1} w_{i,j} \right)}{\sum_{k=1}^{N} f_{Z_{i,j}|W_{i,j}} \left( Z_{i,k} \middle| \sum_{k=1}^{j-1} w_{i,j} \right)} $$

hence, resample the set of particles:

$$ \begin{gathered} \sum_{k=1}^{j-1} w_{i,1}, \dots, \sum_{k=1}^{j-1} w_{i,j-1}, \sum_{k=1}^{j-1} w_{i,j} \\ \vdots \\ \sum_{k=1}^{j-1} w_{i,1}, \dots, \sum_{k=1}^{j-1} w_{i,j-1}, \sum_{k=1}^{j-1} w_{i,j} \end{gathered} $$

according to their importance weights and rename the new particles (i.e., the vectors) as:

$$ \begin{gathered} \sum_{k=1}^{j-1} w_{i,1}, \dots, \sum_{k=1}^{j-1} w_{i,j-1}, \sum_{k=1}^{j-1} w_{i,j} \\ \vdots \\ \sum_{k=1}^{j-1} w_{i,1}, \dots, \sum_{k=1}^{j-1} w_{i,j-1}, \sum_{k=1}^{j-1} w_{i,j} \end{gathered} $$

For $j = 1$, in the first prediction step, to initialize the algorithm, draw a pseudorandom sample of size N from $W(t_{i,1})$, denote its elements as $w_{i,1}, ..., w_{i,N}$, and define the particles as:

$$ \begin{gathered} \sum_{k=1}^{j-1} w_{i,1} = \sum_{k=1}^{j-1} w_{i,j} \\ \vdots \\ \sum_{k=1}^{j-1} w_{i,1}, ..., \sum_{k=1}^{j-1} w_{i,j} \end{gathered} $$
The particle $\stackrel{j-1}{k}w_{i,1}, \dots, \stackrel{j-1}{k}w_{i,j-1}, \stackrel{j-1}{k}w_{i,j}$ should be intended as a pseudorandom realization of $W_{i,j}$ given $Z_{i,j-1} = z_{i,j-1}$, and $\stackrel{j}{k}w_{i,1}, \dots, \stackrel{j}{k}w_{i,j-1}, \stackrel{j}{k}w_{i,j}$ should be intended as a pseudorandom realization of $W_{i,j}$ given $Z_{i,j} = z_{i,j}$, where $W_{i,j} = \{W(t_{i,1}), \dots, W(t_{i,j})\}, Z_{i,j} = \{Z(t_{i,1}), \dots, Z(t_{i,j})\}$, and $z_{i,j} = \{z(t_{i,1}), \dots, z(t_{i,j})\}$ is the realization of $Z_{i,j}$.

The conditional pdfs that are needed to compute the likelihood function (24) can be approximated as:

$$f_{Z_{i,j}|Z_{i,j-1}}(z_{i,j}|z_{i,j-1}) \approx \frac{\sum_{k=1}^{N} f_{Z_{i,j}|W_{i,j}}(z_{i,j}|\stackrel{j-1}{k}w_{i,j})}{N},$$

where $\stackrel{j-1}{k}w_{i,j}$ is the last component of the particle vector $\stackrel{j-1}{k}w_{i,1}, \dots, \stackrel{j-1}{k}w_{i,j-1}, \stackrel{j-1}{k}w_{i,j}$ generated at jth prediction step. Likewise, for example, the conditional mean of a function $g(W_i)$ of $W_i$, given $Z_i = z_i$ and $\xi$, can be computed as:

$$E\{g(W_i)|Z_i = z_i, \xi\} \approx \frac{\sum_{k=1}^{N} g(\stackrel{n_i}{k}w_i)}{N},$$

where $\stackrel{n_i}{k}w_i$ is the particle $\stackrel{n_i}{k}w_{i,1}, \dots, \stackrel{n_i}{k}w_{i,n_i-1}, \stackrel{n_i}{k}w_{i,n_i}$ generated at the n_i-th update step, under a perturbed model whose parameter vector is set to $\xi$. Obviously, the quality of these approximations improves with N.

This particle filter algorithm is also used to compute the ML estimate of the complementary cdf of $RUL(t) \bar{F}_{RUL(t)}(\tau|z(t))$ (22) and the ML estimate of the $MRUL(t)$ (23). In fact, more specifically, by using the notations introduced in Section 4, if $t_l \le t < t_{l+1}$, so that the set $z(t) = \{Z(t_j); j \ge 1, t_j \le t\}$ contains l perturbed measurements of the degradation level of a certain (selected) unit, given a pseudorandom sample of size N from $W(t_l)|Z(t) = z(t)$, say $\stackrel{l}{1}w_l, \dots, \stackrel{l}{N}w_l$, the ML estimate of the cdf $F_{RUL(t)}(\tau|z(t))$ (22) can be computed as follows:

$$F_{RUL(t)}(\tau|z(t)) \approx 1 - \frac{\sum_{k=1}^{N} F_{W(t+\tau)|W(t_l)}(\stackrel{l}{k}w_l)}{N},$$

where both the parameters of the perturbed model used to generate the particles and the parameters of the conditional cdf $F_{W(t+\tau)|W(t_l)}(\cdot|\cdot)$ are set to their ML estimates. Similarly, the ML estimate of the $MRUL(t)$ (23) is computed as follows:

$$MRUL(t) \approx \frac{\sum_{k=1}^{N} \int_{0}^{\infty} F_{W(t+\tau)|W(t_l)}(\stackrel{l}{k}w_l) \cdot d\tau}{N} = \int_{0}^{\infty} \frac{\sum_{k=1}^{N} F_{W(t+\tau)|W(t_l)}(\stackrel{l}{k}w_l) \cdot d\tau}{N},$$

where the integral is calculated numerically.